Scqubits: a Python package for superconducting qubits

$\textbf{scqubits}$ is an open-source Python package for simulating and analyzing superconducting circuits. It provides convenient routines to obtain energy spectra of common superconducting qubits, such as the transmon, fluxonium, flux, cos(2$\phi$) and the 0-$\pi$ qubit. $\textbf{scqubits}$ also features a number of options for visualizing the computed spectral data, including plots of energy levels as a function of external parameters, display of matrix elements of various operators as well as means to easily plot qubit wavefunctions. Many of these tools are not limited to single qubits, but extend to composite Hilbert spaces consisting of coupled superconducting qubits and harmonic (or weakly anharmonic) modes. The library provides an extensive suite of methods for estimating qubit coherence times due to a variety of commonly considered noise channels. While all functionality of $\textbf{scqubits}$ can be accessed programatically, the package also implements GUI-like widgets that, with a few clicks can help users both create relevant Python objects, as well as explore their properties through various plots. When applicable, the library harnesses the computing power of multiple cores via multiprocessing. $\textbf{scqubits}$ further exposes a direct interface to the Quantum Toolbox in Python (QuTiP) package, allowing the user to efficiently leverage QuTiP's proven capabilities for simulating time evolution.


Introduction
Superconducting qubits [7,10,17,27] have secured the rank of one of the most promising and widely researched hardware architectures for quantum information processing. All devices in this category are relatively simple circuits, and display genuine quantum properties such as discrete energy spectra and quantum-coherent time evolution. Coupling superconducting qubits to external electromagnetic fields enables the realization of gate operations as well as quantum-state readout using the framework of circuit quantum electrodynamics (circuit QED) [3,4,26].
The computation of energy spectra, eigenstates, and matrix elements of relevant operators is a key prerequisite for the design and fabrication of superconducting qubits, as well as for the quantitative analysis of experimental data collected in state-ofthe-art experiments. Circuit quantization [6,9,25] provides the systematic framework for deriving the Hamiltonian operator which describes a given circuit mathematically. However, with the exception of the simple LC oscillator and the limiting behavior of nonlinear circuits in specific parameter regimes, computation of qubit spectra cannot be accomplished analytically, but rather requires numerical solution of Hermitian eigenvalue problems.
The scqubits package provides a user-friendly, object-oriented Python library of the most common superconducting qubits. It facilitates automatic construction of circuit Hamiltonians in an appropriate basis, provides high-level routines for finding eigenenergies, eigenstates, and matrix elements, and allows the user to quickly visualize these quantities as a function of external parameters. While the scqubits routines for numerics and plotting rely heavily on NumPy, SciPy, and Matplotlib, no detailed knowledge of these internals is required from the user to employ scqubits efficiently. In this way, the package helps lower initial barriers encountered by new researchers entering the field, and can be easily integrated for educational purposes. At the same time, the library aims to fulfill the role of a unifying tool for expert workers in the field, enabling quick and efficient investigations of a wide variety of circuit-QED systems, and the interactive exploration of their behavior in different parameter regimes.
The analysis of interesting circuit-QED systems invariably involves the consideration of coupled systems, composed of qubits on one hand, and harmonic modes on the other hand (realized, for example, as on-chip transmission-line resonators or 3d cavities). The scqubits library simplifies the setup of such composite Hilbert spaces and grants a seamless interface to the well-established QuTiP [14,15] framework which can be leveraged for the simulation of time evolution.
This paper is organized as follows. In the next section, we provide an overview of the library and its core functionality. In Sec. 3 we discuss how to construct and analyze composite Hilbert space systems that may consist of multiple qubits and/or resonators. Next, in Sec. 4 we present how scqubits lets users perform parameter sweeps and easily explore how properties of a composite system vary as circuit parameters or control fields change. In Sec. 5 we review how scqubits can be used to estimate coherence times of different qubits, and how those coherence times can be visualized. After that, in Sec. 6 we provide a brief overview of the interactive exploration capabilities of the library, that allow users to study properties of various systems in ways that require very little actual programming. Finally, we summarize and conclude in Sec. 8.

