Cosine transmon model: from charge basis to strong-drive dynamics¶

This tutorial follows Qubex's model="cosine" path from the finite Cooper-pair charge basis to a truncated energy-basis model. It then compares the resulting Hamiltonian, interaction operator, strong-drive dynamics, and runtime with the Duffing approximation.

All calculations are offline. The saved numerical results are simulator results for the parameters below, not device measurements.

Goal¶

The notebook follows four steps:

1. construct the cosine Hamiltonian in the finite charge basis;
2. transform and truncate the Hamiltonian and charge operator in the energy basis;
3. compare the resulting cosine and Duffing models under a strong resonant drive; and
4. separate model-compilation cost from time-evolution cost.

The main comparison uses the same measured inputs, \(\omega_{01}/(2\pi)=5.0\,\mathrm{GHz}\) and \(\alpha/(2\pi)=(\omega_{12}-\omega_{01})/(2\pi)=-0.25\,\mathrm{GHz}\), for both local models.

In [ ]:

from time import perf_counter

import numpy as np
import plotly.graph_objects as go
from IPython.display import display
from numpy.testing import assert_allclose
from plotly.subplots import make_subplots
from qutip import Qobj

from qubex.simulator import (
    CompiledCosineTransmon,
    Control,
    QuantumSimulator,
    QuantumSystem,
    Transmon,
)

from time import perf_counter import numpy as np import plotly.graph_objects as go from IPython.display import display from numpy.testing import assert_allclose from plotly.subplots import make_subplots from qutip import Qobj from qubex.simulator import ( CompiledCosineTransmon, Control, QuantumSimulator, QuantumSystem, Transmon, )

1. Define matched local models¶

We use \(\hbar=1\) throughout, so Hamiltonians and the coefficients \(\omega\), \(\alpha\), and \(\Omega\) have angular-frequency units of rad/ns. The Transmon API accepts cyclic-frequency values in GHz: \(x\,\mathrm{GHz}\) corresponds to \(2\pi x\,\mathrm{rad/ns}\).

For model="cosine", Qubex fits positive charging and Josephson energies so that the finite-charge-basis spectrum reproduces the requested \(\omega_{01}\) and signed anharmonicity \(\alpha\). With \(\hbar=1\), \(E_C\) and \(E_J\) are expressed in angular-frequency units; the public charge_basis provenance reports \(E_C/(2\pi)\) and \(E_J/(2\pi)\) in GHz.

Both models below retain seven energy levels. The cosine model additionally uses a charge cutoff \(N_{\mathrm{c}}=15\), so its numerical charge basis contains \(2N_{\mathrm{c}}+1=31\) states before energy truncation.

In [ ]:

FREQUENCY = 5.0
ANHARMONICITY = -0.25
RETAINED_DIMENSION = 7
CHARGE_CUTOFF = 15
OFFSET_CHARGE = 0.0

duffing = Transmon(
    label="Q0",
    dimension=RETAINED_DIMENSION,
    frequency=FREQUENCY,
    anharmonicity=ANHARMONICITY,
    model="duffing",
    offset_charge=OFFSET_CHARGE,
)
cosine = Transmon(
    label="Q0",
    dimension=RETAINED_DIMENSION,
    frequency=FREQUENCY,
    anharmonicity=ANHARMONICITY,
    model="cosine",
    charge_cutoff=CHARGE_CUTOFF,
    offset_charge=OFFSET_CHARGE,
)

compiled_duffing = duffing.compile()
compiled_cosine = cosine.compile()
assert isinstance(compiled_cosine, CompiledCosineTransmon)
charge_basis = compiled_cosine.charge_basis

FREQUENCY = 5.0 ANHARMONICITY = -0.25 RETAINED_DIMENSION = 7 CHARGE_CUTOFF = 15 OFFSET_CHARGE = 0.0 duffing = Transmon( label="Q0", dimension=RETAINED_DIMENSION, frequency=FREQUENCY, anharmonicity=ANHARMONICITY, model="duffing", offset_charge=OFFSET_CHARGE, ) cosine = Transmon( label="Q0", dimension=RETAINED_DIMENSION, frequency=FREQUENCY, anharmonicity=ANHARMONICITY, model="cosine", charge_cutoff=CHARGE_CUTOFF, offset_charge=OFFSET_CHARGE, ) compiled_duffing = duffing.compile() compiled_cosine = cosine.compile() assert isinstance(compiled_cosine, CompiledCosineTransmon) charge_basis = compiled_cosine.charge_basis

2. Build the cosine Hamiltonian in the charge basis¶

Let \(\lvert n\rangle\) denote a Cooper-pair charge state with integer \(n\). With \(\hbar=1\), the angular-frequency Hamiltonian is

\[ \hat H_{\mathrm{charge}} =4E_C(\hat n-n_g)^2 -E_J\cos\hat\varphi. \]

The phase operator shifts charge, so

\[ \cos\hat\varphi =\frac{1}{2}\sum_n \left( \lvert n+1\rangle\langle n\rvert +\lvert n\rangle\langle n+1\rvert \right). \]

Therefore, the charging term is diagonal and the Josephson term is nearest-neighbor hopping with matrix element \(-E_J/2\). The finite range \(n=-N_{\mathrm{c}},\ldots,+N_{\mathrm{c}}\) is a numerical cutoff, not a physical bound on the number of Cooper pairs. The public charge_basis.hamiltonian is stored in cyclic GHz, so the code multiplies it by \(2\pi\) before checking this equation.

In [ ]:

charge_numbers = charge_basis.charge_numbers
relative_charge_values = charge_numbers - charge_basis.offset_charge
charging_energy = 2 * np.pi * charge_basis.charging_energy
josephson_energy = 2 * np.pi * charge_basis.josephson_energy
charge_hamiltonian = 2 * np.pi * charge_basis.hamiltonian
shift = np.diag(np.ones(charge_numbers.size - 1), k=1)
reconstructed_charge_hamiltonian = 4 * charging_energy * np.diag(
    relative_charge_values**2
) - 0.5 * josephson_energy * (shift + shift.T)

