Manipulating States and Operators

Introduction

In the previous guide section Basic Operations on Quantum Objects, we saw how to create states and operators, using the functions built into QuantumToolbox. In this portion of the guide, we will look at performing basic operations with states and operators. For more detailed demonstrations on how to use and manipulate these objects, see the examples given in the tutorial section.

State Vectors (kets or bras)

Here we begin by creating a Fock basis (or fock) vacuum state vector $|0⟩$ with in a Hilbert space with 5 number states, from 0 to 4:

vac = basis(5, 0)

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

and then create a lowering operator $\hat{a}$ corresponding to 5 number states using the destroy function:

a = destroy(5)

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=false
5×5 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
     ⋅      1.0+0.0im          ⋅              ⋅          ⋅    
     ⋅          ⋅      1.41421+0.0im          ⋅          ⋅    
     ⋅          ⋅              ⋅      1.73205+0.0im      ⋅    
     ⋅          ⋅              ⋅              ⋅      2.0+0.0im
     ⋅          ⋅              ⋅              ⋅          ⋅

Now lets apply the lowering operator $\hat{a}$ to our vacuum state vac:

a * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

We see that, as expected, the vacuum is transformed to the zero vector. A more interesting example comes from using the adjoint of the lowering operator $\hat{a}$, the raising operator ${\hat{a}}^{†}$:

a' * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

The raising operator has in indeed raised the state vac from the vacuum to the $|1⟩$ state. Instead of using the adjoint method to raise the state, we could have also used the built-in create function to make a raising operator:

ad = create(5)
ad * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

which does the same thing. We can raise the vacuum state more than once by successively apply the raising operator:

ad * ad * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
                0.0 + 0.0im
                0.0 + 0.0im
 1.4142135623730951 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

or just taking the square of the raising operator ${\left({\hat{a}}^{†}\right)}^2$:

ad^2 * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
                0.0 + 0.0im
                0.0 + 0.0im
 1.4142135623730951 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

Applying the raising operator twice gives the expected $\sqrt{n+1}$ dependence. We can use the product of ${\hat{a}}^{†}\hat{a}$ to also apply the number operator to the state vector vac:

ad * a * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

or on the $|1⟩$ state:

ad * a * (ad * vac)

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

or on the $|2⟩$ state:

ad * a * (ad^2 * vac)

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
                0.0 + 0.0im
                0.0 + 0.0im
 2.8284271247461907 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

Notice how in this last example, application of the number operator does not give the expected value $n=2$, but rather $2\sqrt{2}$. This is because this last state is not normalized to unity as ${\hat{a}}^{†}|n⟩=\sqrt{n+1}|n+1⟩$. Therefore, we should normalize (or use unit) our vector first:

ad * a * normalize(ad^2 * vac)

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
                0.0 + 0.0im
                0.0 + 0.0im
 2.0000000000000004 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

Since we are giving a demonstration of using states and operators, we have done a lot more work than we should have. For example, we do not need to operate on the vacuum state to generate a higher number Fock state. Instead we can use the basis (or fock) function to directly obtain the required state:

ket = basis(5, 2)

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

Notice how it is automatically normalized. We can also use the built in number operator num:

n = num(5)

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=true
5×5 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 5 stored entries:
 0.0+0.0im      ⋅          ⋅          ⋅          ⋅    
     ⋅      1.0+0.0im      ⋅          ⋅          ⋅    
     ⋅          ⋅      2.0+0.0im      ⋅          ⋅    
     ⋅          ⋅          ⋅      3.0+0.0im      ⋅    
     ⋅          ⋅          ⋅          ⋅      4.0+0.0im

Therefore, instead of ad * a * normalize(ad^2 * vac), we have:

n * ket

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 2.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

We can also create superpositions of states:

ket = normalize(basis(5, 0) + basis(5, 1))

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
 0.7071067811865475 + 0.0im
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

where we have used the normalize function again to normalize the state. Apply the number operator again:

n * ket

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
                0.0 + 0.0im
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im

We can also create coherent states and squeezed states by applying the displace and squeeze functions to the vacuum state:

vac = basis(5, 0)

d = displace(5, 1im)

s = squeeze(5, 0.25 + 0.25im)

d * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
  0.6065568176126114 + 0.0im
                 0.0 + 0.6062813254779008im
 -0.4303873979781239 + 0.0im
                 0.0 - 0.24104350624628343im
 0.14552146626026782 + 0.0im

d * s * vac

Quantum Object:   type=Ket()   dims=[5]   size=(5,)
5-element Vector{ComplexF64}:
   0.6589378628979528 + 0.08139380927428255im
  0.10779461991734264 + 0.5157973476443225im
   -0.375672173778824 - 0.013268528813115979im
 -0.02688063342547209 - 0.2382877475293735im
    0.263528135717665 + 0.11512177609766486im

