# Source code for pennylane.pulse.transmon

# Copyright 2018-2023 Xanadu Quantum Technologies Inc.

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""This module contains the classes/functions specific for simulation of superconducting transmon hardware systems"""
import warnings

from dataclasses import dataclass
from typing import Callable, List, Union

import pennylane as qml
import pennylane.numpy as np
from pennylane.pulse import HardwareHamiltonian
from pennylane.pulse.hardware_hamiltonian import HardwarePulse
from pennylane.typing import TensorLike
from pennylane.wires import Wires

# TODO ladder operators once there is qudit support
# pylint: disable=unused-argument
def a(wire, d=2):
"""creation operator"""
return qml.s_prod(0.5, qml.PauliX(wire)) + qml.s_prod(0.5j, qml.PauliY(wire))

"""annihilation operator"""
return qml.s_prod(0.5, qml.PauliX(wire)) + qml.s_prod(-0.5j, qml.PauliY(wire))

# pylint: disable=too-many-arguments
[docs]def transmon_interaction(
qubit_freq: Union[float, list],
connections: list,
coupling: Union[float, list],
wires: list,
anharmonicity=None,
d=2,
):
r"""Returns a :class:ParametrizedHamiltonian representing the circuit QED Hamiltonian of a
superconducting transmon system.

The Hamiltonian is given by

.. math::

H = \sum_{q\in \text{wires}} \omega_q b^\dagger_q b_q
+ \sum_{(i, j) \in \mathcal{C}} g_{ij} \left(b^\dagger_i b_j + b_j^\dagger b_i \right)
+ \sum_{q\in \text{wires}} \alpha_q b^\dagger_q b^\dagger_q b_q b_q

where :math:[b_p, b_q^\dagger] = \delta_{pq} are creation and annihilation operators.
The first term describes the effect of the dressed qubit frequencies qubit_freq :math:= \omega_q/ (2\pi),
the second term their coupling :math:= g_{ij}/(2\pi) and the last the
anharmonicity :math:= \alpha_q/(2\pi), which all can vary for
different qubits. In practice, these operators are restricted to a finite dimension of the
local Hilbert space (default d=2 corresponds to qubits).
In that case, the anharmonicity is set to :math:\alpha=0 and ignored.

The values of :math:\omega and :math:\alpha are typically around :math:5 \times 2\pi \text{GHz}
and :math:0.3 \times 2\pi \text{GHz}, respectively.
It is common for different qubits to be out of tune with different energy gaps. The coupling strength
:math:g typically varies between :math:[0.001, 0.1] \times 2\pi \text{GHz}. For some example parameters,
see e.g. arXiv:1804.04073 <https://arxiv.org/abs/1804.04073>_,
arXiv:2203.06818 <https://arxiv.org/abs/2203.06818>_, or arXiv:2210.15812 <https://arxiv.org/abs/2210.15812>_.

.. note:: Currently only supporting d=2 with qudit support planned in the future. For d=2, we have :math:b:=\frac{1}{2}(\sigma^x + i \sigma^y).

.. seealso::

:func:~.transmon_drive

Args:
qubit_freq (Union[float, list[float], Callable]): List of dressed qubit frequencies. This should be in units
of frequency (GHz), and will be converted to angular frequency :math:\omega internally where
needed, i.e. multiplied by :math:2 \pi. When passing a single float all qubits are assumed to
have that same frequency. When passing a parametrized function, it must have two
arguments, the first one being the trainable parameters and the second one being time.
connections (list[tuple(int)]): List of connections (i, j) between qubits i and j.
When the wires in connections are not contained in wires, a warning is raised.
coupling (Union[float, list[float]]): List of coupling strengths. This should be in units
of frequency (GHz), and will be converted to angular frequency internally where
needed, i.e. multiplied by :math:2 \pi. Needs to match the length of connections.
When passing a single float need explicit wires.
anharmonicity (Union[float, list[float]]): List of anharmonicities. This should be in units
of frequency (GHz), and will be converted to angular frequency internally where
needed, i.e. multiplied by :math:2 \pi. Ignored when d=2.
When passing a single float all qubits are assumed to have that same anharmonicity.
wires (list): Needs to be of the same length as qubit_freq. Note that there can be additional
wires in the resulting operator from the connections, which are treated independently.
d (int): Local Hilbert space dimension. Defaults to d=2 and is currently the only supported value.

