Quantum process tomography

Quantum process tomography (QPT) is a useful technique for characterizing experimental implementations of quantum gates involving a small number of qubits. Recently, QPT has been used in a number of beautiful experiments with for example superconducting qubits [Bialczak2010][Dewes2012], nuclear spins [Steger2012], NV-centers [Fuchs2012], ion traps [Riebe2006], just to mention a few.

QPT can also be a useful theoretical tool, which can give insight in how a gate transforms states, and be used to study how noise or other imperfections deteriorate a gate. Whereas a fidelity or distance measure can give a single number that indicates how far from ideal a gate is, a quantum process tomography analysis can give detailed information about exactly what kind of errors various imperfections introduce.

In QuTiP, we have recently added support for computing and visualizing quantum process tomography matrices. Here we give a simple example that demonstrates these new functions for a few common one, two and three qubit gates. The gates in this example are all available in QuTiP's gate module.

"""
Plot the process tomography matrices for some 1, 2, and 3-qubit qubit gates.
"""
from qutip import *

gate_list = [['C-NOT', cnot()],
             ['SWAP', swap()],
             ['$i$SWAP', iswap()],
             ['$\sqrt{i\mathrm{SWAP}}$', sqrtiswap()],
             ['S-NOT', snot()],
             ['$\pi/2$ phase gate', phasegate(pi/2)],
             ['Toffoli', toffoli()],
             ['Fredkin', fredkin()]]

# loop though the gate list and plot the gates
for gate in gate_list:

 name  = gate[0]
 U_psi = gate[1]

 N = len(U_psi.dims[0]) # number of qubits

 # create a superoperator for the density matrix
 # transformation rho = U_psi * rho_0 * U_psi.dag()
 U_rho = spre(U_psi) * spost(U_psi.dag())

 # operator basis for the process tomography
 op_basis = [[qeye(2), sigmax(), sigmay(), sigmaz()] for i in range(N)]

 # labels for operator basis
 op_label = [["i", "x", "y", "z"] for i in range(N)]

 # calculate the chi matrix
 chi = qpt(U_rho, op_basis)

 # visualize the chi matrix
 qpt_plot_combined(chi, op_label, name)

show()

This program produces the following figures:

These figures are all examples of ideal quantum gates. To illustrate how the quality of a gate can deteriorate due to interaction between the qubits and their environment, let's consider a simple but experimentally realizable implementation of the two-qubit iSWAP gate. One way of implementing this gate is to produce the interaction Hamiltonian $H = g\left(\sigma_x\otimes\sigma_x + \sigma_y\otimes\sigma_y\right)$, and turn on this interaction for a period of time $T=\pi/4g$ [Schuch2003] [Steffen2006]. In the absence of noise and other imperfections, this procedure results in an ideal iSWAP gate to be applied to the two qubits. With QuTiP, we can easily simulate this process, and additionally study how and to what extent the gate deteriorates due to dissipation, imperfect switching of the interaction, detuning, etc. The following code demonstrates how to do this kind of analysis using QuTiP.

"""
Plot the process tomography matrix for a 2-qubit iSWAP gate.
"""
from qutip import *

g = 1.0 * 2 * pi # coupling strength
g1 = 0.75        # relaxation rate
g2 = 0.25        # dephasing rate
n_th = 1.5       # bath temperature

T = pi/(4*g)
H = g * (tensor(sigmax(), sigmax()) + tensor(sigmay(), sigmay()))
psi0 = tensor(basis(2,1), basis(2,0))

c_ops = []
# qubit 1 collapse operators
sm1 = tensor(sigmam(), qeye(2))
sz1 = tensor(sigmaz(), qeye(2))
c_ops.append(sqrt(g1 * (1+n_th)) * sm1)
c_ops.append(sqrt(g1 * n_th) * sm1.dag())
c_ops.append(sqrt(g2) * sz1)
# qubit 2 collapse operators
sm2 = tensor(qeye(2), sigmam())
sz2 = tensor(qeye(2), sigmaz())
c_ops.append(sqrt(g1 * (1+n_th)) * sm2)
c_ops.append(sqrt(g1 * n_th) * sm2.dag())
c_ops.append(sqrt(g2) * sz2)

# process tomography basis
op_basis = [[qeye(2), sigmax(), sigmay(), sigmaz()]] * 2
op_label = [["i", "x", "y", "z"]] * 2

# dissipative gate
U_diss = propagator(H, T, c_ops)
chi = qpt(U_diss, op_basis)
qpt_plot_combined(chi, op_label)

# ideal gate 
U_psi = (-1j * H * T).expm()
U_ideal = spre(U_psi) * spost(U_psi.dag())
chi = qpt(U_ideal, op_basis)
qpt_plot_combined(chi, op_label)

show()

The code produces the following figures (the ideal iSWAP gate to the left, the dissipative one to the right):

Note: the quantum process tomography functions were added after the 2.0 release, so to run these examples you need to get the latest development version of QuTiP from the SVN repository. Alternatively, if you prefer to stay on the official release (recommended for stability), you can download this file and put it in the same directory as the example code and add "from qpt import *" at the top of the example files (just below the line "from qutip import *").

QuTiP talk at scipy2012

Slides for a talk about QuTiP at the conference Scientific Computing with Python, at Austin 2012, can be downloaded from here.

How large quantum systems can QuTiP handle?