Of course, displacing the vacuum gives a coherent state, which can also be generated using the built in coherent function.

Density matrices

One of the main purpose of QuantumToolbox is to explore the dynamics of open quantum systems, where the most general state of a system is no longer a state vector, but rather a density matrix. Since operations on density matrices operate identically to those of vectors, we will just briefly highlight creating and using these structures.

The simplest density matrix is created by forming the outer-product $|\psi ⟩⟨\psi |$ of a ket vector:

ket = basis(5, 2)
ket * ket'

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=true
5×5 Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

A similar task can also be accomplished via the fock_dm or ket2dm functions:

fock_dm(5, 2)

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=true
5×5 Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

ket2dm(ket)

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=true
5×5 Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

If we want to create a density matrix with equal classical probability of being found in the $|2⟩$ or $|4⟩$ number states, we can do the following:

0.5 * fock_dm(5, 2) + 0.5 * fock_dm(5, 4) # with fock_dm
0.5 * ket2dm(basis(5, 2)) + 0.5 * ket2dm(basis(5, 4)) # with ket2dm

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=true
5×5 Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.5+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.5+0.0im

There are also several other built-in functions for creating predefined density matrices, for example coherent_dm and thermal_dm which create coherent state and thermal state density matrices, respectively.

coherent_dm(5, 1.25)

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=true
5×5 Matrix{ComplexF64}:
 0.209807+0.0im  0.261411+0.0im  …   0.155726+0.0im   0.133908+0.0im
 0.261411+0.0im  0.325707+0.0im      0.194028+0.0im   0.166843+0.0im
 0.235097+0.0im  0.292921+0.0im      0.174497+0.0im   0.150049+0.0im
 0.155726+0.0im  0.194028+0.0im      0.115585+0.0im  0.0993908+0.0im
 0.133908+0.0im  0.166843+0.0im     0.0993908+0.0im  0.0854655+0.0im

thermal_dm(5, 1.25)

Quantum Object:   type=Operator()   dims=[5]   size=(5, 5)   ishermitian=true
5×5 Matrix{Float64}:
 0.46928  0.0       0.0       0.0        0.0
 0.0      0.260711  0.0       0.0        0.0
 0.0      0.0       0.144839  0.0        0.0
 0.0      0.0       0.0       0.0804663  0.0
 0.0      0.0       0.0       0.0        0.0447035

QuantumToolbox also provides a set of distance metrics for determining how close two density matrix distributions are to each other. For example, fidelity, and trace distance (tracedist) are included. For more metric functions, see section Entropy and Metrics in the API page.

x = coherent_dm(5, 1.25)

y = coherent_dm(5, 1.25im)

z = thermal_dm(5, 0.125)

fidelity(x, y)

0.2125206772965313

Note that the definition of fidelity here is from Nielsen and Chuang (2012). It is the square root of the fidelity defined in Jozsa (1994). We also know that for two pure states, the trace distance ($T$) and the fidelity ($F$) are related by $T=\sqrt{1-F^2}$:

tracedist(x, y) ≈ sqrt(1 - (fidelity(x, y))^2)

true

For a pure state and a mixed state, $1-F\le T$ which can also be verified:

1 - fidelity(x, z) < tracedist(x, z)

true

Two-level systems (Qubits)

Having spent a fair amount of time on basis states that represent harmonic oscillator states, we now move on to qubit, or two-level quantum systems (for example a spin-$1/2$). To create a state vector corresponding to a qubit system, we use the same basis, or fock, function with only two levels:

spin = basis(2, 0)

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im

Now at this point one may ask how this state is different than that of a harmonic oscillator in the vacuum state truncated to two energy levels?

vac = basis(2, 0)

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im

At this stage, there is no difference. This should not be surprising as we called the exact same function twice. The difference between the two comes from the action of the spin operators sigmax, sigmay, sigmaz, sigmap, and sigmam on these two-level states. For example, if vac corresponds to the vacuum state of a harmonic oscillator, then, as we have already seen, we can use the raising operator (create) to get the $|1⟩$ state:

create(2) * vac

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im

For a spin system, the operator analogous to the raising operator is the ${\hat{\sigma }}_{+}$ operator sigmap. Applying on the spin state gives:

sigmap() * spin

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im

Now we see the difference! The sigmap operator acting on the spin state returns the zero vector. Why is this? To see what happened, let us use the ${\hat{\sigma }}_z$ (sigmaz) operator:

sigmaz()

Quantum Object:   type=Operator()   dims=[2]   size=(2, 2)   ishermitian=true
2×2 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
 1.0+0.0im       ⋅    
     ⋅      -1.0+0.0im

sigmaz() * spin

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im

spin2 = basis(2, 1)

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im

sigmaz() * spin2