Returns:
:class:~.ParametrizedHamiltonian: a :class:~.ParametrizedHamiltonian representing the transmon interaction

**Example**

We can set up the transmon interaction Hamiltonian with uniform coefficients by passing float values.

.. code-block::

connections = [[0, 1], [1, 3], [2, 1], [4, 5]]
H = qml.pulse.transmon_interaction(qubit_freq=0.5, connections=connections, coupling=1., wires=range(6))

The resulting :class:~.HardwareHamiltonian: consists of 4 coupling terms and 6 qubits
because there are six different wire indices in connections.

>>> print(H)
HardwareHamiltonian: terms=10

We can also provide individual values for each of the qubit energies and coupling strengths,
here of order :math:0.1 \times 2\pi\text{GHz} and :math:1 \times 2\pi\text{GHz}, respectively.

.. code-block::

qubit_freqs = [0.5, 0.4, 0.3, 0.2, 0.1, 0.]
couplings= [1., 2., 3., 4.]
H = qml.pulse.transmon_interaction(qubit_freq=qubit_freqs,
connections=connections,
coupling=couplings,
wires=range(6))

The interaction term is dependent only on the typically fixed transmon energies and coupling strengths.
Executing this as a pulse program via :func:~.evolve would correspond to all driving fields being turned off.
To add a driving field, see :func:~.transmon_drive.

"""
if d != 2:
raise NotImplementedError(
"Currently only supporting qubits. Qutrits and qudits are planned in the future."
)

# if wires is None and qml.math.ndim(omega) == 0:
#     raise ValueError(
#         f"Cannot instantiate wires automatically. Either need specific wires or a list of omega."
#         f"Received wires {wires} and omega of type {type(omega)}"
#     )

# wires = wires or list(range(len(omega)))

n_wires = len(wires)

if not Wires(wires).contains_wires(Wires(np.unique(connections).tolist())):
warnings.warn(
f"Caution, wires and connections do not match. "
f"I.e., wires in connections {connections} are not contained in the wires {wires}"
)

# Prepare coefficients
if anharmonicity is None:
anharmonicity = [0.0] * n_wires

# TODO: make coefficients callable / trainable. Currently not supported
if callable(qubit_freq) or qml.math.ndim(qubit_freq) == 0:
qubit_freq = [qubit_freq] * n_wires
elif len(qubit_freq) != n_wires:
raise ValueError(
f"Number of qubit frequencies in {qubit_freq} does not match the provided wires = {wires}"
)

if qml.math.ndim(coupling) == 0:
coupling = [coupling] * len(connections)
if len(coupling) != len(connections):
raise ValueError(
f"Number of coupling terms {coupling} does not match the provided connections = {connections}"
)

settings = TransmonSettings(connections, qubit_freq, coupling, anharmonicity=anharmonicity)

omega = [callable_freq_to_angular(f) if callable(f) else (2 * np.pi * f) for f in qubit_freq]
g = [callable_freq_to_angular(c) if callable(c) else (2 * np.pi * c) for c in coupling]

# qubit term
coeffs = list(omega)
observables = [ad(i, d) @ a(i, d) for i in wires]

# coupling term
coeffs += list(g)
observables += [ad(i, d) @ a(j, d) + ad(j, d) @ a(i, d) for (i, j) in connections]

# TODO Qudit support. Currently not supported but will be in the future.
# if d>2:
#     if qml.math.ndim(anharmonicity)==0:
#         anharmonicity = [anharmonicity] * n_wires
#     if len(anharmonicity) != n_wires:
#         raise ValueError(f"Number of qubit anharmonicities anharmonicity = {anharmonicity} does not match the provided wires = {wires}")
#     # anharmonicity term
#     alpha = [2 * np.pi * a for a in anharmonicity]
#     coeffs += list(alpha)
#     observables += [ad(i, d) @ ad(i, d) @ a(i, d) @ a(i, d) for i in wires]

return HardwareHamiltonian(
coeffs, observables, settings=settings, reorder_fn=_reorder_AmpPhaseFreq
)

def callable_freq_to_angular(fn):
"""Add a factor of 2pi to a callable result to convert from Hz to rad/s"""

def angular_fn(p, t):
return 2 * np.pi * fn(p, t)

return angular_fn

@dataclass
class TransmonSettings:
"""Dataclass that contains the information of a Transmon setup.