assert_allclose(
    reconstructed_charge_hamiltonian,
    charge_hamiltonian,
    rtol=0.0,
    atol=1e-12,
)

display(Qobj(charge_hamiltonian) / (2 * np.pi))

charge_numbers = charge_basis.charge_numbers relative_charge_values = charge_numbers - charge_basis.offset_charge charging_energy = 2 * np.pi * charge_basis.charging_energy josephson_energy = 2 * np.pi * charge_basis.josephson_energy charge_hamiltonian = 2 * np.pi * charge_basis.hamiltonian shift = np.diag(np.ones(charge_numbers.size - 1), k=1) reconstructed_charge_hamiltonian = 4 * charging_energy * np.diag( relative_charge_values**2 ) - 0.5 * josephson_energy * (shift + shift.T) assert_allclose( reconstructed_charge_hamiltonian, charge_hamiltonian, rtol=0.0, atol=1e-12, ) display(Qobj(charge_hamiltonian) / (2 * np.pi))

Transform first, then truncate¶

Qubex first diagonalizes the finite charge-basis Hamiltonian and then retains its lowest \(d\) eigenstates. Collect these eigenvectors as the columns of \(\hat V^{(d)}\):

\[ \hat H_{\mathrm{charge}}\hat V^{(d)} =\hat V^{(d)}\operatorname{diag}(\omega_0,\ldots,\omega_{d-1}), \]

After shifting the ground-state energy to zero, the retained energy-basis Hamiltonian is

\[ \hat H_{\cos}^{(d)} =\operatorname{diag} \left( \omega_0-\omega_0,\ldots, \omega_{d-1}-\omega_0 \right). \]

The relative charge operator must be projected with the same \(\hat V^{(d)}\):

\[ \hat n_{\mathrm{proj}}^{(d)} =\hat V^{(d)\dagger}(\hat n-n_g)\hat V^{(d)}, \qquad \hat{\widetilde n}^{(d)} =\frac{\hat n_{\mathrm{proj}}^{(d)}} {\left\lvert\langle0\rvert\hat n_{\mathrm{proj}}^{(d)}\lvert1\rangle\right\rvert}. \]

Qubex exposes \(\hat{\widetilde n}^{(d)}\) as interaction_operator. Its normalization keeps the weak-drive \(0\leftrightarrow1\) rate and adjacent exchange-strength convention aligned with the Duffing model. It is dimensionless; it is not the unnormalized physical charge \(2e(\hat n-n_g)\).

There are two distinct numerical truncations:

• Charge cutoff \(N_{\mathrm{c}}\): truncate the infinite charge ladder to \(2N_{\mathrm{c}}+1\) states for diagonalization.
• Retained dimension \(d\): after diagonalization, keep the lowest \(d\) energy eigenstates and project every operator into that subspace.

In [ ]:

eigenvectors = charge_basis.eigenvectors
eigenenergies = 2 * np.pi * charge_basis.eigenenergies
projected_hamiltonian = eigenvectors.conj().T @ charge_hamiltonian @ eigenvectors
relative_charge_operator = np.diag(relative_charge_values)
projected_charge = eigenvectors.conj().T @ relative_charge_operator @ eigenvectors
projected_charge = 0.5 * (projected_charge + projected_charge.conj().T)
normalized_projected_charge = projected_charge / abs(projected_charge[0, 1])

compiled_cosine_angular_energies = np.asarray(
    compiled_cosine.hamiltonian.diag(), dtype=float
)
assert_allclose(
    projected_hamiltonian - eigenenergies[0] * np.eye(RETAINED_DIMENSION),
    np.diag(compiled_cosine_angular_energies),
    rtol=0.0,
    atol=2e-12,
)
assert_allclose(
    normalized_projected_charge,
    compiled_cosine.interaction_operator.full(),
    rtol=0.0,
    atol=2e-13,
)

eigenvectors = charge_basis.eigenvectors eigenenergies = 2 * np.pi * charge_basis.eigenenergies projected_hamiltonian = eigenvectors.conj().T @ charge_hamiltonian @ eigenvectors relative_charge_operator = np.diag(relative_charge_values) projected_charge = eigenvectors.conj().T @ relative_charge_operator @ eigenvectors projected_charge = 0.5 * (projected_charge + projected_charge.conj().T) normalized_projected_charge = projected_charge / abs(projected_charge[0, 1]) compiled_cosine_angular_energies = np.asarray( compiled_cosine.hamiltonian.diag(), dtype=float ) assert_allclose( projected_hamiltonian - eigenenergies[0] * np.eye(RETAINED_DIMENSION), np.diag(compiled_cosine_angular_energies), rtol=0.0, atol=2e-12, ) assert_allclose( normalized_projected_charge, compiled_cosine.interaction_operator.full(), rtol=0.0, atol=2e-13, )

Energy eigenvectors in the charge basis¶

The signed amplitudes \(\langle n\vert m\rangle\) provide a direct check of the charge cutoff: the retained wavefunctions should decay before \(n=\pm N_{\mathrm{c}}\). In the low-energy, nearly harmonic regime, each successive state also gains one node. The overall sign of an eigenvector is arbitrary.

In [ ]:

charge_basis_amplitudes = np.real(eigenvectors.T).copy()
for level in range(RETAINED_DIMENSION):
    largest_component = np.argmax(np.abs(charge_basis_amplitudes[level]))
    if charge_basis_amplitudes[level, largest_component] < 0:
        charge_basis_amplitudes[level] *= -1

