Using Tensor Products and Partial Traces¶

Tensor products¶

To describe the states of multipartite quantum systems - such as two coupled qubits, a qubit coupled to an oscillator, etc. - we need to expand the Hilbert space by taking the tensor product of the state vectors for each of the system components. Similarly, the operators acting on the state vectors in the combined Hilbert space (describing the coupled system) are formed by taking the tensor product of the individual operators.

In QuTiP the function qutip.tensor.tensor is used to accomplish this task. This function takes as argument a collection:

>>> tensor(op1, op2, op3)

or a list:

>>> tensor([op1, op2, op3])

of state vectors or operators and returns a composite quantum object for the combined Hilbert space. The function accepts an arbitray number of states or operators as argument. The type returned quantum object is the same as that of the input(s).

For example, the state vector describing two qubits in their ground states is formed by taking the tensor product of the two single-qubit ground state vectors:

In [1]: tensor(basis(2, 0), basis(2, 0))
Out[1]: 
Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket
Qobj data =
[[1.]
 [0.]
 [0.]
 [0.]]

or equivalently using the list format:

In [2]: tensor([basis(2, 0), basis(2, 0)])
Out[2]: 
Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket
Qobj data =
[[1.]
 [0.]
 [0.]
 [0.]]

This is straightforward to generalize to more qubits by adding more component state vectors in the argument list to the qutip.tensor.tensor function, as illustrated in the following example:

In [3]: tensor((basis(2, 0) + basis(2, 1)).unit(),
   ...:       (basis(2, 0) + basis(2, 1)).unit(), basis(2, 0))
   ...: 
Out[3]: 
Quantum object: dims = [[2, 2, 2], [1, 1, 1]], shape = (8, 1), type = ket
Qobj data =
[[0.5]
 [0. ]
 [0.5]
 [0. ]
 [0.5]
 [0. ]
 [0.5]
 [0. ]]

This state is slightly more complicated, describing two qubits in a superposition between the up and down states, while the third qubit is in its ground state.

To construct operators that act on an extended Hilbert space of a combined system, we similarly pass a list of operators for each component system to the qutip.tensor.tensor function. For example, to form the operator that represents the simultaneous action of the \(\sigma_x\) operator on two qubits:

In [4]: tensor(sigmax(), sigmax())
Out[4]: 
Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
Qobj data =
[[0. 0. 0. 1.]
 [0. 0. 1. 0.]
 [0. 1. 0. 0.]
 [1. 0. 0. 0.]]

To create operators in a combined Hilbert space that only act only on a single component, we take the tensor product of the operator acting on the subspace of interest, with the identity operators corresponding to the components that are to be unchanged. For example, the operator that represents \(\sigma_z\) on the first qubit in a two-qubit system, while leaving the second qubit unaffected:

In [5]: tensor(sigmaz(), identity(2))
Out[5]: 
Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
Qobj data =
[[ 1.  0.  0.  0.]
 [ 0.  1.  0.  0.]
 [ 0.  0. -1.  0.]
 [ 0.  0.  0. -1.]]

Example: Constructing composite Hamiltonians¶

The qutip.tensor.tensor function is extensively used when constructing Hamiltonians for composite systems. Here we’ll look at some simple examples.

In [6]: H = tensor(sigmaz(), identity(2)) + tensor(identity(2),
   ...:           sigmaz()) + 0.05 * tensor(sigmax(), sigmax())
   ...: 

In [7]: H
Out[7]: 
Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
Qobj data =
[[ 2.    0.    0.    0.05]
 [ 0.    0.    0.05  0.  ]
 [ 0.    0.05  0.    0.  ]
 [ 0.05  0.    0.   -2.  ]]

In [8]: H = (tensor(sigmaz(), identity(2), identity(2)) +
   ...:     tensor(identity(2), sigmaz(), identity(2)) +
   ...:     tensor(identity(2), identity(2), sigmaz()) +
   ...:     0.5 * tensor(sigmax(), sigmax(), identity(2)) +
   ...:     0.25 * tensor(identity(2), sigmax(), sigmax()))
   ...: 

In [9]: H
Out[9]: 
Quantum object: dims = [[2, 2, 2], [2, 2, 2]], shape = (8, 8), type = oper, isherm = True
Qobj data =
[[ 3.    0.    0.    0.25  0.    0.    0.5   0.  ]
 [ 0.    1.    0.25  0.    0.    0.    0.    0.5 ]
 [ 0.    0.25  1.    0.    0.5   0.    0.    0.  ]
 [ 0.25  0.    0.   -1.    0.    0.5   0.    0.  ]
 [ 0.    0.    0.5   0.    1.    0.    0.    0.25]
 [ 0.    0.    0.    0.5   0.   -1.    0.25  0.  ]
 [ 0.5   0.    0.    0.    0.    0.25 -1.    0.  ]
 [ 0.    0.5   0.    0.    0.25  0.    0.   -3.  ]]

In [10]: N = 10

In [11]: omega_a = 1.0

In [12]: omega_c = 1.25

In [13]: g = 0.05

In [14]: a = tensor(identity(2), destroy(N))