Overview of the scqubits library
This section gives a broad overview and introduces the main building blocks of scqubits based on concrete examples of their usage. We stress that a jupyter notebook containing all of the source code in this manuscript can be found in the github examples repository (see Sec. 7). As a regular Python package, scqubits is imported via 1 import scqubits as scq The simplest way to explore spectral properties of individual superconducting qubits is to invoke the dedicated graphical user interface via 1 scq.GUI () which outputs the display shown in Fig. 1. With the exception of the initial call, usage of this interface requires no further knowledge of Python and is particularly well-suited for beginners. In the following we describe the more powerful and flexible programmatic usage of scqubits.
As an illustration of basic scqubits functionality, we consider the transmon qubit. Like all qubit types implemented in scqubits, Transmon and its flux-tunable variant TunableTransmon are realized as Python classes. Each class instance stores all relevant circuit parameters and control-field values, and provides a collection of methods used for common computations and visualization. (See Table 3 for a summary of the most commonly used qubit class methods.) An instance of the TunableTransmon class is created by the following call which provides all necessary system parameters for initialization: 1 tmon = scq.TunableTransmon( The initialization arguments include the relevant circuit parameters: for a flux-tunable transmon, these are the maximum Josephson energy E Jmax from the SQUID loop, charging energy E C , and offset charge n g (further details are provided in appendix A.III). If scqubits is used inside a jupyter notebook, then 1 transmon = scq.TunableTransmon.create() is an alternative way to create and initialize a TunableTransmon instance. The resulting graphical interface offers simple widgets to enter required parameters and displays the transmon circuit for reference. By default, energies are assumed to be   given as frequencies in units of GHz, although this can be easily modified -see Sec. 2.4. Finding eigenenergies and eigenstates of superconducting qubits invariably involves truncation of the infinite-dimensional Hilbert space. In each qubit class, this can be controlled by setting the appropriate truncation parameter, such as ncut for the transmon qubit. While typical values are suggested in each widget, convergence with respect to this cutoff (and similar cutoffs in other qubit classes) must be established by the user (see 2.1.1 below for a more detailed discussion).

Computing and plotting energy spectra
The energy eigenvalues of the transmon Hamiltonian are obtained by calling the eigenvals method. The optional parameter evals _ count specifies the desired number of eigenenergies: 1 tmon.eigenvals(evals _ count=12) Execution of this line yields a NumPy array of the lowest twelve energy eigenvalues. To plot eigenenergies as a function of one of the qubit parameters (here, EJmax, EC, flux, or ng), we generate an array of parameter values of interest, and then call the method plot _ evals _ vs _ paramvals. The latter takes as arguments the name of the parameter to be varied, an array of parameter values, and optionally the eigenvalue number. The following is an example for plotting the lowest six eigenenergies as a function of offset charge ng for 220 equally spaced points in the interval n g ∈ [−2, 2]: 1 ng _ list = np.linspace(-2, 2, 220) 2 tmon.plot _ evals _ vs _ paramvals('ng', ng _ list,  The full eigensystem consisting of both eigenvalues and eigenvectors is obtained through the method eigensys: 1 evals, evecs = tmon.eigensys() For the transmon qubit, this calculation is based on 1 scipy.linalg.eigh which performs diagonalization of the Hamiltonian matrix expressed in the charge basis. As dictated by the SciPy package, the eigenvector corresponding to the j-th eigenvalue is the Numpy array evecs.T[j].

Hilbert space truncation and convergence
As mentioned above, obtaining qubit spectra requires truncating the Hilbert space to some finite dimension. Users can control the Hilbert space dimension used by scqubits, by passing appropriate values for ncut and/or grid parameters to qubit class constructors 2 . These parameters effectively define the number of basis states that are used during diagonalization. In qubits with multiple degrees of freedom, users need to set a cutoff independently for each one.
Choosing cutoffs or grid sizes that are too small may lead to inaccurate results. The specific parameter values required for convergence naturally depend on the qubit type, circuit energies, as well as how many eigenenergies and/or eigenvectors the user wishes to obtain. While in some simple cases (e.g., transmon qubit -see appendix A.III) one can establish stringent cutoff requirements, for most circuits convergence must be established by trial and error.
A heuristic approach to this end is to repeat calculations with successive increases in cutoffs and grid sizes until results are essentially unchanged (within the desired accuracy). Once this is achieved, errors are typically limited by the default tolerances set by the scipy routines which scqubits internally uses for matrix diagonalization.

Plotting wavefunctions
The transmon qubit is a circuit with a single degree of freedom. Hence, its wavefunctions are readily plotted in the two bases natural for the transmon: the discrete charge basis ψ j (n) = n|ψ j , and 1 In some cases (e.g., 0 − π or the cos 2φ qubit), scqubits, uses scipy.sparse.linalg.eigsh for diagonalization. 2 In some qubits with multiple degrees of freedom, these variable names may be slightly modified to reflect which degrees of freedom are being addressed. See the API documentation [1] as well as the Appendices for details. the continuous phase basis ψ j (ϕ) = ϕ|ψ j . The TunableTransmon class offers two methods for this purpose: 1 tmon.plot _ n _ wavefunction(which=0, mode='real'); 2 tmon.plot _ phi _ wavefunction(which=[0,2,6], 3 mode='abs _ sqr'); Figure 2(b-c) shows the generated graphs. Above, the keyword argument which specifies the selection of wavefunctions to plot. In the first case, which =0 indicates that only the ground state wavefunction (index 0) should be displayed. In the second case, which is assigned a list consisting of multiple such indices, resulting in a plot showing several wavefunctions at once. The mode option determines how to plot the wavefunction which is generally complex-valued. The self-explanatory options are 'real', 'imag', 'abs', and 'abs _ sqr'.
For one-dimensional wavefunctions in ϕ basis ("position" basis), the plot mimics the textbook format widely used for 1d spectra: wavefunctions of the eigenstates are shown alongside the potentialenergy function, and are offset vertically by their corresponding eigenenergies.

Evaluating and visualizing matrix elements
Matrix elements of qubit operators play an important role in determining coupling strengths between a qubit and another quantum system. They are also critical for the sensitivity of the qubit to various sources of noise. In the case of the transmon qubit, for example, matrix elements of the charge operator are of frequent interest. This charge operator,n, is accessible through the class method n _ operator. Evaluation of the matrix elements for the current set of parameters is implemented by the method matrixelement _ table. An overview plot of matrix elements with respect to the lowest ten transmon eigenstates is obtained by 1 tmon.plot _ matrixelements('n _ operator', where the reference to the operator is provided in string format. The resulting plot is shown in Fig. 3.
Occasionally, it is further useful to plot matrix elements as a function of an external parameter.  This is accomplished by plot _ matelem _ vs _ paramvals Applied to the transmon charge matrix elements, we can easily visualize the dependence on the offset charge n g : 1 ng _ list = np.linspace(-2, 2, 220) 2 tmon.plot _ matelem _ vs _ paramvals('n _ operator', 3 'ng', ng _ list, Here, the names of both the operator and the external parameter are given as string arguments, followed by the array of values for the external parameter. Finally, select _ elems=4 specifies that all matrix elements ψ i |n|ψ j with 0 ≤ i, j ≤ 3 are requested in the plot. The output is shown in Fig.  4.

Units
scqubits allows the user to associate specific units with the energy values they provide when instantiating various qubit classes. The knowledge by scqubits of implied units is necessary for calculations that involve coherence time estimations (see Sec. 5), but also come in handy for automatic labeling of plot axes. All energies are assumed to be expressed in terms of frequencies (not angular frequencies), with the default being GHz.

Composite Hilbert spaces and interface with QuTiP
An important aspect of modeling superconducting circuits is the ability to study composite systems. scqubits provides an easy mechanism to explore setups that may consist of multiple qubits as well as harmonic (and weakly anharmonic) oscillators. Along with providing an easy way of constructing composite-system Hilbert spaces, and calculating and visualizing their many properties, scqubits also allows for easy exporting of effective Hamiltonians to QuTiP, an established toolbox for studying stationary and dynamical properties of closed and open quantum systems. At the heart of this functionality is the HilbertSpace class, which provides the data structures and methods for handling composite Hilbert space objects, and which we briefly explore in the sections below.
3.1 Example system: two transmons coupled to a harmonic mode Transmon qubits can be capacitively coupled to a common harmonic mode, realized by an LC oscillator or a transmission-line resonator. The Hamiltonian describing such a composite system is given byĤ where j = 1, 2 enumerates the two transmon qubits, E osc is the single-photon energy for the resonator. Furthermore,n j is the charge number operator for qubit j, and g j is the coupling strength between qubit j and the resonator. The first step consists of creating the objects describing the individual building blocks of the full Hilbert space. Here, these will be the two transmons and one oscillator: The significance of truncated _ dim lies in a simple hierarchical diagonalization scheme. Specifically, each subsystem is diagonalized separately in step one. Subsequently, the lowest few bare subsystem eigenstates up to truncation level truncated _ dim are fed forward into HilbertSpace.

Creating the HilbertSpace object
The desired HilbertSpace object can be created in two ways: either by utilizing the GUI displayed via hs = scq.HilbertSpace.create(), or programmatically by initializing the HilbertSpace object with a list of all subsystems and then specifying individual interaction terms:

Specifying interactions
Interaction terms describing the coupling between subsystems (or modifying a single subsystem itself) can be specified via the method add _ interaction in three different ways.

Operator-product based interface:
Interaction terms involving multiple subsystems S = 1, 2, 3, . . . are often of the form where the operatorsÂ j ,B j act on subsystem j.
(In the first case, the operatorsÂ j are expected to be hermitian.) This structure is captured in the following way: In this operator-product interface, op1, op2,. . . are either of the form (<array>, <subsystem>), i.e., a tuple with an array or sparse matrix in the first position, and the corresponding subsystem in the second position; or of the form <callable>, i.e., the operator is provided as a callable method which will automatically yield the subsystem the operator function is bound to. (These two choices can be mixed and matched.) The option add _ hc=True adds the hermitian conjugate to the specified interaction term.
Note that interactions based on only one operator are possible (simply drop all but the op1 entry). One example use case of this is the creation of a higher-order non-linearity a † a † a † aaa in a Kerr oscillator.

2.
String-based interface: The add _ interaction method can also be used to define the interaction in string form, by providing an expression that can be evaluated by the Python interpreter. Here, expr is a string used to define the interaction as a Python expression. It may use variables that are already defined globally, and operators given by the names provided in op1, op2, . . ..

Obtaining the Hamiltonian and spectrum
With the interactions specified, the full Hamiltonian of the coupled system can be obtained via the method hamiltonian, 1 dressed _ hamiltonian = hs.hamiltonian () which represents H in the basis of the bare product states composed of subsystem eigenstates. Since the Hamiltonian obtained this way is a proper Qobj, it can easily be handed over to QuTiP's time evolution routines such as mesolve. Eigenenergies and eigenstates can now either be obtained via the usual scqubits methods, 1 evals = hs.eigenvals() # or 2 evals, evecs = hs.eigensys () or by invoking QuTiP methods on the Hamiltonian itself, e.g., hs.hamiltonian.eigenstates().

Sweeping over external parameters
Determining the dependence of physical observables on one or multiple external parameter(s) is a common way to gain intuition for the properties and behavior of a system. Such parameter sweeps can be performed with scqubits on multiple levels: (1) at the level of a single qubit, and (2) at the level of a composite quantum system. At the single-qubit level, each qubit class provides several methods that enable producing parameter sweep data and plots. Central quantities of interest, in this case, are energy eigenvalues and matrix elements -in particular, their dependence on parameters like flux or offset charge. The relevant methods available for every implemented qubit class are: • For flexible parameter scans, scqubits provides the ParameterSweep class. To illustrate its usage, we reuse a composite Hilbert space akin to the one presented above: two tunable transmon qubits capacitively coupled to an oscillator. The ParameterSweep class facilitates computation of spectra as function of one or multiple external parameter(s). For efficiency in computing a variety of derived quantities and creating plots, the computed bare and dressed spectral data are stored internally. A ParameterSweep object is initialized by providing the following parameters: 1. hilbertspace: a HilbertSpace object that describes the quantum system of interest 2. paramvals _ by _ name: a dictionary that maps each parameter name (string) to an array of parameter values 3. update _ hilbertspace: a function that defines how parameter changes affect the system 4. subsys _ update _ info: (optional) for potential speed-up, specify which subsystems undergo changes as each of the parameters is varied 5. deepcopy: (optional) determines whether the HilbertSpace object and all constituents should be duplicated and disconnected from the global objects 6. num _ cpus: (optional) number of CPU cores requested for the sweep evaluation These ingredients all enter as initialization arguments of the ParameterSweep object. Once initialized, spectral data is generated and stored.
In our example, we consider the strength of a global magnetic field as the parameter to be changed. This field determines the magnetic fluxes for both qubits, in proportions according to their SQUID loop areas. We will reference the flux for transmon 1, and express the flux for transmon 2 in terms of it via an area ratio. In addition, we will vary the offset charge of transmon 2.
The following code illustrates this functionality: In the code above, update _ hilbertspace directly manipulates transmon instances via their global references. Alternatively, HilbertSpace constituents can be accessed via sweep.hilbertspace [<id _ str>] where id _ str is a string identifier either provided explicitly at initialization of object instances, or autogenerated by scqubits. (Details and examples of this functionality are available in the documentation and example notebooks -see Sec. 7.)

Generated ParameterSweep data
Much of the computed data that is stored and immediately retrievable after this sweep. The stored data is accessed as if ParameterSweep were a Python dict with the following string keys: 1. "evals", "evecs": dressed eigenenergies and eigenstates 2. "bare _ evals", "bare _ evecs": bare eigenenergies and eigenstates for each subsystem 3. "lamb", "chi", "kerr": dispersive energy coefficients Data are returned as a NamedSlotsNdarray, a subclass of the regular NumPy ndarray. It features several convenient new slicing options, such as slicing by parameter name, or slicing by reference to a particular parameter value. See [2] for a more more detailed discussion as well as examples.

Transition plots
Energy spectra obtained in single-tone or two-tone spectroscopy always represent transition energies, rather than absolute energies of individual eigenstates. To generate data mimicking this, appropriate differences between eigenenergies must be taken.
The methods for generating transition energy data and plotting them are transitions and plot _ transitions. To create a plot, we first "preslice" the ParameterSweep instance which specifies a sweep along a single axis.
Then, the plot _ transitions method can be called 3 : When coloring transitions according to the dispersive-limit product state labels, one potential artifact is nearly unavoidable: whenever states undergo avoided crossings and the dispersive limit breaks down, coloring must discontinuously switch from one branch to another. scqubits attempts to interrupt coloring in such regions. However, if the avoided crossing occurs over a range comparable to the parameter value spacing, then discontinuities from connecting separate branches will remain visible.

Transition plot options:
The generated transition plot above is based on a number of default settings, including: (i) the origin of each transition is the system's ground state, (ii) single-photon transitions are plotted in light grey and (iii) transitions within each individual subsystem are marked separately in color and accompanied by a legend. This is possible in regions where the dispersive approximation holds, i.e., hybridization between subsystems remains weak. Labels in the legend are excitation levels of individual systems: ((0,0,0), (1,0,0)) denotes a transition from the ground state to the state with subsystem 1 in the first excited state, and subsystems 2 and 3 in their respective ground states.
Many aspects of transition energy plots can be changed. The following illustrates a subset of options that change which transitions are plotted and how. The computed data is subsequently accessible via 1 <ParameterSweep>["custom name"]

Estimation of coherence times
scqubits provides extensive functionality that allows users to estimate coherence times of various qubits. Very general methods t1 and tphi are implemented for each noisy qubit, which can be used for calculations of depolarization as well as pure dephasing times for almost arbitrary noise processes (see sections 5.1 and 5.2 for more details). Furthermore, a large variety of predefined methods corresponding to more specific noise channels (e.g., dephasing due to 1/f charge noise, or depolarization from coupling to a transmission line) are implemented as well (for a list, see Table 1). Due to different qubit properties and circuit topologies, each qubit is affected by some specific subset of these predefined noisy processes. To see which channels are currently implemented for any given qubit, one can run the command (here shown for the case of a transmon) This returned list contains methods with selfexplanatory names that either start with t1 or tphi and represent the depolarization or pure dephasing processes, respectively. Each of these methods can then be called directly to obtain the corresponding coherence time (which will take into account the current units setting -see section 2.4). For example, in order to calculate the pure dephasing time due to 1/f flux noise, one can simply execute 1 tmon.tphi _ 1 _ over _ f _ flux() Each available method can be provided with a set of custom parameters that give the user a means to fine-tune various noise process properties (e.g., bath temperature, 1/f flux-noise strength, etc.). There are also common method options that specify which qubit levels should be considered (levels 0 and 1 are assumed by default), whether a rate instead of a time should be returned, or whether depolarization rates should include both upwards and downwards transitions. A more complicated method call could therefore take the form 1 tmon.t1 _ charge _ impedance(i=3, j=1, This returns a depolarization rate (not time) between levels 3 and 1, due to the qubit coupling to a 50 Ω transmission line, at a temperature of 100 mK. In this case, the impedance (parameter Z) can be a constant, or alternatively an angular frequencydependent Python function that will be evaluated at the frequency difference between levels i and j.
In the next sections we briefly outline the basic physics and assumptions underlying the estimation of the various pure dephasing and depolarization times.

Dephasing due to 1/f noise
Dephasing noise leads to loss of coherence, i.e., the relative phases relevant in superpositions of multiple states are lost over time. One of the most important kinds of noise affecting superconducting qubits is 1/f noise, which leads to slow fluctuations of the energy-level spacing. The spectral density function characterizing this noise is given by Here, A λ corresponds to the amplitude or strength of the particular noise channel λ, such as charge or flux. scqubits uses sensible default values for this quantity based on the literature. Alternative values, however, can be set when provided by the user. The pure dephasing time due to a noise channel labeled λ (away from sweet spots 4 ) is given by [12,16] with t exp representing the measurement time, and ω low the low-frequency cutoff. (If not provided by the user then defaults of t exp = 10 µs and ω low /2π = 1 Hz are used.) As already hinted above, some qubits provide predefined methods for estimating effects due to 1/f flux, charge as well as Josephson junction critical-current noise channels, named tphi _ 1 _ over _ f _ flux, tphi _ 1 _ over _ f _ charge and tphi _ 1 _ over _ f _ cc respectively. Each qubit also implements a more general method tphi _ 1 _ over _ f which accepts a user-defined operator ∂ λĤ (for an arbitrary, user-defined noise channel λ), that is then internally used by scqubits to implement Eq. 3.

Depolarization
Noise may also cause depolarization of the qubit by inducing spontaneous transitions among eigenstates. scqubits uses the standard perturbative approach (Fermi's Golden Rule) to approximate the resulting transition rates due to different noise channels. The rate of a transition from state i to state j can be expressed as whereB λ is the noise operator, and S(ω ij ) the spectral density function evaluated at the angular frequency associated with the transition frequency, ω ij = ω j − ω i . ω ij , which is positive in the case of decay (the qubit emits energy to the bath), and negative in case of excitations (the qubit absorbs energy from the bath). Unless stated otherwise (see channel-specific documentation [2]), it is assumed that the depolarizing noise channels satisfy detailed balance, which implies where T is the bath temperature, and k B Boltzmann's constant. By default, all t1 methods estimate the depolarization times based on the sum of the upward and downward rates. This behavior is controlled by the argument total, which can be modified by the user. For example, setting total= False will calculate only a single-directional transition rate from the state indexed i to the state indexed j. As in the case of 1/f dephasing, each qubit implements a subset of predefined depolarization methods such as t1 _ capacitive or t1 _ flux _ bias _ line for example (see [2] for details on what methods are implemented for each qubit). Coherence times due to arbitrary depolarization processes can also be readily calculated. This is done using a general method t1, which accepts an arbitrary, user-constructedB λ operator, as well as a custom spectral density function S(ω).

Effective coherence times
Coherence times observed in experiments will typically be due to the combined effect of multiple contributing noise channels. scqubits can easily combine channels and compute effective coherence times or rates. In the case of depolarization, the effective coherence time is obtained from where the sum runs over all default noise channels, i.e., those methods with names beginning with t1). For each qubit, the included default channels can be listed by calling the qubit's effective _ noise _ channels method. To calculate the effective T eff 1 time based on the default noise channels and their parameters, one simply executes 1 tmon.t1 _ effective() Similarly, scqubits can calculate an effective dephasing time T eff 2 (using the method t2 _ effective). This time scale includes contributions from both pure dephasing as well as depolarization and is defined as Once again the k index cycles over the set of default noise channels. Both t1 _ effective as well as t2 _ effective can be easily customized, so that only select noise channels are incorporated, or calculations be based on specific user-defined parameters. For example, a smaller set of channels can be specified by passing a list of channel methods. Further, options shared by all noise channels can be set via the common _ noise _ options keyword argument, which accepts a dictionary of options. This is illustrated in the following example, where the temperature is set to T = 0.050 K: In addition to common _ noise _ options, channelspecific noise options can be provided. This is accomplished by replacing the name of the noise-channel method with a tuple of the form ( channel _ name, noise _ options) with noise _ options a Python dictionary. In the example below, we calculate an effective T eff 2 , using a non-default value for the 1/f flux-noise strength A _ flux that internally gets passed to the qubit's tphi _ 1 _ over _ f _ flux method:  Capacitive loss due to dielectric dissipation [22,24] t1 _ charge _ impedance Loss due to charge-coupling to an impedance (e.g., open transmission line) [13,23] t1 _ flux _ bias _ line Loss due to current fluctuations in the flux-bias line [11,16] t1 _ inductive Inductive loss due to quasiparticle tunneling in Josephson junction chains that are used to implement superinductances [22,24] t1 _ quasiparticle _ tunneling Loss due to quasiparticle tunneling across a single Josephson junction [22,24] T φ pure-dephasing processes tphi _ 1 _ over _ f _ cc Dephasing due to 1/f critical-current noise (fluctuations of E J in a Josephson junction) [11,13,16] tphi _ 1 _ over _ f _ charge Dephasing due to 1/f charge noise (fluctuations in charge offset) [11,13,16] tphi _ 1 _ over _ f _ flux Dephasing due to 1/f flux noise (fluctuations in the applied magnetic flux) [11,13,16]

Coherence visualization
A common way to understand and visualize how noise affects a given qubit, is to plot decoherence times as a function of one of the external parameters, such as flux, charge or one of the qubit internal energy parameters, say E J , for example. Each qubit provides a flexible method called plot _ coherence _ vs _ paramvals, which facilitates this functionality. To provide an overview of the dependence of decoherence properties on, say, flux, the effect of all supported noise channels (as defined by each quibit's supported _ noise _ channels method), can be visualized in individual plots as follows: 1 tmon.plot _ coherence _ vs _ paramvals( Here, param _ vals is an array of flux values for which coherence data is generated and plotted. The set of plots to be included can be determined by the user and further customized by either passing various plot options directly to the plot _ coherence _ vs _ paramvals method, or by manipulating the properties of the matplotlib Axes object returned by plot _ coherence _ vs _ paramvals: The resulting graphical output (given a tmon object, as defined in Sec. 2), is shown in Fig. 5. In this example, only plots for noise channels tphi _ 1 _ over _ f _ flux and t1 _ capacitive are included. The optional scale argument defines a custom yaxis scale factor. Note that the default units of GHz will automatically yield coherence times in units  of ns. Setting scale=1e-3 and adjusting the ylabel of the plot allows us to switch to µs on the fly. Customizing plots in such ways helps making plots visually appealing and publication-ready, without globally changing unit setting (see section 2.4). Finally, scqubits streamlines visualization of the effective noise introduced in section 5.3 above through the methods plot _ t1 _ effective _ vs _ paramvals and plot _ t2 _ effective _ vs _ paramvals.
Both methods work in a way analogous to plot _ coherence _ vs _ paramvals, except that now a single plot of the effective coherence time is presented. Provision of an optional noise _ channels list now sets the specific noise channels included in the calculation of T eff 1 and T eff 2 .

Interactive exploration
Exploring the properties of coupled quantum systems can benefit from visual aides, such as inspection how observables change when system parameters are modified. The Explorer class in scqubits provides multiple interactive viewgraphs collecting an important set of information regarding the userdefined system of interest. The Explorer is based on a HilbertSpace object describing the composite circuit-QED system of interest. As explained in Section 4, ParameterSweep can then be used to record a sweep of an external tuning parameter. This could be, for example, an external magnetic flux, an offset charge, or even a circuit parameter such as a capacitance (difficult to change in-situ in an experiment).

Example: fluxonium coupled to a resonator
As a concrete example, we consider a system composed of a fluxonium qubit, coupled through its charge operator to the voltage inside a resonator. The initialization of the composite Hilbert space proceeds as usual: we first define the individual two subsystems that will make up the Hilbert space, Here, the truncated _ dim parameters are for the hierarchical diagonalization of the composite Hilbert space. For the fluxonium subsystem, cutoff fixes the internal Hilbert space dimension to 110. Once diagonalized, however, only a few eigenstates are usually meant to be retained and included in the composite Hilbert space description. In the example above, the lowest nine states are selected. Similarly, we retain five levels of the resonator, i.e., photon states n = 0, 1, . . . , 4 are included.
Next, the two subsystems are declared as the two components of a joint Hilbert space: The interaction between fluxonium and resonator is of the formĤ int = gn(a + a † ), wheren is the fluxonium's charge operator, fluxonium.n _ operator: At this point, we may start the interactive Explorer class; sample output is shown in Fig. 6.
scqubits is an open source package, and all of its source code is freely available online under the BSD-3 license. The package is divided into three separate repositories: the first contains the source code of the core package, the second the Sphinx code for the online documentation, and finally the third collects example jupyter notebooks that illustrate how various features of scqubits can be used. To install scqubits, users can either clone the main github repository directly and install via pip install ., or alternatively download and install the package through the PyPI or Anaconda package repositories (see [2]). For a list of online links to the various github pages containing all the source code, online documentation, live example jupyter notebooks, as well as PyPi and Anaconda package index repository pages, see table 2. Users, are welcome and encouraged to file bug reports, post comments and suggestions, as well as initiate pull requests on the relevant github pages. Finally, users who find scqubits useful, are encouraged to cite this paper. The relevant information can be readily accessed by executing the cite function, like so:

Conclusions
With this paper, we have introduced the scqubits library: an open-source Python toolbox enabling the simulation of superconducting qubits -both at the single-qubit level, and at the level of composite quantum systems consisting of multiple qubits and resonators. The current functionality of the library encompasses a broad range of computational Bare spectra of the individual qubits, 2) Wavefunctions of the bare qubits, 3) Dressed spectrum of the composite Hilbert space, 4) Spectrum for n-photon qubit transitions, starting from a given initial state, 5) AC Stark shift χ 01 for any of the qubits, and 6) Charge matrix elements for any of the qubits, using the same initial state as in point 4) tasks and visualization tools commonly used in research involving superconducting qubits. While maintaining the existing interface, future work will aim to extend the scope of the library, for example by including a systematic workflow for fitting experimental two-tone spectroscopy data, analyzing custom circuits defined by the user, implementing newer qubit designs, as well adding to the list of predefined noise channels in order to extend coherence time estimations.

