# Source code for scqubits.core.transmon

# transmon.py
#
# This file is part of scqubits: a Python package for superconducting qubits,
# arXiv:2107.08552 (2021). https://arxiv.org/abs/2107.08552
#
#    Copyright (c) 2019 and later, Jens Koch and Peter Groszkowski
#
#    LICENSE file in the root directory of this source tree.
############################################################################

import math
import os

from typing import Any, Dict, List, Optional, Tuple, Union

import numpy as np

from matplotlib.axes import Axes
from matplotlib.figure import Figure
from numpy import ndarray

import scqubits.core.constants as constants
import scqubits.core.descriptors as descriptors
import scqubits.core.discretization as discretization
import scqubits.core.qubit_base as base
import scqubits.core.storage as storage
import scqubits.io_utils.fileio_serializers as serializers
import scqubits.utils.plot_defaults as defaults
import scqubits.utils.plotting as plot

from scqubits.core.discretization import Grid1d
from scqubits.core.noise import NoisySystem
from scqubits.core.storage import WaveFunction

LevelsTuple = Tuple[int, ...]
Transition = Tuple[int, int]
TransitionsTuple = Tuple[Transition, ...]

# —Cooper pair box / transmon——————————————————————————————————————————————

[docs]class Transmon(base.QubitBaseClass1d, serializers.Serializable, NoisySystem): r"""Class for the Cooper-pair-box and transmon qubit. The Hamiltonian is represented in dense form in the number basis, :math:H_\text{CPB}=4E_\text{C}(\hat{n}-n_g)^2+\frac{E_\text{J}}{2}(|n\rangle\langle n+1|+\text{h.c.}). Initialize with, for example:: Transmon(EJ=1.0, EC=2.0, ng=0.2, ncut=30) Parameters ---------- EJ: Josephson energy EC: charging energy ng: offset charge ncut: charge basis cutoff, n = -ncut, ..., ncut truncated_dim: desired dimension of the truncated quantum system; expected: truncated_dim > 1 id_str: optional string by which this instance can be referred to in HilbertSpace and ParameterSweep. If not provided, an id is auto-generated. """ EJ = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE") EC = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE") ng = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE") ncut = descriptors.WatchedProperty(int, "QUANTUMSYSTEM_UPDATE") def __init__( self, EJ: float, EC: float, ng: float, ncut: int, truncated_dim: int = 6, id_str: Optional[str] = None, ) -> None: base.QuantumSystem.__init__(self, id_str=id_str) self.EJ = EJ self.EC = EC self.ng = ng self.ncut = ncut self.truncated_dim = truncated_dim self._evec_dtype = np.float_ self._default_grid = discretization.Grid1d(-np.pi, np.pi, 151) self._default_n_range = (-5, 6) self._image_filename = os.path.join( os.path.dirname(os.path.abspath(__file__)), "qubit_img/fixed-transmon.jpg" )
[docs] @staticmethod def default_params() -> Dict[str, Any]: return {"EJ": 15.0, "EC": 0.3, "ng": 0.0, "ncut": 30, "truncated_dim": 10}
[docs] @classmethod def supported_noise_channels(cls) -> List[str]: """Return a list of supported noise channels""" return [ "tphi_1_over_f_cc", "tphi_1_over_f_ng", "t1_capacitive", "t1_charge_impedance", ]
[docs] @classmethod def effective_noise_channels(cls) -> List[str]: """Return a default list of channels used when calculating effective t1 and t2 noise.""" noise_channels = cls.supported_noise_channels() noise_channels.remove("t1_charge_impedance") return noise_channels
[docs] def n_operator(self) -> ndarray: """Returns charge operator n in the charge basis""" diag_elements = np.arange(-self.ncut, self.ncut + 1, 1) return np.diag(diag_elements)
[docs] def exp_i_phi_operator(self) -> ndarray: """Returns operator :math:e^{i\\varphi} in the charge basis""" dimension = self.hilbertdim() entries = np.repeat(1.0, dimension - 1) exp_op = np.diag(entries, -1) return exp_op
[docs] def cos_phi_operator(self) -> ndarray: """Returns operator :math:\\cos \\varphi in the charge basis""" cos_op = 0.5 * self.exp_i_phi_operator() cos_op += cos_op.T return cos_op
[docs] def sin_phi_operator(self) -> ndarray: """Returns operator :math:\\sin \\varphi in the charge basis""" sin_op = -1j * 0.5 * self.exp_i_phi_operator() sin_op += sin_op.conjugate().T return sin_op
[docs] def hamiltonian(self) -> ndarray: """Returns Hamiltonian in charge basis""" dimension = self.hilbertdim() hamiltonian_mat = np.diag( [ 4.0 * self.EC * (ind - self.ncut - self.ng) ** 2 for ind in range(dimension) ] ) ind = np.arange(dimension - 1) hamiltonian_mat[ind, ind + 1] = -self.EJ / 2.0 hamiltonian_mat[ind + 1, ind] = -self.EJ / 2.0 return hamiltonian_mat
[docs] def d_hamiltonian_d_ng(self) -> ndarray: """Returns operator representing a derivative of the Hamiltonian with respect to charge offset ng.""" return -8 * self.EC * self.n_operator()
[docs] def d_hamiltonian_d_EJ(self) -> ndarray: """Returns operator representing a derivative of the Hamiltonian with respect to EJ.""" return -self.cos_phi_operator()
[docs] def hilbertdim(self) -> int: """Returns Hilbert space dimension""" return 2 * self.ncut + 1
[docs] def potential(self, phi: Union[float, ndarray]) -> ndarray: """Transmon phase-basis potential evaluated at phi. Parameters ---------- phi: phase variable value """ return -self.EJ * np.cos(phi)
[docs] def plot_n_wavefunction( self, esys: Tuple[ndarray, ndarray] = None, mode: str = "real", which: int = 0, nrange: Tuple[int, int] = None, **kwargs ) -> Tuple[Figure, Axes]: """Plots transmon wave function in charge basis Parameters ---------- esys: eigenvalues, eigenvectors mode: 'abs_sqr', 'abs', 'real', 'imag' which: index or indices of wave functions to plot (default value = 0) nrange: range of n to be included on the x-axis (default value = (-5,6)) **kwargs: plotting parameters """ if nrange is None: nrange = self._default_n_range n_wavefunc = self.numberbasis_wavefunction(esys, which=which) amplitude_modifier = constants.MODE_FUNC_DICT[mode] n_wavefunc.amplitudes = amplitude_modifier(n_wavefunc.amplitudes) kwargs = { **defaults.wavefunction1d_discrete(mode), **kwargs, } # if any duplicates, later ones survive return plot.wavefunction1d_discrete(n_wavefunc, xlim=nrange, **kwargs)
[docs] def plot_phi_wavefunction( self, esys: Tuple[ndarray, ndarray] = None, which: int = 0, phi_grid: Grid1d = None, mode: str = "abs_sqr", scaling: float = None, **kwargs ) -> Tuple[Figure, Axes]: """Alias for plot_wavefunction""" return self.plot_wavefunction( esys=esys, which=which, phi_grid=phi_grid, mode=mode, scaling=scaling, **kwargs )
[docs] def numberbasis_wavefunction( self, esys: Tuple[ndarray, ndarray] = None, which: int = 0 ) -> WaveFunction: """Return the transmon wave function in number basis. The specific index of the wave function to be returned is which. Parameters ---------- esys: if None, the eigensystem is calculated on the fly; otherwise, the provided eigenvalue, eigenvector arrays as obtained from .eigensystem(), are used (default value = None) which: eigenfunction index (default value = 0) """ if esys is None: evals_count = max(which + 1, 3) esys = self.eigensys(evals_count=evals_count) evals, evecs = esys n_vals = np.arange(-self.ncut, self.ncut + 1) return storage.WaveFunction(n_vals, evecs[:, which], evals[which])
[docs] def wavefunction( self, esys: Optional[Tuple[ndarray, ndarray]] = None, which: int = 0, phi_grid: Grid1d = None, ) -> WaveFunction: """Return the transmon wave function in phase basis. The specific index of the wavefunction is which. esys can be provided, but if set to None then it is calculated on the fly. Parameters ---------- esys: if None, the eigensystem is calculated on the fly; otherwise, the provided eigenvalue, eigenvector arrays as obtained from .eigensystem() are used which: eigenfunction index (default value = 0) phi_grid: used for setting a custom grid for phi; if None use self._default_grid """ if esys is None: evals_count = max(which + 1, 3) evals, evecs = self.eigensys(evals_count=evals_count) else: evals, evecs = esys n_wavefunc = self.numberbasis_wavefunction(esys, which=which) phi_grid = phi_grid or self._default_grid phi_basis_labels = phi_grid.make_linspace() phi_wavefunc_amplitudes = np.empty(phi_grid.pt_count, dtype=np.complex_) for k in range(phi_grid.pt_count): phi_wavefunc_amplitudes[k] = (1j ** which / math.sqrt(2 * np.pi)) * np.sum( n_wavefunc.amplitudes * np.exp(1j * phi_basis_labels[k] * n_wavefunc.basis_labels) ) return storage.WaveFunction( basis_labels=phi_basis_labels, amplitudes=phi_wavefunc_amplitudes, energy=evals[which], )
def _compute_dispersion( self, dispersion_name: str, param_name: str, param_vals: ndarray, transitions_tuple: TransitionsTuple = ((0, 1),), levels_tuple: Optional[LevelsTuple] = None, point_count: int = 50, num_cpus: Optional[int] = None, ) -> Tuple[ndarray, ndarray]: if dispersion_name != "ng": return super()._compute_dispersion( dispersion_name, param_name, param_vals, transitions_tuple=transitions_tuple, levels_tuple=levels_tuple, point_count=point_count, num_cpus=num_cpus, ) max_level = ( np.max(transitions_tuple) if levels_tuple is None else np.max(levels_tuple) ) previous_ng = self.ng self.ng = 0.0 specdata_ng_0 = self.get_spectrum_vs_paramvals( param_name, param_vals, evals_count=max_level + 1, get_eigenstates=False, num_cpus=num_cpus, ) self.ng = 0.5 specdata_ng_05 = self.get_spectrum_vs_paramvals( param_name, param_vals, evals_count=max_level + 1, get_eigenstates=False, num_cpus=num_cpus, ) self.ng = previous_ng if levels_tuple is not None: dispersion = np.asarray( [ [ np.abs( specdata_ng_0.energy_table[param_index, j] - specdata_ng_05.energy_table[param_index, j] ) for param_index, _ in enumerate(param_vals) ] for j in levels_tuple ] ) return specdata_ng_0.energy_table, dispersion dispersion_list = [] for i, j in transitions_tuple: list_ij = [] for param_index, _ in enumerate(param_vals): ei_0 = specdata_ng_0.energy_table[param_index, i] ei_05 = specdata_ng_05.energy_table[param_index, i] ej_0 = specdata_ng_0.energy_table[param_index, j] ej_05 = specdata_ng_05.energy_table[param_index, j] list_ij.append( np.max([np.abs(ei_0 - ej_0), np.abs(ei_05 - ej_05)]) - np.min([np.abs(ei_0 - ej_0), np.abs(ei_05 - ej_05)]) ) dispersion_list.append(list_ij) return specdata_ng_0.energy_table, np.asarray(dispersion_list)
# — Flux-tunable Cooper pair box / transmon———————————————————————————————————————————
[docs]class TunableTransmon(Transmon, serializers.Serializable, NoisySystem): r"""Class for the flux-tunable transmon qubit. The Hamiltonian is represented in dense form in the number basis, :math:H_\text{CPB}=4E_\text{C}(\hat{ n}-n_g)^2+\frac{\mathcal{E}_\text{J}(\Phi)}{2}(|n\rangle\langle n+1|+\text{ h.c.}), Here, the effective Josephson energy is flux-tunable: :math:\mathcal{ E}_J(\Phi) = E_{J,\text{max}} \sqrt{\cos^2(\pi\Phi/\Phi_0) + d^2 \sin^2( \pi\Phi/\Phi_0)} and :math:d=(E_{J2}-E_{J1})(E_{J1}+E_{J2}) parametrizes the junction asymmetry. Initialize with, for example:: TunableTransmon(EJmax=1.0, d=0.1, EC=2.0, flux=0.3, ng=0.2, ncut=30) Parameters ---------- EJmax: maximum effective Josephson energy (sum of the Josephson energies of the two junctions) d: junction asymmetry parameter EC: charging energy flux: flux threading the SQUID loop, in units of the flux quantum ng: offset charge ncut: charge basis cutoff, n = -ncut, ..., ncut truncated_dim: desired dimension of the truncated quantum system; expected: truncated_dim > 1 id_str: optional string by which this instance can be referred to in HilbertSpace and ParameterSweep. If not provided, an id is auto-generated. """ EJmax = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE") d = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE") flux = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE") def __init__( self, EJmax: float, EC: float, d: float, flux: float, ng: float, ncut: int, truncated_dim: int = 6, id_str: Optional[str] = None, ) -> None: base.QuantumSystem.__init__(self, id_str=id_str) self.EJmax = EJmax self.EC = EC self.d = d self.flux = flux self.ng = ng self.ncut = ncut self.truncated_dim = truncated_dim self._evec_dtype = np.float_ self._default_grid = discretization.Grid1d(-np.pi, np.pi, 151) self._default_n_range = (-5, 6) self._image_filename = os.path.join( os.path.dirname(os.path.abspath(__file__)), "qubit_img/tunable-transmon.jpg" ) @property def EJ(self) -> float: # type: ignore """This is the effective, flux dependent Josephson energy, playing the role of EJ in the parent class Transmon""" return self.EJmax * np.sqrt( np.cos(np.pi * self.flux) ** 2 + self.d ** 2 * np.sin(np.pi * self.flux) ** 2 )
[docs] @staticmethod def default_params() -> Dict[str, Any]: return { "EJmax": 20.0, "EC": 0.3, "d": 0.01, "flux": 0.0, "ng": 0.0, "ncut": 30, "truncated_dim": 10, }
[docs] @classmethod def supported_noise_channels(cls) -> List[str]: """Return a list of supported noise channels""" return [ "tphi_1_over_f_flux", "tphi_1_over_f_cc", "tphi_1_over_f_ng", "t1_capacitive", "t1_flux_bias_line", "t1_charge_impedance", ]
[docs] def d_hamiltonian_d_flux(self) -> ndarray: """Returns operator representing a derivative of the Hamiltonian with respect to flux.""" return ( np.pi * self.EJmax * np.cos(np.pi * self.flux) * np.sin(np.pi * self.flux) * (self.d ** 2 - 1) / np.sqrt( np.cos(np.pi * self.flux) ** 2 + self.d ** 2 * np.sin(np.pi * self.flux) ** 2 ) * self.cos_phi_operator() )
def _compute_dispersion( self, dispersion_name: str, param_name: str, param_vals: ndarray, transitions_tuple: TransitionsTuple = ((0, 1),), levels_tuple: Optional[LevelsTuple] = None, point_count: int = 50, num_cpus: Optional[int] = None, ) -> Tuple[ndarray, ndarray]: if dispersion_name != "flux": return super()._compute_dispersion( dispersion_name, param_name, param_vals, transitions_tuple=transitions_tuple, levels_tuple=levels_tuple, point_count=point_count, num_cpus=num_cpus, ) max_level = ( np.max(transitions_tuple) if levels_tuple is None else np.max(levels_tuple) ) previous_flux = self.flux self.flux = 0.0 specdata_flux_0 = self.get_spectrum_vs_paramvals( param_name, param_vals, evals_count=max_level + 1, get_eigenstates=False, num_cpus=num_cpus, ) self.flux = 0.5 specdata_flux_05 = self.get_spectrum_vs_paramvals( param_name, param_vals, evals_count=max_level + 1, get_eigenstates=False, num_cpus=num_cpus, ) self.flux = previous_flux if levels_tuple is not None: dispersion = np.asarray( [ [ np.abs( specdata_flux_0.energy_table[param_index, j] # type:ignore - specdata_flux_05.energy_table[ param_index, j ] # type:ignore ) for param_index, _ in enumerate(param_vals) ] for j in levels_tuple ] ) return specdata_flux_0.energy_table, dispersion # type:ignore dispersion_list = [] for i, j in transitions_tuple: list_ij = [] for param_index, _ in enumerate(param_vals): ei_0 = specdata_flux_0.energy_table[param_index, i] # type:ignore ei_05 = specdata_flux_05.energy_table[param_index, i] # type:ignore ej_0 = specdata_flux_0.energy_table[param_index, j] # type:ignore ej_05 = specdata_flux_05.energy_table[param_index, j] # type:ignore list_ij.append( np.max([np.abs(ei_0 - ej_0), np.abs(ei_05 - ej_05)]) - np.min([np.abs(ei_0 - ej_0), np.abs(ei_05 - ej_05)]) ) dispersion_list.append(list_ij) return specdata_flux_0.energy_table, np.asarray(dispersion_list) # type:ignore