amplitude_figure = make_subplots(
    rows=7,
    cols=1,
    shared_xaxes=True,
    shared_yaxes=True,
    subplot_titles=tuple(
        f"Energy level m = {level}" for level in range(RETAINED_DIMENSION)
    ),
    vertical_spacing=0.045,
)
for level in range(RETAINED_DIMENSION):
    row = level + 1
    amplitude_figure.add_trace(
        go.Scatter(
            x=charge_numbers,
            y=charge_basis_amplitudes[level],
            mode="lines+markers",
            line={"color": "#4C78A8", "width": 2.0},
            marker={"size": 4},
            showlegend=False,
            hovertemplate=(f"m={level}<br>n=%{{x}}<br>⟨n|m⟩=%{{y:.5f}}<extra></extra>"),
        ),
        row=row,
        col=1,
    )
amplitude_figure.update_xaxes(range=[-CHARGE_CUTOFF, CHARGE_CUTOFF], matches="x")
amplitude_figure.update_xaxes(
    title_text="Charge number n",
    row=7,
    col=1,
)
amplitude_figure.update_yaxes(
    title_text="Amplitude",
    range=[-0.65, 0.65],
    zeroline=True,
    zerolinecolor="#6B7280",
)
amplitude_figure.update_layout(
    title={
        "text": "Low-energy eigenvectors in the charge basis<br><sup>Signed amplitudes ⟨n|m⟩; one additional node per level</sup>",
        "x": 0.5,
    },
    template="plotly_white",
    font={"family": "Arial, sans-serif", "color": "#252A34"},
    margin={"l": 70, "r": 35, "t": 85, "b": 65},
    width=1000,
    height=1200,
)
display(amplitude_figure)

charge_basis_amplitudes = np.real(eigenvectors.T).copy() for level in range(RETAINED_DIMENSION): largest_component = np.argmax(np.abs(charge_basis_amplitudes[level])) if charge_basis_amplitudes[level, largest_component] < 0: charge_basis_amplitudes[level] *= -1 amplitude_figure = make_subplots( rows=7, cols=1, shared_xaxes=True, shared_yaxes=True, subplot_titles=tuple( f"Energy level m = {level}" for level in range(RETAINED_DIMENSION) ), vertical_spacing=0.045, ) for level in range(RETAINED_DIMENSION): row = level + 1 amplitude_figure.add_trace( go.Scatter( x=charge_numbers, y=charge_basis_amplitudes[level], mode="lines+markers", line={"color": "#4C78A8", "width": 2.0}, marker={"size": 4}, showlegend=False, hovertemplate=(f"m={level}<br>n=%{{x}}<br>⟨n|m⟩=%{{y:.5f}}<extra></extra>"), ), row=row, col=1, ) amplitude_figure.update_xaxes(range=[-CHARGE_CUTOFF, CHARGE_CUTOFF], matches="x") amplitude_figure.update_xaxes( title_text="Charge number n", row=7, col=1, ) amplitude_figure.update_yaxes( title_text="Amplitude", range=[-0.65, 0.65], zeroline=True, zerolinecolor="#6B7280", ) amplitude_figure.update_layout( title={ "text": "Low-energy eigenvectors in the charge basis<br>^{Signed amplitudes ⟨n|m⟩; one additional node per level}", "x": 0.5, }, template="plotly_white", font={"family": "Arial, sans-serif", "color": "#252A34"}, margin={"l": 70, "r": 35, "t": 85, "b": 65}, width=1000, height=1200, ) display(amplitude_figure)

3. Compare the retained Hamiltonians and operators¶

The Duffing model uses an oscillator basis and the analytic Hamiltonian

\[ \hat H_{\mathrm{Duffing}} =\omega_{01}\hat N+\frac{\alpha}{2}\hat N(\hat N-1), \qquad \hat X_{\mathrm{D}}=\hat a+\hat a^\dagger. \]

The cosine fit matches \(\omega_{01}\) and \(\alpha\), so levels 0, 1, and 2 agree within the fit tolerance. The models first separate at level 3, where a single constant Kerr coefficient no longer reproduces the cosine spectrum. The adjacent interaction matrix elements likewise depart from the Duffing values \(\sqrt{m+1}\).

In [ ]:

duffing_angular_energies = np.asarray(compiled_duffing.hamiltonian.diag(), dtype=float)
cosine_angular_energies = np.asarray(compiled_cosine.hamiltonian.diag(), dtype=float)
duffing_transition_frequencies = np.diff(duffing_angular_energies)
cosine_transition_frequencies = np.diff(cosine_angular_energies)
duffing_adjacent_elements = np.abs(
    np.diag(compiled_duffing.interaction_operator.full(), k=1)
)
cosine_adjacent_elements = np.abs(
    np.diag(compiled_cosine.interaction_operator.full(), k=1)
)
assert_allclose(duffing_adjacent_elements, np.sqrt(np.arange(1, RETAINED_DIMENSION)))

duffing_angular_energies = np.asarray(compiled_duffing.hamiltonian.diag(), dtype=float) cosine_angular_energies = np.asarray(compiled_cosine.hamiltonian.diag(), dtype=float) duffing_transition_frequencies = np.diff(duffing_angular_energies) cosine_transition_frequencies = np.diff(cosine_angular_energies) duffing_adjacent_elements = np.abs( np.diag(compiled_duffing.interaction_operator.full(), k=1) ) cosine_adjacent_elements = np.abs( np.diag(compiled_cosine.interaction_operator.full(), k=1) ) assert_allclose(duffing_adjacent_elements, np.sqrt(np.arange(1, RETAINED_DIMENSION)))

In [ ]:

comparison_figure = make_subplots(
    rows=2,
    cols=2,
    subplot_titles=(
        "Adjacent transition angular frequencies",
        "Retained angular-frequency difference",
        "Adjacent interaction-operator magnitudes",
        "Adjacent-element difference",
    ),
    horizontal_spacing=0.14,
    vertical_spacing=0.16,
)
transition_labels = [f"{level}→{level + 1}" for level in range(RETAINED_DIMENSION - 1)]
comparison_figure.add_trace(
    go.Scatter(
        x=transition_labels,
        y=duffing_transition_frequencies / (2 * np.pi),
        mode="lines+markers",
        name="Duffing",
        line={"color": "#4C78A8", "width": 2.5},
        marker={"size": 7},
    ),
    row=1,
    col=1,
)
comparison_figure.add_trace(
    go.Scatter(
        x=transition_labels,
        y=cosine_transition_frequencies / (2 * np.pi),
        mode="lines+markers",
        name="Cosine",
        line={"color": "#F58518", "dash": "dash", "width": 2.5},
        marker={"size": 7},
    ),
    row=1,
    col=1,
)
comparison_figure.add_trace(
    go.Bar(
        x=[str(level) for level in range(RETAINED_DIMENSION)],
        y=1e3 * (cosine_angular_energies - duffing_angular_energies) / (2 * np.pi),
        name="Cosine - Duffing",
        marker={"color": "#6B7280"},
        showlegend=False,
    ),
    row=1,
    col=2,
)
comparison_figure.add_trace(
    go.Scatter(
        x=transition_labels,
        y=duffing_adjacent_elements,
        mode="lines+markers",
        name="Duffing",
        line={"color": "#4C78A8", "width": 2.5},
        marker={"size": 7},
        showlegend=False,
    ),
    row=2,
    col=1,
)
comparison_figure.add_trace(
    go.Scatter(
        x=transition_labels,
        y=cosine_adjacent_elements,
        mode="lines+markers",
        name="Cosine",
        line={"color": "#F58518", "dash": "dash", "width": 2.5},
        marker={"size": 7},
        showlegend=False,
    ),
    row=2,
    col=1,
)
comparison_figure.add_trace(
    go.Bar(
        x=transition_labels,
        y=cosine_adjacent_elements - duffing_adjacent_elements,
        name="Cosine - Duffing",
        marker={"color": "#6B7280"},
        showlegend=False,
    ),
    row=2,
    col=2,
)
comparison_figure.update_xaxes(title_text="Transition", row=1, col=1)
comparison_figure.update_yaxes(title_text="ω / 2π (GHz)", row=1, col=1)
comparison_figure.update_xaxes(title_text="Energy level m", row=1, col=2)
comparison_figure.update_yaxes(
    title_text="Angular-frequency difference / 2π (MHz)",
    zeroline=True,
    zerolinecolor="#252A34",
    row=1,
    col=2,
)
comparison_figure.update_xaxes(title_text="Transition", row=2, col=1)
comparison_figure.update_yaxes(
    title_text="Normalized magnitude",
    rangemode="tozero",
    row=2,
    col=1,
)
comparison_figure.update_xaxes(title_text="Transition", row=2, col=2)
comparison_figure.update_yaxes(
    title_text="Cosine - Duffing",
    zeroline=True,
    zerolinecolor="#252A34",
    row=2,
    col=2,
)
comparison_figure.update_layout(
    title={
        "text": "Cosine and Duffing retained spectra and interaction elements<br><sup>Matched ω₀₁ and α; adjacent elements normalized at 0↔1; d = 7</sup>",
        "x": 0.5,
    },
    legend={"orientation": "h", "y": 1.13, "x": 0.0},
)
comparison_figure.update_layout(
    template="plotly_white",
    font={"family": "Arial, sans-serif", "color": "#252A34"},
    margin={"l": 70, "r": 35, "t": 85, "b": 65},
    width=1000,
    height=800,
)
display(comparison_figure)

comparison_figure = make_subplots( rows=2, cols=2, subplot_titles=( "Adjacent transition angular frequencies", "Retained angular-frequency difference", "Adjacent interaction-operator magnitudes", "Adjacent-element difference", ), horizontal_spacing=0.14, vertical_spacing=0.16, ) transition_labels = [f"{level}→{level + 1}" for level in range(RETAINED_DIMENSION - 1)] comparison_figure.add_trace( go.Scatter( x=transition_labels, y=duffing_transition_frequencies / (2 * np.pi), mode="lines+markers", name="Duffing", line={"color": "#4C78A8", "width": 2.5}, marker={"size": 7}, ), row=1, col=1, ) comparison_figure.add_trace( go.Scatter( x=transition_labels, y=cosine_transition_frequencies / (2 * np.pi), mode="lines+markers", name="Cosine", line={"color": "#F58518", "dash": "dash", "width": 2.5}, marker={"size": 7}, ), row=1, col=1, ) comparison_figure.add_trace( go.Bar( x=[str(level) for level in range(RETAINED_DIMENSION)], y=1e3 * (cosine_angular_energies - duffing_angular_energies) / (2 * np.pi), name="Cosine - Duffing", marker={"color": "#6B7280"}, showlegend=False, ), row=1, col=2, ) comparison_figure.add_trace( go.Scatter( x=transition_labels, y=duffing_adjacent_elements, mode="lines+markers", name="Duffing", line={"color": "#4C78A8", "width": 2.5}, marker={"size": 7}, showlegend=False, ), row=2, col=1, ) comparison_figure.add_trace( go.Scatter( x=transition_labels, y=cosine_adjacent_elements, mode="lines+markers", name="Cosine", line={"color": "#F58518", "dash": "dash", "width": 2.5}, marker={"size": 7}, showlegend=False, ), row=2, col=1, ) comparison_figure.add_trace( go.Bar( x=transition_labels, y=cosine_adjacent_elements - duffing_adjacent_elements, name="Cosine - Duffing", marker={"color": "#6B7280"}, showlegend=False, ), row=2, col=2, ) comparison_figure.update_xaxes(title_text="Transition", row=1, col=1) comparison_figure.update_yaxes(title_text="ω / 2π (GHz)", row=1, col=1) comparison_figure.update_xaxes(title_text="Energy level m", row=1, col=2) comparison_figure.update_yaxes( title_text="Angular-frequency difference / 2π (MHz)", zeroline=True, zerolinecolor="#252A34", row=1, col=2, ) comparison_figure.update_xaxes(title_text="Transition", row=2, col=1) comparison_figure.update_yaxes( title_text="Normalized magnitude", rangemode="tozero", row=2, col=1, ) comparison_figure.update_xaxes(title_text="Transition", row=2, col=2) comparison_figure.update_yaxes( title_text="Cosine - Duffing", zeroline=True, zerolinecolor="#252A34", row=2, col=2, ) comparison_figure.update_layout( title={ "text": "Cosine and Duffing retained spectra and interaction elements<br>^{Matched ω₀₁ and α; adjacent elements normalized at 0↔1; d = 7}", "x": 0.5, }, legend={"orientation": "h", "y": 1.13, "x": 0.0}, ) comparison_figure.update_layout( template="plotly_white", font={"family": "Arial, sans-serif", "color": "#252A34"}, margin={"l": 70, "r": 35, "t": 85, "b": 65}, width=1000, height=800, ) display(comparison_figure)