Superconducting qubit and oscillator classes
scqubits currently implements several common types of superconducting qubits along with a linear and non-linear oscillators, as well as a generic qubit (i.e., a simple two-level system), each realized as a Python class 5 . A brief summary is shown in the following table [2]: class description

3-junction flux qubit [19]
ZeroPi 0-π qubit (symmetric) [5,8] FullZeroPi 0-π qubit (coupled to ζ-mode) [8] Cos2PhiQubit cos 2φ qubit [24] Oscillator Quantum harmonic oscillator KerrOscillator Nonlinear Kerr oscillator GenericQubit A two-level system All the qubit classes (except for GenericQubit) define a number of important methods that can be used for diagonalization, computation of matrix elements and spectral data, as well as for plotting. A few of these are summarized in Table 3.
Besides these methods, each superconducting qubit class also implements a predefined number of quantum operators which can simplify doing various calculations. These may include the phaseφ or numbern operators, but also more specialized ones such as cos(φ) and the like. See the API documentation [1], for a comprehensive list.
In the following sections we describe the qubit and oscillator classes that scqubits implements in more detail: present their circuit diagrams (where applicable), give definitions of the respective Hamiltonians, and briefly discuss their numerical implementation in the library.

A.III Transmon qubit
The Cooper pair box [20] and transmon qubit [16] share the same underlying circuit composed of a single Josephson junction and a parallel capacitance which may either be the pure junction capacitance, or include an external shunt capacitance: The circuit Hamiltonian can be written aŝ wheren (φ) is the charge (phase) operator 7 Here E C = e 2 /2C Σ is the charging energy associated with the combined capacitances of the junction, the shunt capacitor, and any additional ground capacitance and/or capacitance to a charge bias line, The Josephson energy of the junction is related to its critical current via E J = I c Φ 0 /2π. The quantity n g is the dimensionless offset charge capturing the capacitive coupling to a bias voltage source as well as electric-potential fluctuations of the environment. Internally, the Transmon class employs the chargebasis representation to construct the Hamiltonian matrix. This matrix is infinite-dimensional, in principle, and must hence be truncated. To this end, a charge-number cutoff ncut is introduced. Given this cutoff, the included charge states |n are within the range -ncut≤ n ≤ncut. The cutoff must be chosen sufficiently large to avoid truncation errors 8 . 7 Since the Hamiltonian is periodic inφ, we stress that only periodic functions ofφ are formally defined. 8 For low-lying transmon wavefunctions (EJ /EC 1), there is a simple cutoff criterion: The ground state is close to Gaussian with standard deviation σ = (8EC /EJ ) 1/4 . Treating n as continuous and assuming an ng of order 1, Fourier transform of the ground state yields ψ0(n) which is also close to a Gaussian, here with standard deviation σ = (EJ /8EC ) 1/4 . Using a 3σ estimate, one finds that ncut should be no smaller than ncutmin ≈ 2(EJ /EC ) 1/4 .

