Composite Hilbert Spaces & QuTiP

The HilbertSpace class provides data structures and methods for handling composite Hilbert spaces which may consist of multiple qubits and/or oscillators coupled to each other. To harness the power of QuTiP (a toolbox for simulating stationary and dynamical properties of closed and open quantum systems), HilbertSpace provides a convenient interface: it generates qutip.qobj objects which can then be handed over directly to QuTiP.

In the following, the basic functionality of the HilbertSpace class is illustrated for the example of 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:

\begin{equation} H=H_\text{tmon,1} + H_\text{tmon,2} + \omega_r a^\dagger a + \sum_{j=1,2}g_j n_j(a+a^\dagger), \end{equation}

where \(j=1,2\) enumerates the two transmon qubits, \(\omega_r\) is the (angular) frequency of 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.

Creating a HilbertSpace instance

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:

tmon1 = scq.Transmon(
    EJ=40.0,
    EC=0.2,
    ng=0.3,
    ncut=40,
    truncated_dim=4     # after diagonalization, we will keep 3 levels
)

tmon2 = scq.Transmon(
    EJ=15.0,
    EC=0.15,
    ng=0.0,
    ncut=30,
    truncated_dim=4
)

resonator = scq.Oscillator(
    E_osc=4.5,
    truncated_dim=4  # up to 3 photons (0,1,2,3)
)

The HilbertSpace object is now created in one of two ways. The first is by utilizing the GUI:

g = 0.1
hilbertspace = scq.HilbertSpace.create()

The second way of creating a HilbertSpace object is through regular programmatic creation: provide a list of all subsystems and then specify individual interaction terms (see the next subsection for the latter):

hilbertspace = scq.HilbertSpace([tmon1, tmon2, resonator])
print(hilbertspace)

HilbertSpace:  subsystems
-------------------------

Transmon------------|
                    | EJ: 40.0
                    | EC: 0.2
                    | ng: 0.3
                    | ncut: 40
                    | truncated_dim: 4
                    |
                    | dim: 81


Transmon------------|
                    | EJ: 15.0
                    | EC: 0.15
                    | ng: 0.0
                    | ncut: 30
                    | truncated_dim: 4
                    |
                    | dim: 61


Oscillator----------|
                    | E_osc: 4.5
                    | l_osc: None
                    | truncated_dim: 4
                    |
                    | dim: 4

Specifying interactions between subsystems

Interaction terms describing the coupling between subsystems can be specified in three different ways.

Option 1: operator-product based interface

Interaction terms involving multiple subsystems \(S=1,2,3,\ldots\) are often of the form

\(V = g A_1 A_2 A_3 \cdots \qquad \text{or} \qquad V = g B_1 B_2 B_3 \cdots + \text{h.c.}\)

where the operators \(A_j\), \(B_j\) act on subsystem \(j\). In the first case, the operators \(A_j\) are expected to be hermitean.

This structure is captured by using the add_interaction method in the following way:

g1 = 0.1  # coupling resonator-tmon1 (without charge matrix elements)
operator1 = tmon1.n_operator()
operator2 = resonator.creation_operator() + resonator.annihilation_operator()

hilbertspace.add_interaction(
    g=g1,
    op1=(operator1, tmon1),
    op2=(operator2, resonator)
)

g2 = 0.2  # coupling resonator-CPB2 (without charge matrix elements)
hilbertspace.add_interaction(
    g=g2,
    op1=tmon2.n_operator,
    op2=resonator.creation_operator,
    add_hc=True
)

In this instance, op1 and op2 take either

1. (<array>, <subsystem>): a tuple with an array or sparse matrix in the first position and the subsystem in the second position, or

2. <callable>: the operator as a callable function (which will automatically yield the subsystem the operator function is bound to).

These two choices can be mixed and matched.

Note:

• Interactions based on only one operator are enabled (simply drop the op2 entry). One use case of this is to create a higher order non-linearity \(a^\dagger a^\dagger a^\dagger a a a\) in a Kerr oscillator.