In [ ]:

duffing_interaction_magnitude = np.abs(compiled_duffing.interaction_operator.full())
cosine_interaction_magnitude = np.abs(compiled_cosine.interaction_operator.full())
interaction_scale_max = max(
    np.max(duffing_interaction_magnitude),
    np.max(cosine_interaction_magnitude),
)
interaction_figure = make_subplots(
    rows=1,
    cols=2,
    shared_yaxes=True,
    subplot_titles=(
        "Duffing interaction operator",
        "Cosine projected charge operator",
    ),
    horizontal_spacing=0.12,
)
interaction_figure.add_trace(
    go.Heatmap(
        z=duffing_interaction_magnitude,
        coloraxis="coloraxis",
        hovertemplate="row=%{y}<br>column=%{x}<br>|X|=%{z:.6f}<extra></extra>",
    ),
    row=1,
    col=1,
)
interaction_figure.add_trace(
    go.Heatmap(
        z=cosine_interaction_magnitude,
        coloraxis="coloraxis",
        hovertemplate="row=%{y}<br>column=%{x}<br>|ñ|=%{z:.6f}<extra></extra>",
    ),
    row=1,
    col=2,
)
interaction_figure.update_xaxes(title_text="Energy-basis column")
interaction_figure.update_yaxes(
    title_text="Energy-basis row",
    autorange="reversed",
)
interaction_figure.update_layout(
    title={
        "text": "Interaction operators after energy truncation<br><sup>Absolute matrix elements; common color scale</sup>",
        "x": 0.5,
    },
    coloraxis={
        "colorscale": "Blues",
        "cmin": 0.0,
        "cmax": interaction_scale_max,
        "colorbar": {"title": "Magnitude"},
    },
    template="plotly_white",
    font={"family": "Arial, sans-serif", "color": "#252A34"},
    margin={"l": 70, "r": 80, "t": 85, "b": 65},
    width=1000,
    height=500,
)
display(interaction_figure)

duffing_interaction_magnitude = np.abs(compiled_duffing.interaction_operator.full()) cosine_interaction_magnitude = np.abs(compiled_cosine.interaction_operator.full()) interaction_scale_max = max( np.max(duffing_interaction_magnitude), np.max(cosine_interaction_magnitude), ) interaction_figure = make_subplots( rows=1, cols=2, shared_yaxes=True, subplot_titles=( "Duffing interaction operator", "Cosine projected charge operator", ), horizontal_spacing=0.12, ) interaction_figure.add_trace( go.Heatmap( z=duffing_interaction_magnitude, coloraxis="coloraxis", hovertemplate="row=%{y}<br>column=%{x}<br>|X|=%{z:.6f}<extra></extra>", ), row=1, col=1, ) interaction_figure.add_trace( go.Heatmap( z=cosine_interaction_magnitude, coloraxis="coloraxis", hovertemplate="row=%{y}<br>column=%{x}<br>|ñ|=%{z:.6f}<extra></extra>", ), row=1, col=2, ) interaction_figure.update_xaxes(title_text="Energy-basis column") interaction_figure.update_yaxes( title_text="Energy-basis row", autorange="reversed", ) interaction_figure.update_layout( title={ "text": "Interaction operators after energy truncation<br>^{Absolute matrix elements; common color scale}", "x": 0.5, }, coloraxis={ "colorscale": "Blues", "cmin": 0.0, "cmax": interaction_scale_max, "colorbar": {"title": "Magnitude"}, }, template="plotly_white", font={"family": "Arial, sans-serif", "color": "#252A34"}, margin={"l": 70, "r": 80, "t": 85, "b": 65}, width=1000, height=500, ) display(interaction_figure)

The cosine interaction_operator retains the full projected-charge matrix, including diagonal and nonadjacent elements. For driven evolution, current Qubex keeps only its first upper off-diagonal and constructs

\[ \hat L=\sum_{m=0}^{d-2}\widetilde n_{m,m+1} \lvert m\rangle\langle m+1\rvert, \]

and uses \(\hat L\) and \(\hat L^\dagger\) in the co-rotating drive Hamiltonian. This is the drive rotating-wave approximation used by QuantumSimulator.

Qubex currently assumes that the error from this drive RWA is smaller than the error introduced by replacing the cosine model with the Duffing approximation, and therefore neglects the omitted terms. The comparison below tests the local-model difference under that assumption; it does not quantify the RWA error itself.

4. Compare idealized strong-drive dynamics¶

We apply the same resonant rectangular drive to both retained models: \(\Omega/(2\pi)=600\,\mathrm{MHz}\), \(\omega_{\mathrm{d}}=\omega_{01}\), and a duration of \(20\,\mathrm{ns}\). This deliberately strong, idealized drive reveals higher-level differences; it is not an AWG-ready pulse.

In the drive rotating frame, the simulated single-transmon Hamiltonian is

\[ \hat H_{\mathrm{sim}}(t) =\hat H_0-\omega_{\mathrm{d}}\hat N +\frac{1}{2} \left[ \Omega(t)\hat L^\dagger+\Omega^*(t)\hat L \right]. \]

Comparison boundary:

• local Hamiltonian: Duffing or full cosine spectrum, both truncated to \(d=7\);
• coupling RWA: not applicable because there is one uncoupled transmon;
• drive RWA: applied; only adjacent elements of \(\hat L\) are retained;
• frame: drive frame rotating at \(\omega_{\mathrm{d}}\);
• decoherence: omitted; and
• amplitude convention: the normalized \(0\leftrightarrow1\) matrix element is one for both models, so the same \(\Omega\) gives the same weak-drive Rabi rate, not necessarily the same laboratory AWG voltage.

Consequently, this is not a full-charge, lab-frame, counter-rotating strong-drive validation.

In [ ]:

DRIVE_OMEGA_OVER_2PI_GHZ = 0.60
DRIVE_OMEGA = 2 * np.pi * DRIVE_OMEGA_OVER_2PI_GHZ
DRIVE_FREQUENCY = FREQUENCY
DRIVE_DURATION_NS = 20.0
OUTPUT_SAMPLES = 401
SOLVER_OPTIONS = {
    "rtol": 1e-9,
    "atol": 1e-11,
    "max_step": 0.02,
}

duffing_system = QuantumSystem(objects=[duffing])
cosine_system = QuantumSystem(objects=[cosine])
duffing_simulator = QuantumSimulator(duffing_system)
cosine_simulator = QuantumSimulator(cosine_system)
duffing_control = Control(
    target=duffing,
    waveform=np.array([DRIVE_OMEGA], dtype=np.complex128),
    durations=np.array([DRIVE_DURATION_NS]),
    frequency=DRIVE_FREQUENCY,
)
cosine_control = Control(
    target=cosine,
    waveform=np.array([DRIVE_OMEGA], dtype=np.complex128),
    durations=np.array([DRIVE_DURATION_NS]),
    frequency=DRIVE_FREQUENCY,
)

duffing_result = duffing_simulator.sesolve(
    [duffing_control],
    initial_state=duffing_system.ground_state,
    n_samples=OUTPUT_SAMPLES,
    options=SOLVER_OPTIONS,
)
cosine_result = cosine_simulator.sesolve(
    [cosine_control],
    initial_state=cosine_system.ground_state,
    n_samples=OUTPUT_SAMPLES,
    options=SOLVER_OPTIONS,
)
duffing_populations = np.array(
    [np.abs(state.full().ravel()) ** 2 for state in duffing_result.states]
)
cosine_populations = np.array(
    [np.abs(state.full().ravel()) ** 2 for state in cosine_result.states]
)
times_ns = duffing_result.times

DRIVE_OMEGA_OVER_2PI_GHZ = 0.60 DRIVE_OMEGA = 2 * np.pi * DRIVE_OMEGA_OVER_2PI_GHZ DRIVE_FREQUENCY = FREQUENCY DRIVE_DURATION_NS = 20.0 OUTPUT_SAMPLES = 401 SOLVER_OPTIONS = { "rtol": 1e-9, "atol": 1e-11, "max_step": 0.02, } duffing_system = QuantumSystem(objects=[duffing]) cosine_system = QuantumSystem(objects=[cosine]) duffing_simulator = QuantumSimulator(duffing_system) cosine_simulator = QuantumSimulator(cosine_system) duffing_control = Control( target=duffing, waveform=np.array([DRIVE_OMEGA], dtype=np.complex128), durations=np.array([DRIVE_DURATION_NS]), frequency=DRIVE_FREQUENCY, ) cosine_control = Control( target=cosine, waveform=np.array([DRIVE_OMEGA], dtype=np.complex128), durations=np.array([DRIVE_DURATION_NS]), frequency=DRIVE_FREQUENCY, ) duffing_result = duffing_simulator.sesolve( [duffing_control], initial_state=duffing_system.ground_state, n_samples=OUTPUT_SAMPLES, options=SOLVER_OPTIONS, ) cosine_result = cosine_simulator.sesolve( [cosine_control], initial_state=cosine_system.ground_state, n_samples=OUTPUT_SAMPLES, options=SOLVER_OPTIONS, ) duffing_populations = np.array( [np.abs(state.full().ravel()) ** 2 for state in duffing_result.states] ) cosine_populations = np.array( [np.abs(state.full().ravel()) ** 2 for state in cosine_result.states] ) times_ns = duffing_result.times

In [ ]:

population_figure = make_subplots(
    rows=7,
    cols=1,
    shared_xaxes=True,
    shared_yaxes=True,
    subplot_titles=tuple(
        f"Energy level m = {level}" for level in range(RETAINED_DIMENSION)
    ),
    vertical_spacing=0.045,
)
for level in range(RETAINED_DIMENSION):
    row = level + 1
    population_figure.add_trace(
        go.Scatter(
            x=times_ns,
            y=duffing_populations[:, level],
            mode="lines",
            name="Duffing",
            legendgroup="duffing",
            showlegend=row == 1,
            line={"color": "#4C78A8", "width": 2.2},
        ),
        row=row,
        col=1,
    )
    population_figure.add_trace(
        go.Scatter(
            x=times_ns,
            y=cosine_populations[:, level],
            mode="lines",
            name="Cosine",
            legendgroup="cosine",
            showlegend=row == 1,
            line={"color": "#F58518", "dash": "dash", "width": 2.2},
        ),
        row=row,
        col=1,
    )
    population_figure.update_yaxes(
        title_text="Population",
        range=[0.0, 1.0],
        row=row,
        col=1,
    )
population_figure.update_xaxes(title_text="Time (ns)", row=7, col=1)
population_figure.update_layout(
    title={
        "text": "Resonant-drive energy-level populations<br><sup>Ω / (2π) = 600 MHz, 20 ns, d = 7, drive RWA</sup>",
        "x": 0.5,
    },
    legend={"orientation": "h", "y": 1.08, "x": 0.0},
)
population_figure.update_layout(
    template="plotly_white",
    font={"family": "Arial, sans-serif", "color": "#252A34"},
    margin={"l": 70, "r": 35, "t": 85, "b": 65},
    width=1000,
    height=1400,
)
display(population_figure)