A.IV Fluxonium qubit
The circuit of the fluxonium qubit [18] consists of a Josephson junction shunted by a large inductor: The resulting qubit Hamiltonian takes the form where E C = e 2 /2C J and E J are the charging and Josephson energies of the junction, and E L = ( Φ 0 2π ) 2 /L is the inductive energy. The Hilbert-space basis used by scqubits for construction of the Hamiltonian is the harmonic-oscillator basis associated with the inductor and junction capacitor. The Hamiltonian can be rewritten in this basis by employing the usual ladder operators a, a † , Here, we have rewritten cos(φ − ϕ ext ) = 1 2 e −iϕext e iφ + h.c., and usedφ = ϕosc √ 2 (â +â † ) with ϕ osc = (8E C /E L ) 1/4 denoting the oscillator length. Numerical evaluation of the matrix exponential in (A.6) is carried out via scipy.linalg.expm().
It must be noted that the harmonic-oscillator basis is not well-adapted to fluxonium eigenstates that localize in individual wells of the cosine contribution to the potential. Consequently, the inevitable truncation in the harmonic-oscillator basis must generally proceed with caution, and the cutoff number cutoff be chosen sufficiently large for convergence.

A.V Flux qubit
The 3-junction flux qubit that scqubits implements, was first proposed in [21]. Its circuit consists of 3 Josephson junctions in a loop that is threaded with an external flux Φ ext .
E Jk (C Jk ) with k ∈ {1, 2, 3} are the Josephson energies (capacitances) associated with each junction. For normal qubit operations one junction is chosen to be smaller than the other two. The effective Hamiltonian of such a flux qubit is described bŷ In the above expression, E C represents a charging energy matrix (which includes effects of small capacitors C g ), while n gj (with j ∈ {1, 2}) are the charge offsets. The two degrees of freedom are represented byφ 1 andφ 2 , with their conjugates charge operatorsn 1 andn 2 respectively. Numerical diagonalization is performed in the charge basis for botĥ ϕ 1 andφ 2 (see discussion in A.III).
A sample initialization code of the Flux qubit, where the third junction is assumed to be smaller than the other two by a factor of α, is shown below