.. see-also:: :func:transmon_interaction

Args:
connections (List): List [[idx_q0, idx_q1], ..] of connected qubits (wires)
qubit_freq (List[float, Callable]):
coupling (List[list, TensorLike, Callable]):
anharmonicity (List[float, Callable]):

"""

connections: List
qubit_freq: Union[float, Callable]
coupling: Union[list, TensorLike, Callable]
anharmonicity: Union[float, Callable]

def __eq__(self, other):
return (
qml.math.all(self.connections == other.connections)
and qml.math.all(self.qubit_freq == other.qubit_freq)
and qml.math.all(self.coupling == other.coupling)
and qml.math.all(self.anharmonicity == other.anharmonicity)
)

if other is not None:
new_connections = list(self.connections) + list(other.connections)
new_qubit_freq = list(self.qubit_freq) + list(other.qubit_freq)
new_coupling = list(self.coupling) + list(other.coupling)
new_anh = list(self.anharmonicity) + list(other.anharmonicity)
return TransmonSettings(
new_connections, new_qubit_freq, new_coupling, anharmonicity=new_anh
)

return self

[docs]def transmon_drive(amplitude, phase, freq, wires, d=2):
r"""Returns a :class:ParametrizedHamiltonian representing the drive term of a transmon qubit.

The Hamiltonian is given by

.. math::

\Omega(t) \sin\left(\phi(t) + \nu t\right) \sum_q Y_q

where :math:\{Y_q\} are the Pauli-Y operators on wires :math:\{q\}.
The arguments amplitude, phase and freq correspond to :math:\Omega / (2\pi), :math:\phi
and :math:\nu / (2\pi), respectively, and can all be either fixed numbers (float) or depend on time
(callable). If they are time-dependent, they need to abide by the restrictions imposed
in :class:ParametrizedHamiltonian and have a signature of two parameters, (params, t).