Quantum Object:   type=Ket()   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
  0.0 + 0.0im
 -1.0 + 0.0im

The answer is now apparent. Since the QuantumToolbox sigmaz function uses the standard $Z$-basis representation of the ${\hat{\sigma }}_z$ spin operator, the spin state corresponds to the $|\uparrow ⟩$ state of a two-level spin system while spin2 gives the $|\downarrow ⟩$ state. Therefore, in our previous example sigmap() * spin, we raised the qubit state out of the truncated two-level Hilbert space resulting in the zero state.

While at first glance this convention might seem somewhat odd, it is in fact quite handy. For one, the spin operators remain in the conventional form. Second, this corresponds nicely with the quantum information definitions of qubit states, where the excited $|\uparrow ⟩$ state is label as $|0⟩$, and the $|\downarrow ⟩$ state by $|1⟩$.

If one wants to create spin operators for higher spin systems, then the jmat function comes in handy.

Expectation values

Some of the most important information about quantum systems comes from calculating the expectation value of operators, both Hermitian and non-Hermitian, as the state or density matrix of the system varies in time. Therefore, in this section we demonstrate the use of the expect function. To begin:

vac = basis(5, 0)

one = basis(5, 1)

c = create(5)

N = num(5)

coh = coherent_dm(5, 1.0im)

cat = normalize(basis(5, 4) + 1.0im * basis(5, 3))

println(expect(N, vac) ≈ 0)
println(expect(N, one) ≈ 1)
println(expect(N, coh) ≈ 0.9970555745806597)
println(expect(c, cat) ≈ 1im)

true
true
true
true

The expect function also accepts lists or arrays of state vectors or density matrices for the second input:

states = [normalize(c^k * vac) for k in 0:4]

expect(N, states)

5-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im
 2.0 + 0.0im
 3.0 + 0.0im
 4.0 + 0.0im

cat_list = [normalize(basis(5, 4) + x * basis(5, 3)) for x in [0, 1.0im, -1.0, -1.0im]]

expect(c, cat_list)

4-element Vector{ComplexF64}:
                 0.0 + 0.0im
                 0.0 + 0.9999999999999998im
 -0.9999999999999998 + 0.0im
                 0.0 - 0.9999999999999998im

Notice how in this last example, all of the return values are complex numbers. This is because the expect function looks to see whether the operator is Hermitian or not. If the operator is Hermitian, then the output will always be real. In the case of non-Hermitian operators, the return values may be complex. Therefore, the expect function will return an array of complex values for non-Hermitian operators when the input is a list/array of states or density matrices.

Of course, the expect function works for spin states and operators:

up = basis(2, 0)

dn = basis(2, 1)

println(expect(sigmaz(), up) ≈ 1)
println(expect(sigmaz(), dn) ≈ -1)

true
true

as well as the composite objects discussed in the next section Tensor Products and Partial Traces:

spin1 = basis(2, 0)

spin2 = basis(2, 1)

two_spins = tensor(spin1, spin2)

sz1 = tensor(sigmaz(), qeye(2))

sz2 = tensor(qeye(2), sigmaz())

println(expect(sz1, two_spins) ≈ 1)
println(expect(sz2, two_spins) ≈ -1)

true
true

Superoperators and Vectorized Operators

In addition to state vectors and density operators, QuantumToolbox allows for representing maps that act linearly on density operators using the Liouville supermatrix formalisms.

This support is based on the correspondence between linear operators acting on a Hilbert space, and vectors in two copies of that Hilbert space (which is also called the Fock-Liouville space),

$$\text{vec}:\mathcal{L}(\mathcal{H})\to \mathcal{H}\otimes \mathcal{H}.$$

Therefore, a given density matrix $\hat{\rho }$ can then be vectorized, denoted as

$$|\hat{\rho }⟩⟩=\text{vec}(\hat{\rho }).$$

QuantumToolbox uses the column-stacking convention for the isomorphism between $\mathcal{L}(\mathcal{H})$ and $\mathcal{H}\otimes \mathcal{H}$. This isomorphism is implemented by the functions mat2vec (or operator_to_vector) and vec2mat (or vector_to_operator):

rho = Qobj([1 2; 3 4])

Quantum Object:   type=Operator()   dims=[2]   size=(2, 2)   ishermitian=false
2×2 Matrix{Int64}:
 1  2
 3  4

vec_rho = mat2vec(rho)

Quantum Object:   type=OperatorKet()   dims=[2]   size=(4,)
4-element Vector{Int64}:
 1
 3
 2
 4

rho2 = vec2mat(vec_rho)

Quantum Object:   type=Operator()   dims=[2]   size=(2, 2)   ishermitian=false
2×2 Matrix{Int64}:
 1  2
 3  4