• Interactions with more than two operators are built by specifying op3, op4, etc.

The option add_hc = True adds the hermitian conjugate to the Hamiltonian, and is available in the present operator-product based interface, as well as in the one described next as option 2.

Option 2: string- based interface

The add_interaction method can also be used to define the interaction in string form, providing an expression that can be evaluated by the Python interpreter (after operator name substitution).

hilbertspace_example = scq.HilbertSpace([tmon1, tmon2, resonator])

g3 = 0.1

hilbertspace_example.add_interaction(
    expr="g3 * cos(n) * adag",
    op1=("n", tmon1.n_operator),
    op2=("adag", resonator.creation_operator),
    add_hc=True
)

expr is a string used to define the interaction as a Python expression. It should use variables that are already defined globally, and operators given by the names provided in op1, op2, op3, etc. When using expr, each argument op1, op2, … should be a tuple of the form (<operator name as string>,  <operator as callable>) or (<operator name as string>,  <operator as array | sparse matrix>, <subsystem>).

Beyond built-in Python functions and operations, acceptable functions to use in expr are:

Function	Function translation
'cos'	'Qobj.cosm'
'sin'	'Qobj.sinm'
'exp'	'Qobj.expm'
'sqrt'	'Qobj.sqrtm'

Option 3: Qobj interface

Finally add_interaction can also been used to directly add a Qobj object that has already been properly identity wrapped:

# Generate a Qobj
g = 0.1
a = qt.destroy(4)
kerr = a.dag() * a.dag() * a * a
id = qt.qeye(4)
V = g * qt.tensor(id, id, kerr)

hilbertspace_example.add_interaction(qobj = V)

In this case, the add_interaction method only requires the argument qobj which must be an object of the type Qobj.

Obtaining the Hamiltonian

With the interactions specified, the full Hamiltonian of the coupled system can be obtained via the method .hamiltonian(). This conveniently results in a qubit.Qobj operator:

dressed_hamiltonian = hilbertspace.hamiltonian()
dressed_hamiltonian

Quantum object: dims = [[4, 4, 4], [4, 4, 4]], shape = (64, 64), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}-48.968 & 0.030 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.030 & -44.468 & 0.042 & 0.0 & 0.261 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.042 & -39.968 & 0.052 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.052 & -35.468 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.261 & 0.0 & 0.0 & -44.881 & \cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & -4.723 & 0.0 & 0.0 & -0.749 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & -14.477 & 0.030 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.030 & -9.977 & 0.042 & 0.0\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & -0.749 & 0.0 & 0.042 & -5.477 & 0.052\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \cdots & 0.0 & 0.0 & 0.0 & 0.052 & -0.977\\\end{array}\right)\end{equation*}

This matrix is the representation of the Hamiltonian in the bare product basis.

The non-interacting (“bare”) Hamiltonian and the interaction Hamiltonian can be accessed similarly by using the methods .bare_hamiltonian() and .interaction_hamiltonian().

Obtaining the eigenspectrum via QuTiP

Since the Hamiltonian obtained this way is a proper qutip.qobj, all QuTiP routines are now available. In the first case, we are still making use of the scqubit HilbertSpace.eigensys() method. In the second, case, we use QuTiP’s method .eigenenergies():

evals, evecs = hilbertspace.eigensys(evals_count=4)
print(evals)

[-48.97770317 -45.02707241 -44.36656205 -41.18438832]

dressed_hamiltonian = hilbertspace.hamiltonian()
dressed_hamiltonian.eigenenergies()