Together with the qubit :math:Z terms in :func:transmon_interaction, driving with this term can generate
:math:X and :math:Y rotations by setting :math:\phi accordingly and driving on resonance
(see eqs. (79) - (92) in 1904.06560 <https://arxiv.org/abs/1904.06560>_).
Further, it can generate entangling gates by driving at cross-resonance with a coupled qubit
(see eqs. (131) - (137) in 1904.06560 <https://arxiv.org/abs/1904.06560>_).
Such a coupling is described in :func:transmon_interaction.

For realistic simulations, one may restrict the amplitude, phase and drive frequency parameters.
For example, the authors in 2008.04302 <https://arxiv.org/abs/2008.04302>_ impose the restrictions of
a maximum amplitude :math:\Omega_{\text{max}} = 20 \text{MHz} and the carrier frequency to deviate at most
:math:\nu - \omega = \pm 1 \text{GHz} from the qubit frequency :math:\omega
(see :func:~.transmon_interaction).
The phase :math:\phi(t) is typically a slowly changing function of time compared to :math:\Omega(t).

.. note:: Currently only supports d=2 with qudit support planned in the future.
For d>2, we have :math:Y \mapsto i (\sigma^- - \sigma^+)
with lowering and raising operators  :math:\sigma^{\mp}.

.. note:: Due to convention in the respective fields, we omit the factor :math:\frac{1}{2} present in the related constructor :func:~.rydberg_drive

.. seealso::

:func:~.rydberg_drive, :func:~.transmon_interaction

Args:
amplitude (Union[float, callable]): Float or callable representing the amplitude of the driving field.
This should be in units of frequency (GHz), and will be converted to angular frequency
:math:\Omega(t) internally where needed, i.e. multiplied by :math:2 \pi.
phase (Union[float, callable]): Float or callable returning phase :math:\phi(t) (in radians).
Can be a fixed number (float) or depend on time (callable)
freq (Union[float, callable]): Float or callable representing the frequency of the driving field.
This should be in units of frequency (GHz), and will be converted to angular frequency
:math:\nu(t) internally where needed, i.e. multiplied by :math:2 \pi.
wires (Union[int, list[int]]): Label of the qubit that the drive acts upon. Can be a list of multiple wires.
d (int): Local Hilbert space dimension. Defaults to d=2 and is currently the only supported value.

Returns:
:class:~.ParametrizedHamiltonian: a :class:~.ParametrizedHamiltonian representing the transmon interaction

**Example**

We can construct a drive term acting on qubit 0 in the following way. We parametrize the amplitude and phase
via :math:\Omega(t)/(2 \pi) = A \times \sin^2(\pi t) and :math:\phi(t) = \phi_0 (t - \frac{1}{2}). The squared
sine ensures that the amplitude will be strictly positive (a requirement for some hardware). For simplicity, we
set the drive frequency to zero :math:\nu=0.

.. code-block:: python3

def amp(A, t):
return A * jnp.exp(-t**2)

def phase(phi0, t):
return phi0

freq = 0

H = qml.pulse.transmon_drive(amp, phase, freq, 0)

t = 0.
A = 1.
phi0 = jnp.pi/2
params = [A, phi0]

Evaluated at :math:t = 0 with the parameters :math:A = 1 and :math:\phi_0 = \pi/2 we obtain
:math:2 \pi A \exp(0) \sin(\pi/2 + 0)\sigma^y = 2 \pi \sigma^y.

>>> H(params, t)
6.283185307179586*(PauliY(wires=))

We can combine transmon_drive() with :func:~.transmon_interaction to create a full driven transmon Hamiltonian.
Let us look at a chain of three transmon qubits that are coupled with their direct neighbors. We provide all
frequencies in GHz (conversion to angular frequency, i.e. multiplication by :math:2 \pi, is taken care of
internally where needed).

We use values around :math:\omega = 5 \times 2\pi \text{GHz} for resonant frequencies, and coupling strenghts
on the order of around :math:g = 0.01 \times 2\pi \text{GHz}.

We parametrize the drive Hamiltonians for the qubits with amplitudes as squared sinusodials of
maximum amplitude :math:A, and constant drive frequencies of value :math:\nu. We set the
phase to zero :math:\phi=0, and we make the parameters :math:A and :math:\nu trainable
for every qubit. We simulate the evolution for a time window of :math:[0, 5]\text{ns}.

.. code-block:: python3

qubit_freqs = [5.1, 5., 5.3]
connections = [[0, 1], [1, 2]]  # qubits 0 and 1 are coupled, as are 1 and 2
g = [0.02, 0.05]
H = qml.pulse.transmon_interaction(qubit_freqs, connections, g, wires=range(3))

def amp(max_amp, t): return max_amp * jnp.sin(t) ** 2
freq = qml.pulse.constant  # Parametrized constant frequency
phase = 0.0
time = 5

for q in range(3):
H += qml.pulse.transmon_drive(amp, phase, freq, q)  # Parametrized drive for each qubit

dev = qml.device("default.qubit.jax", wires=range(3))

@jax.jit
@qml.qnode(dev, interface="jax")
def qnode(params):
qml.evolve(H)(params, time)
return qml.expval(qml.PauliZ(0) + qml.PauliZ(1) + qml.PauliZ(2))

We evaluate the Hamiltonian with some arbitrarily chosen maximum amplitudes (here on the order of :math:0.5 \times 2\pi \text{GHz})
and set the drive frequency equal to the qubit frequencies. Note how the order of the construction
of H determines the order with which the parameters need to be passed to
:class:~.ParametrizedHamiltonian and :func:~.evolve. We made the drive frequencies
trainable parameters by providing constant callables through :func:~.pulse.constant instead of fixed values (like the phase).
This allows us to differentiate with respect to both the maximum amplitudes and the frequencies and optimize them.

>>> max_amp0, max_amp1, max_amp2 = [0.5, 0.3, 0.6]
>>> fr0, fr1, fr2 = qubit_freqs
>>> params = [max_amp0, fr0, max_amp1, fr1, max_amp2, fr2]
>>> qnode(params)
Array(-1.57851962, dtype=float64)
[Array(-13.50193649, dtype=float64),
Array(3.1112141, dtype=float64),
Array(16.40286521, dtype=float64),
Array(-4.30485667, dtype=float64),
Array(4.75813949, dtype=float64),
Array(3.43272354, dtype=float64)]

"""
if d != 2:
raise NotImplementedError(
"Currently only supports qubits (d=2). Qutrits and qudits support is planned in the future."
)

wires = Wires(wires)

# TODO: use creation and annihilation operators when introducing qutrits
# Note that exp(-iw)a* + exp(iw)a = cos(w)X - sin(w)Y for a=1/2(X+iY)
# We compute the coeffs and observables of the EM field
coeffs = [AmplitudeAndPhaseAndFreq(qml.math.sin, amplitude, phase, freq)]

drive_y_term = sum(qml.PauliY(wire) for wire in wires)

observables = [drive_y_term]

pulses = [HardwarePulse(amplitude, phase, freq, wires)]

return HardwareHamiltonian(coeffs, observables, pulses=pulses, reorder_fn=_reorder_AmpPhaseFreq)

# pylint:disable = too-few-public-methods,too-many-return-statements
class AmplitudeAndPhaseAndFreq:
"""Class storing combined amplitude, phase and freq callables"""

def __init__(self, trig_fn, amp, phase, freq, hz_to_rads=2 * np.pi):
self.amp_is_callable = callable(amp)
self.phase_is_callable = callable(phase)
self.freq_is_callable = callable(freq)

# all 3 callable

if self.amp_is_callable and self.phase_is_callable and self.freq_is_callable:

def callable_amp_and_phase_and_freq(params, t):
return (
* amp(params, t)
* trig_fn(phase(params, t) + hz_to_rads * freq(params, t) * t)
)

self.func = callable_amp_and_phase_and_freq
return

# 2 out of 3 callable

if self.amp_is_callable and self.phase_is_callable:

def callable_amp_and_phase(params, t):
return (
* amp(params, t)
* trig_fn(phase(params, t) + hz_to_rads * freq * t)
)

self.func = callable_amp_and_phase
return

if self.amp_is_callable and self.freq_is_callable:

def callable_amp_and_freq(params, t):
return (
* amp(params, t)
* trig_fn(phase + hz_to_rads * freq(params, t) * t)
)

self.func = callable_amp_and_freq
return

if self.phase_is_callable and self.freq_is_callable:

def callable_phase_and_freq(params, t):
return (
* amp
* trig_fn(phase(params, t) + hz_to_rads * freq(params, t) * t)
)

self.func = callable_phase_and_freq
return

# 1 out of 3 callable

if self.amp_is_callable:

def callable_amp(params, t):
return hz_to_rads * amp(params, t) * trig_fn(phase + hz_to_rads * freq * t)

self.func = callable_amp
return

if self.phase_is_callable:

def callable_phase(params, t):
return hz_to_rads * amp * trig_fn(phase(params, t) + hz_to_rads * freq * t)

self.func = callable_phase
return

if self.freq_is_callable:

def callable_freq(params, t):
return hz_to_rads * amp * trig_fn(phase + hz_to_rads * freq(params, t) * t)

self.func = callable_freq
return

# 0 out of 3 callable
# (the remaining coeff is still callable due to explicit time dependence)

def no_callable(_, t):
return hz_to_rads * amp * trig_fn(phase + hz_to_rads * freq * t)

self.func = no_callable

def __call__(self, params, t):
return self.func(params, t)

def _reorder_AmpPhaseFreq(params, coeffs_parametrized):
"""Takes params, and reorganizes it based on whether the Hamiltonian has
callable phase and/or callable amplitude and/or callable freq.

Consolidates amplitude, phase and freq parameters if they are callable,
and duplicates parameters since they will be passed to two operators in the Hamiltonian"""

reordered_params = []

coeff_idx = 0
params_idx = 0

for i, coeff in enumerate(coeffs_parametrized):
if i == coeff_idx:
if isinstance(coeff, AmplitudeAndPhaseAndFreq):
is_callables = [
coeff.phase_is_callable,
coeff.amp_is_callable,
coeff.freq_is_callable,
]

num_callables = sum(is_callables)

# package parameters according to how many coeffs are callable
reordered_params.extend([params[params_idx : params_idx + num_callables]])

coeff_idx += 1
params_idx += num_callables

else:
reordered_params.append(params[params_idx])
coeff_idx += 1
params_idx += 1

return reordered_params


Using PennyLane

Development

API

Internals