People are often curious to know how many qubits can be simulated on a modern computer, using for example QuTiP. The primary focus of QuTiP is not to manage as large quantum systems as possible, but we do try to make QuTiP efficient, and it is useful to have some feeling for how large systems we can manage.
The answer to the question posed in the title depends strongly on the type of the problem. For example, whether the Hamiltonian is constant or time-dependent, the system is closed or open, states are pure or mixed. Ultimately it is the RAM memory that limits the size of the systems we can manage, but in some cases you are likely to run out of patience long before the RAM memory is exhausted (in particular for time-dependent systems).
Computationally, the simplest problem that QuTiP will typically be used for is the evolution of a closed quantum system with a constant Hamiltonian. To explore how large quantum systems of this kind we can simulate with QuTiP, let’s consider a chain of spins coupled to their nearest neighbors through a sigma-x interaction. In QuTiP it is easy to program this model in a way that parameterize the number of spins in the chain, allowing us to simulate systems with increasing number of spins until we run out of RAM memory (resulting in a MemoryError exception being raised by python) or until we run out of patience. You can find the program code here:

"""
QuTiP example: Dynamics of a spin-1/2 chain
"""
from qutip import *
from pylab import *
rcParams['text.usetex'] = True
rcParams['font.size'] = 16
rcParams['font.family'] = 'serif'

def system_ops(N, h, J, gamma):

    si = qeye(2); sx = sigmax(); sz = sigmaz()

    sx_list = []; sz_list = []
    for n in range(N):
        op_list = [si for m in range(N)]
        op_list[n] = sx
        sx_list.append(tensor(op_list))
        op_list[n] = sz
        sz_list.append(tensor(op_list))
        
    # construct the hamiltonian
    H = 0
    for n in range(N):
        H += - 0.5 * h[n] * sz_list[n]                # energy splitting terms
    for n in range(N-1):
        H += - 0.5 * J[n] * sx_list[n] * sx_list[n+1] # interaction terms

    # collapse operators for spin dephasing
    c_op_list = []
    for n in range(N):
        if gamma[n] > 0.0:
            c_op_list.append(sqrt(gamma[n]) * sz_list[n])

    # intital state, first and last spin in state |1>, the rest in state |0>
    psi_list = [basis(2,0) for n in range(N-2)]
    psi_list.insert(0, basis(2,1)) # first
    psi_list.append(basis(2,1))    # last
    psi0 = tensor(psi_list)

    return H, c_op_list, sz_list, psi0

def evolve(N, h, J, tlist, gamma, solver):

    H, c_ops, e_ops, psi0 = system_ops(N, h, J, gamma)
    output = solver(H, psi0, tlist, c_ops, e_ops)
    return (1 - array(output.expect))/2.0 # <sigma_z> to excitation prob.

# ---
solver = mesolve              # select solver mesolve or mcsolve
N = 16                        # number of spins
h =  1.0 * 2 * pi * ones(N)   # energy splittings
J = 0.05 * 2 * pi * ones(N)   # coupling
gamma = 0.0 * ones(N)         # dephasing rate
tlist = linspace(0, 150, 150)

sn_expt = evolve(N, h, J, tlist, gamma, solver)

# plot excitation probabilities for each spin
yvals = []
ylabels = []
for n in range(N):
    x = 0.1 # spacing between stacked plots
    p = real(sn_expt[n]) * (1-2*x)
    base = n + zeros(shape(p)) + 0.5 + x
    plot(tlist, base, 'k', lw=1)        
    plot(tlist, base+(1-2*x), 'k', lw=1)        
    fill_between(tlist, base, base + p, color='b', alpha=0.75) 
    fill_between(tlist, base + p, base + (1-2*x), color='g', alpha=0.2)
    yvals.append(n + 1)
    ylabels.append("Qubit %d" % (n+1,))

autoscale(tight=True)
xlabel(r'Time [ns]',fontsize=14)
title(r'Occupation probabilities of spins in a chain', fontsize=16)
yticks(yvals, ylabels, fontsize=12)
savefig('spinchain-%d.pdf' % N)

In this simulation, we initially put excitations at the spins at the end of the chain and evolve the state vector. The result shows how the excitations propagate in the spin chain and eventually disperse.

My patience ran out after 21 spins. This corresponds to a very large quantum state space, 2’097’152 states to be exact, which is a quite impressive number. But admittedly, this simulation took a long time to run: About 7 hours on my laptop, and the memory usage peaked at about 4.3 Gb of RAM.
Open quantum systems are significantly more computationally demanding to evolve. First, we need to use the density matrix rather than a state vector to describe the system, which immediately squares the problem size. Second, the equations of motion (Master equation) are more complicated, including for example superoperators, requiring us to rearrange the equations of motion to obtain an ODE on standard form (of squared system size). Alternatively, we can use a stochastic equation of motion, such as the quantum trajectory Monte-Carlo method, and circumvent the requirement to square the system size, but this also comes at a steep price as we then need to average a large number of trajectories with the original system size. Regardless if we use a master equation or the Monte-Carlo quantum trajectory method, we cannot hope to evolve as large systems as we could for a closed quantum system.

Including dissipation in the example program above (by selecting gamma=0.0025), I could simulate systems up to 12 spins before the RAM memory was exhausted (on a 16Gb desktop). The 12 spin simulation took a weekend to run, and the peak memory usage was about 20 Gb (temporarily swapping 4 Gb). On this computer, further increasing the system size would result in memory error when using the master equation solver. Using the Monto-Carlo solver (by setting solver=mcsolve in the code) I could simulate open systems up to 16 spins in a similar time period (a few days).
If we include a time-dependence in for example the Hamiltonian, the memory foot-print is not significantly altered, but the problem becomes more CPU demanding, increasing the simulation time by sometimes quite a lot. However, at least in principle, we can simulate time-dependent quantum systems of similar size as time-independent problems.