The QuantumObject.type attribute indicates whether a quantum object is a vector corresponding to an OperatorKet, or its Hermitian conjugate OperatorBra. One can also use isoper, isoperket, and isoperbra to check the type:

println(isoper(vec_rho))
println(isoperket(vec_rho))
println(isoperbra(vec_rho))
println(isoper(vec_rho'))
println(isoperket(vec_rho'))
println(isoperbra(vec_rho'))

false
true
false
false
false
true

Because Julia is a column-oriented languages (like Fortran and MATLAB), in QuantumToolbox, we define the spre (left), spost (right), and sprepost (left-and-right) multiplication superoperators as follows:

$$\begin{array}{rl}\hat{A}\hat{\rho }& \to \text{spre}(\hat{A})\ast \text{vec}(\hat{\rho })=\hat{\mathbb{1}}\otimes \hat{A}|\hat{\rho }⟩⟩,\\ \hat{\rho }\hat{B}& \to \text{spost}(\hat{B})\ast \text{vec}(\hat{\rho })={\hat{B}}^T\otimes \hat{\mathbb{1}}|\hat{\rho }⟩⟩,\\ \hat{A}\hat{\rho }\hat{B}& \to \text{sprepost}(\hat{A},\hat{B})\ast \text{vec}(\hat{\rho })={\hat{B}}^T\otimes \hat{A}|\hat{\rho }⟩⟩,\end{array}$$

where $\hat{\mathbb{1}}$ represents the identity operator with Hilbert space dimension equal to $\hat{\rho }$.

A = Qobj([1 2; 3 4])
S_A = spre(A)

Quantum Object:   type=SuperOperator()   dims=[2]   size=(4, 4)
4×4 SparseArrays.SparseMatrixCSC{Int64, Int64} with 8 stored entries:
 1  2  ⋅  ⋅
 3  4  ⋅  ⋅
 ⋅  ⋅  1  2
 ⋅  ⋅  3  4

B = Qobj([5 6; 7 8])
S_B = spost(B)

Quantum Object:   type=SuperOperator()   dims=[2]   size=(4, 4)
4×4 SparseArrays.SparseMatrixCSC{Int64, Int64} with 8 stored entries:
 5  ⋅  7  ⋅
 ⋅  5  ⋅  7
 6  ⋅  8  ⋅
 ⋅  6  ⋅  8

S_AB = sprepost(A, B)

Quantum Object:   type=SuperOperator()   dims=[2]   size=(4, 4)
4×4 SparseArrays.SparseMatrixCSC{Int64, Int64} with 16 stored entries:
  5  10   7  14
 15  20  21  28
  6  12   8  16
 18  24  24  32

S_AB ≈ S_A * S_B ≈ S_B * S_A

true

One can also use issuper to check the type:

println(isoper(S_AB))
println(issuper(S_AB))

false
true

With the above definitions, the following equalities hold in Julia:

$$\text{vec}(\hat{A}\hat{\rho }\hat{B})=\text{spre}(\hat{A})\ast \text{spost}(\hat{B})\ast \text{vec}(\hat{\rho })=\text{sprepost}(\hat{A},\hat{B})\ast \text{vec}(\hat{\rho })\mathrm{\forall }\hat{A},\hat{B},\hat{\rho }$$

N  = 10
A = Qobj(rand(ComplexF64, N, N))
B = Qobj(rand(ComplexF64, N, N))
ρ = rand_dm(N) # random density matrix
mat2vec(A * ρ * B) ≈ spre(A) * spost(B) * mat2vec(ρ) ≈ sprepost(A, B) * mat2vec(ρ)

true

In addition, dynamical generators on this extended space, often called Liouvillian superoperators, can be created using the liouvillian function. Each of these takes a Hamiltonian along with a list of collapse operators, and returns a QuantumObject of type SuperOperator that can be exponentiated to find the superoperator for that evolution.

H = 10 * sigmaz()

c = destroy(2)

L = liouvillian(H, [c])

Quantum Object:   type=SuperOperator()   dims=[2]   size=(4, 4)
4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
     ⋅           ⋅            ⋅        1.0+0.0im
     ⋅      -0.5+20.0im       ⋅            ⋅    
     ⋅           ⋅       -0.5-20.0im       ⋅    
     ⋅           ⋅            ⋅       -1.0+0.0im

t = 0.8
exp(L * t)

Quantum Object:   type=SuperOperator()   dims=[2]   size=(4, 4)
4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 5 stored entries:
 1.0+0.0im            ⋅                     ⋅           0.550671+0.0im
     ⋅      -0.641938-0.192987im            ⋅                    ⋅    
     ⋅                ⋅           -0.641938+0.192987im           ⋅    
     ⋅                ⋅                     ⋅           0.449329+0.0im

See the section Lindblad Master Equation Solver for more details.