A.VI 0-π qubit
The 0 − π qubit was first proposed in Ref. [5]. Its circuit consists of two Josephson junctions, two capacitors as well as two superinductors arranged in the following form: The behavior of this qubit has been shown to strongly depend on the presence of parameter disorder [8,11]. In particular, when the two (nonjunction) capacitive or inductive energies are not identical (i.e., C 1 = C 2 or L 1 = L 2 ), the core qubit degrees of freedomθ andφ end up coupling to a low-frequency harmonic modeζ, which would stay uncoupled when there is no disorder in in these parameters. For this reason, scqubits implements two separate classes that can be used for modeling of 0 − π qubits. The first, ZeroPi, assumes both of the inductors as well as the (non-junction) capacitors are identical and hence only includes theθ andφ degrees of freedom. The second, FullZeroPi, allows for small parameter disorder in the superinductnaces L k as well as the (non-junction) capacitors C k , resulting in a three degrees of freedom system, which now also includes the harmonicζ mode. Note that the ZeroPi class does allow disorder in the Josephson junctions (i.e., E J1 and C J1 do not have to be identical to E J2 and C J2 ), as this kind of disorder still leaves theζ mode decoupled fromθ andφ.
To quantify disorder in any parameter X (with X ∈ {C, E L , C J , E J }, we define Then, in the limit of small disorder dX [8,11], a general 0 − π Hamiltonian can be approximated bŷ H tot =Ĥ 0−π +Ĥ int +Ĥ ζ (A.9) withĤ 0−π =2E CJn 2 φ + 2E CΣ (n θ + n g ) 2 (A.10) along witĥ The quantity n g is the charge offset ofn θ , while ϕ ext = 2πΦ ext /Φ 0 . The class ZeroPi only imple-mentsĤ 0−π , as it clear from the above description, that when dE L = dC = 0,Ĥ int = 0, , while FullZeroPi includes all three terms ofĤ tot . Internally, we use charge basis to describe theθ degree of freedom (see A.III), phase basis for theφ degree of freedom, and finallyζ is modeled using harmonic basis (see A.IV). For problems where the disorder dC and dE L can be neglected, it is strongly recommended that ZeroPi is used, as the performance penalty from including the physics of theζ mode can be substantial. Below is sample code showing an initialization of a 0−π qubit. Here, we include a 5% disorder in the Josephson junction energies, but assume there is no disorder in the superinductors and (non-junction) capacitors, allowing us to use the ZeroPi class: A.VII cos 2φ qubit The cos 2φ qubit was proposed in [24]. Its circuit includes a superconducting loop, threaded by an external flux Φ ext , which consists of two superinductors and two Josephson junctions, with an appropriately placed shunt capacitance: Such topology can be engineered to only allow pairs of cooper pairs to tunnel, leading a level of intrinsic protection from noise [24]. The Hamiltonian of a cos 2φ qubit can be written aŝ H cos 2φ = 2E CJn 2 φ + 2E CJ (n θ − n g −n ζ ) 2 + 4E Cn Parameters of the form dX represent disorder 9 . In particular, we have dX = X 1 − X 2 X 1 + X 2 , (A.14) for X ∈ {L, C J }, with L (C J ) being the superinductance (junction capacitance). We define the charging (inductive) energy of the Josephson junction capacitor (superinductor) as E CJk = e 2 /2C Jk (E Lk = ( Φ 0 2π ) 2 /L k ). These expressions let us further write with Y ∈ {E L , E CJ }. Finally, the Josephson energy satisfies Similarly to the 0 − π qubit, cos 2φ circuit consists of three degrees of freedom: θ, φ and ζ, with their respective conjugates defined asn θ ,n φ andn ζ . The quantity n g represents the charge offset of then θ variable. We stress that the labels and notation utilized by scqubits for the degrees of freedom and some of the circuit parameters differs slightly from [24]. The following table outlines how to convert between the two: 9 Note the difference in the definition of disorder here versus the definition for the 0 − π qubit. The definition of parameter disorder in the cos 2φ qubit used by scqubits follows the notation in [24].

scqubits
Reference [24] ζ θ θ ϕ To numerically diagonalize the Hamiltonian of the cos 2φ qubit, the harmonic basis are used for both theφ andζ variables (see the discussion of the Fluxonium qubit in A.IV for more details), while the charge basis are used for theθ variable (see A.III). When instantiating an scqubits object corresponding to this qubit, the user needs to specify cutoffs for basis states described above (this is done using phi _ cut, zeta _ cut, and ncut), which need to be chosen large enough so that convergence is achieved.
An instance of the cos 2φ qubit can be initialized as follows 1 cos2phi _ qubit = scq.Cos2PhiQubit(