population_figure = make_subplots( rows=7, cols=1, shared_xaxes=True, shared_yaxes=True, subplot_titles=tuple( f"Energy level m = {level}" for level in range(RETAINED_DIMENSION) ), vertical_spacing=0.045, ) for level in range(RETAINED_DIMENSION): row = level + 1 population_figure.add_trace( go.Scatter( x=times_ns, y=duffing_populations[:, level], mode="lines", name="Duffing", legendgroup="duffing", showlegend=row == 1, line={"color": "#4C78A8", "width": 2.2}, ), row=row, col=1, ) population_figure.add_trace( go.Scatter( x=times_ns, y=cosine_populations[:, level], mode="lines", name="Cosine", legendgroup="cosine", showlegend=row == 1, line={"color": "#F58518", "dash": "dash", "width": 2.2}, ), row=row, col=1, ) population_figure.update_yaxes( title_text="Population", range=[0.0, 1.0], row=row, col=1, ) population_figure.update_xaxes(title_text="Time (ns)", row=7, col=1) population_figure.update_layout( title={ "text": "Resonant-drive energy-level populations<br>^{Ω / (2π) = 600 MHz, 20 ns, d = 7, drive RWA}", "x": 0.5, }, legend={"orientation": "h", "y": 1.08, "x": 0.0}, ) population_figure.update_layout( template="plotly_white", font={"family": "Arial, sans-serif", "color": "#252A34"}, margin={"l": 70, "r": 35, "t": 85, "b": 65}, width=1000, height=1400, ) display(population_figure)

The two models initially follow the same motion because their \(0\leftrightarrow1\) transition and weak-drive normalization are matched. At this drive strength, however, levels \(m\geq2\) acquire substantial population. Differences in the higher transition frequencies and adjacent charge matrix elements then accumulate as coherent phase shifts and redistribute population across the ladder. The alternating separation of the curves is therefore a multilevel interference effect, not a monotonic increase or decrease of leakage in one model.

5. Compare execution speed¶

Separate the runtime into two stages:

1. Model compilation: Duffing constructs small analytic matrices directly. Cosine fitting repeatedly diagonalizes a \(31\times31\) charge Hamiltonian, then performs the final projection.
2. Precompiled evolution: QuantumSystem compiles each object once during initialization. Both models then evolve matrices with the same retained dimension \(d=7\) and the same solver settings.

The benchmark alternates model order, excludes the first warm-up call, and reports the median and interquartile range. Absolute values depend on the machine and software environment; the useful comparison is the separation between compilation and propagation. A lab-frame non-RWA simulation would be substantially slower because the integrator must resolve the carrier oscillation with much finer time steps. That cost is outside this drive-RWA benchmark.

In [ ]:

duffing.compile()
cosine.compile()
duffing_compile_samples_ms = np.empty(15)
cosine_compile_samples_ms = np.empty(15)
for repeat in range(15):
    order = ("duffing", "cosine") if repeat % 2 == 0 else ("cosine", "duffing")
    for model in order:
        start = perf_counter()
        Transmon(
            label="Q0",
            dimension=RETAINED_DIMENSION,
            frequency=FREQUENCY,
            anharmonicity=ANHARMONICITY,
            model=model,
            charge_cutoff=CHARGE_CUTOFF if model == "cosine" else None,
            offset_charge=OFFSET_CHARGE,
        ).compile()
        elapsed_ms = 1e3 * (perf_counter() - start)
        if model == "duffing":
            duffing_compile_samples_ms[repeat] = elapsed_ms
        else:
            cosine_compile_samples_ms[repeat] = elapsed_ms

duffing_simulator.sesolve(
    [duffing_control],
    initial_state=duffing_system.ground_state,
    n_samples=OUTPUT_SAMPLES,
    options=SOLVER_OPTIONS,
)
cosine_simulator.sesolve(
    [cosine_control],
    initial_state=cosine_system.ground_state,
    n_samples=OUTPUT_SAMPLES,
    options=SOLVER_OPTIONS,
)
duffing_evolution_samples_ms = np.empty(9)
cosine_evolution_samples_ms = np.empty(9)
for repeat in range(9):
    order = ("duffing", "cosine") if repeat % 2 == 0 else ("cosine", "duffing")
    for model in order:
        start = perf_counter()
        if model == "duffing":
            duffing_simulator.sesolve(
                [duffing_control],
                initial_state=duffing_system.ground_state,
                n_samples=OUTPUT_SAMPLES,
                options=SOLVER_OPTIONS,
            )
            duffing_evolution_samples_ms[repeat] = 1e3 * (perf_counter() - start)
        else:
            cosine_simulator.sesolve(
                [cosine_control],
                initial_state=cosine_system.ground_state,
                n_samples=OUTPUT_SAMPLES,
                options=SOLVER_OPTIONS,
            )
            cosine_evolution_samples_ms[repeat] = 1e3 * (perf_counter() - start)

duffing.compile() cosine.compile() duffing_compile_samples_ms = np.empty(15) cosine_compile_samples_ms = np.empty(15) for repeat in range(15): order = ("duffing", "cosine") if repeat % 2 == 0 else ("cosine", "duffing") for model in order: start = perf_counter() Transmon( label="Q0", dimension=RETAINED_DIMENSION, frequency=FREQUENCY, anharmonicity=ANHARMONICITY, model=model, charge_cutoff=CHARGE_CUTOFF if model == "cosine" else None, offset_charge=OFFSET_CHARGE, ).compile() elapsed_ms = 1e3 * (perf_counter() - start) if model == "duffing": duffing_compile_samples_ms[repeat] = elapsed_ms else: cosine_compile_samples_ms[repeat] = elapsed_ms duffing_simulator.sesolve( [duffing_control], initial_state=duffing_system.ground_state, n_samples=OUTPUT_SAMPLES, options=SOLVER_OPTIONS, ) cosine_simulator.sesolve( [cosine_control], initial_state=cosine_system.ground_state, n_samples=OUTPUT_SAMPLES, options=SOLVER_OPTIONS, ) duffing_evolution_samples_ms = np.empty(9) cosine_evolution_samples_ms = np.empty(9) for repeat in range(9): order = ("duffing", "cosine") if repeat % 2 == 0 else ("cosine", "duffing") for model in order: start = perf_counter() if model == "duffing": duffing_simulator.sesolve( [duffing_control], initial_state=duffing_system.ground_state, n_samples=OUTPUT_SAMPLES, options=SOLVER_OPTIONS, ) duffing_evolution_samples_ms[repeat] = 1e3 * (perf_counter() - start) else: cosine_simulator.sesolve( [cosine_control], initial_state=cosine_system.ground_state, n_samples=OUTPUT_SAMPLES, options=SOLVER_OPTIONS, ) cosine_evolution_samples_ms[repeat] = 1e3 * (perf_counter() - start)