array([-48.97770317, -45.02707241, -44.36656205, -41.18438832,
       -41.1776098 , -40.46448065, -39.76202396, -37.45533167,
       -37.22705488, -36.65156212, -36.56694139, -35.88617024,
       -35.14255357, -33.58960343, -33.38439485, -33.12061816,
       -32.66494366, -32.04065582, -31.96284558, -30.93536847,
       -29.65538466, -29.63912664, -28.97946785, -28.85287223,
       -28.64748897, -28.09427075, -27.35745603, -26.97461415,
       -26.21586637, -25.79648462, -25.3216198 , -25.07747135,
       -24.37567532, -24.24760531, -23.66618527, -23.15199748,
       -22.2580349 , -22.06752989, -21.57303944, -21.26626732,
       -20.85288717, -20.51511596, -19.78521899, -19.1938092 ,
       -18.41228666, -17.73463904, -17.66080329, -16.93618597,
       -16.66699841, -15.89010361, -15.58181323, -14.68288545,
       -13.84725844, -13.27037775, -13.08388811, -12.34366698,
       -11.62655644, -10.308913  ,  -9.23168892,  -8.32791939,
        -8.13883089,  -5.82301882,  -4.18180351,  -0.88957289])

Spectrum lookup: converting between bare and dressed indices

We can refer to eigenstates and eigenenergies of the interacting quantum systems in two ways:

1. dressed indices: enumerating eigenenergies and -states as j=0,1,2,... starting from the ground state

2. bare product state indices: for a Hilbert space composed of three subsystems, specify tuples of the form (l,m,n) where l, m, n denote the excitation levels of each bare subsystem.

To use lookup functions converting back and forth between these indexing schemes, first generate the spectrum lookup table via:

hilbertspace.generate_lookup()

This makes lookup functions available via hilbertspace.lookup.<lookup method> where

lookup method	purpose
bare_index	yields bare product state index for given dressed index
dressed_index	yields dressed index for given bare product state index
energy_dressed_index	return eigenenergy for given dressed index
energy_bare_index	return eigenenergy for given bare product-state index
bare_eigenenergies	return bare eigenenergies for specified subsystem
bare_eigenstates	return bare eigenstates for specified subsystem
bare_productstate	return bare product state for given bare index

Here are the bare energies of the first Cooper pair box system:

hilbertspace.lookup.bare_eigenenergies(tmon1)

array([-36.05064983, -28.25601136, -20.67410141, -13.31487799])

The dressed state with index j=8 corresponds to following bare product state index:

hilbertspace.lookup.bare_index(8)

(1, 1, 0)

The bare product state (2,1,0) most closely matches the following dressed state:

hilbertspace.lookup.dressed_index((2, 1, 0))

21

This is the eigenenergy for dressed index j=10:

hilbertspace.lookup.energy_dressed_index(10)

-36.56694138527015

Helper function: identity wrapping

The abbreviated notation commonplace in physics usually omits tensor products with identities. As an example, in a composite Hilbert space consisting of three subsystems \(j=1,2,3\) with individual Hamiltonians \(H_j\), we typically use the notation

\(H=H_1+H_2+H_3\)

which is meant to stand for

\(H = H_1\otimes\mathbb{I}\otimes \mathbb{I} + \mathbb{I}\otimes H_2\otimes \mathbb{I} + \mathbb{I}\otimes\mathbb{I}\otimes H_3\)

In most places of the HilbertSpace class, scqubits manages this “identity wrapping” automatically. In some use cases, it may become necessary to invoke this functionality manually. For that purpose, scqubits implements the helper routine identity_wrap. For instance, the term \(\mathbb{I}\otimes H_2\otimes \mathbb{I}\) would be generated via:

scq.identity_wrap(H2, subsystem2, [subsystem1, subsystem2, subsystem3])

The complete interface is:

scq.identity_wrap(operator, subsystem, subsys_list, op_in_eigenbasis=False, evecs=None)

Wrap given operator in subspace subsystem in identity operators to form full Hilbert-space operator.

Parameters

• operator: operator acting in Hilbert space of subsystem; if str, then this should be an operator name in the subsystem, typically not in eigenbasis

• subsystem (QuantumSys): subsystem to which operator belongs

• subsys_list: list of all subsystems (in fixed order!)

• op_in_eigenbasis (bool, optional): whether operator is given in the subsystem eigenbasis; otherwise, the internal QuantumSys basis is assumed

• evecs (ndarray, optional): internal QuantumSys eigenstates, used to convert operator into eigenbasis

Return type: Qobj