In [15]: sm = tensor(destroy(2), identity(N))

In [16]: sz = tensor(sigmaz(), identity(N))

In [17]: H = 0.5 * omega_a * sz + omega_c * a.dag() * a + g * (a.dag() * sm + a * sm.dag())

Partial trace¶

The partial trace is an operation that reduces the dimension of a Hilbert space by eliminating some degrees of freedom by averaging (tracing). In this sense it is therefore the converse of the tensor product. It is useful when one is interested in only a part of a coupled quantum system. For open quantum systems, this typically involves tracing over the environment leaving only the system of interest. In QuTiP the class method qutip.Qobj.ptrace is used to take partial traces. qutip.Qobj.ptrace acts on the qutip.Qobj instance for which it is called, and it takes one argument sel, which is a list of integers that mark the component systems that should be kept. All other components are traced out.

For example, the density matrix describing a single qubit obtained from a coupled two-qubit system is obtained via:

In [18]:    psi = tensor(basis(2, 0), basis(2, 1))

In [19]:    psi.ptrace(0)
Out[19]: 
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[1. 0.]
 [0. 0.]]

In [20]:    psi.ptrace(1)
Out[20]: 
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[0. 0.]
 [0. 1.]]

Note that the partial trace always results in a density matrix (mixed state), regardless of whether the composite system is a pure state (described by a state vector) or a mixed state (described by a density matrix):

In [21]:    psi = tensor((basis(2, 0) + basis(2, 1)).unit(), basis(2, 0))

In [22]:    psi
Out[22]: 
Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket
Qobj data =
[[0.70710678]
 [0.        ]
 [0.70710678]
 [0.        ]]

In [23]:    psi.ptrace(0)
Out[23]: 
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[0.5 0.5]
 [0.5 0.5]]

In [24]:    rho = tensor(ket2dm((basis(2, 0) + basis(2, 1)).unit()), fock_dm(2, 0))

In [25]:    rho
Out[25]: 
Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True
Qobj data =
[[0.5 0.  0.5 0. ]
 [0.  0.  0.  0. ]
 [0.5 0.  0.5 0. ]
 [0.  0.  0.  0. ]]

In [26]:    rho.ptrace(0)
Out[26]: 
Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[0.5 0.5]
 [0.5 0.5]]

Superoperators and Tensor Manipulations¶

As described in Superoperators and Vectorized Operators, superoperators are operators that act on Liouville space, the vectorspace of linear operators. Superoperators can be represented using the isomorphism \(\mathrm{vec} : \mathcal{L}(\mathcal{H}) \to \mathcal{H} \otimes \mathcal{H}\) [Hav03], [Wat13]. To represent superoperators acting on \(\mathcal{L}(\mathcal{H}_1 \otimes \mathcal{H}_2)\) thus takes some tensor rearrangement to get the desired ordering \(\mathcal{H}_1 \otimes \mathcal{H}_2 \otimes \mathcal{H}_1 \otimes \mathcal{H}_2\).

In particular, this means that qutip.tensor does not act as one might expect on the results of qutip.to_super:

In [27]: A = qeye([2])

In [28]: B = qeye([3])

In [29]: to_super(tensor(A, B)).dims
Out[29]: [[[2, 3], [2, 3]], [[2, 3], [2, 3]]]

In [30]: tensor(to_super(A), to_super(B)).dims
Out[30]: [[[2], [2], [3], [3]], [[2], [2], [3], [3]]]

In the former case, the result correctly has four copies of the compound index with dims [2, 3]. In the latter case, however, each of the Hilbert space indices is listed independently and in the wrong order.

The qutip.super_tensor function performs the needed rearrangement, providing the most direct analog to qutip.tensor on the underlying Hilbert space. In particular, for any two type="oper" Qobjs A and B, to_super(tensor(A, B)) == super_tensor(to_super(A), to_super(B)) and operator_to_vector(tensor(A, B)) == super_tensor(operator_to_vector(A), operator_to_vector(B)). Returning to the previous example:

In [31]: super_tensor(to_super(A), to_super(B)).dims
Out[31]: [[[2, 3], [2, 3]], [[2, 3], [2, 3]]]

The qutip.composite function automatically switches between qutip.tensor and qutip.super_tensor based on the type of its arguments, such that composite(A, B) returns an appropriate Qobj to represent the composition of two systems.

In [32]: composite(A, B).dims
Out[32]: [[2, 3], [2, 3]]

In [33]: composite(to_super(A), to_super(B)).dims
Out[33]: [[[2, 3], [2, 3]], [[2, 3], [2, 3]]]

QuTiP also allows more general tensor manipulations that are useful for converting between superoperator representations [WBC11]. In particular, the tensor_contract function allows for contracting one or more pairs of indices. As detailed in the channel contraction tutorial, this can be used to find superoperators that represent partial trace maps. Using this functionality, we can construct some quite exotic maps, such as a map from \(3 \times 3\) operators to \(2 \times 2\) operators:

In [34]: tensor_contract(composite(to_super(A), to_super(B)), (1, 3), (4, 6)).dims
Out[34]: [[[2], [2]], [[3], [3]]]