In [ ]:

timing_figure = make_subplots(
    rows=1,
    cols=2,
    shared_yaxes=True,
    subplot_titles=("Model compilation", "Precompiled time evolution"),
    horizontal_spacing=0.16,
)
for column, duffing_samples, cosine_samples in (
    (1, duffing_compile_samples_ms, cosine_compile_samples_ms),
    (2, duffing_evolution_samples_ms, cosine_evolution_samples_ms),
):
    medians = [np.median(duffing_samples), np.median(cosine_samples)]
    q1_values = [np.quantile(duffing_samples, 0.25), np.quantile(cosine_samples, 0.25)]
    q3_values = [np.quantile(duffing_samples, 0.75), np.quantile(cosine_samples, 0.75)]
    timing_figure.add_trace(
        go.Bar(
            x=["Duffing", "Cosine"],
            y=medians,
            marker={"color": ["#4C78A8", "#F58518"]},
            text=[f"{value:.3f} ms" for value in medians],
            textposition="outside",
            cliponaxis=False,
            error_y={
                "type": "data",
                "array": [
                    q3 - median for q3, median in zip(q3_values, medians, strict=True)
                ],
                "arrayminus": [
                    median - q1 for median, q1 in zip(medians, q1_values, strict=True)
                ],
                "color": "#252A34",
            },
            showlegend=False,
        ),
        row=1,
        col=column,
    )
    timing_figure.update_xaxes(title_text="Local model", row=1, col=column)
    timing_figure.update_yaxes(rangemode="tozero", row=1, col=column)
timing_figure.update_yaxes(title_text="Median runtime (ms)", row=1, col=1)
timing_figure.update_layout(
    title={
        "text": "Runtime comparison by execution phase<br><sup>Warm-up excluded; bars show median and IQR on this machine</sup>",
        "x": 0.5,
    }
)
timing_figure.update_layout(
    template="plotly_white",
    font={"family": "Arial, sans-serif", "color": "#252A34"},
    margin={"l": 70, "r": 35, "t": 85, "b": 65},
    width=1000,
    height=450,
)
display(timing_figure)

timing_figure = make_subplots( rows=1, cols=2, shared_yaxes=True, subplot_titles=("Model compilation", "Precompiled time evolution"), horizontal_spacing=0.16, ) for column, duffing_samples, cosine_samples in ( (1, duffing_compile_samples_ms, cosine_compile_samples_ms), (2, duffing_evolution_samples_ms, cosine_evolution_samples_ms), ): medians = [np.median(duffing_samples), np.median(cosine_samples)] q1_values = [np.quantile(duffing_samples, 0.25), np.quantile(cosine_samples, 0.25)] q3_values = [np.quantile(duffing_samples, 0.75), np.quantile(cosine_samples, 0.75)] timing_figure.add_trace( go.Bar( x=["Duffing", "Cosine"], y=medians, marker={"color": ["#4C78A8", "#F58518"]}, text=[f"{value:.3f} ms" for value in medians], textposition="outside", cliponaxis=False, error_y={ "type": "data", "array": [ q3 - median for q3, median in zip(q3_values, medians, strict=True) ], "arrayminus": [ median - q1 for median, q1 in zip(medians, q1_values, strict=True) ], "color": "#252A34", }, showlegend=False, ), row=1, col=column, ) timing_figure.update_xaxes(title_text="Local model", row=1, col=column) timing_figure.update_yaxes(rangemode="tozero", row=1, col=column) timing_figure.update_yaxes(title_text="Median runtime (ms)", row=1, col=1) timing_figure.update_layout( title={ "text": "Runtime comparison by execution phase<br>^{Warm-up excluded; bars show median and IQR on this machine}", "x": 0.5, } ) timing_figure.update_layout( template="plotly_white", font={"family": "Arial, sans-serif", "color": "#252A34"}, margin={"l": 70, "r": 35, "t": 85, "b": 65}, width=1000, height=450, ) display(timing_figure)

Takeaways¶

• Qubex fits \(E_C\) and \(E_J\), diagonalizes the finite charge-basis cosine Hamiltonian, and only then retains the lowest \(d\) energy eigenstates.
• The same retained eigenvectors must transform the relative charge operator. Qubex normalizes the projected operator by its \(0\leftrightarrow1\) matrix element to preserve the existing weak-drive and exchange-strength convention.
• Matching \(\omega_{01}\) and \(\alpha\) makes levels 0–2 agree with the Duffing model. The first spectral difference appears at level 3, together with increasing differences in adjacent charge matrix elements.
• A sufficiently strong drive populates these higher levels, where the spectral and matrix-element differences accumulate as coherent phase shifts and change the level populations.
• The displayed dynamics isolate this local-model difference under drive RWA in the \(\omega_{01}\) rotating frame. They do not validate a full lab-frame drive or a laboratory amplitude.
• Cosine compilation is slower because it includes parameter fitting and repeated charge-basis diagonalization. After compilation, equal retained dimensions make the propagation cost much closer.
• Choose \(N_{\mathrm{c}}\) so the retained charge-basis amplitudes are negligible at \(n=\pm N_{\mathrm{c}}\), and increase \(d\) until the driven observable is converged.