Visualization is often an important complement to a simulation of a quantummechanical system. The first method of visualization that come to mind might beto plot the expectation values of a few selected operators. But on top of that,it can often be instructive to visualize for example the state vectors ordensity matices that describe the state of the system, or how the state istransformed as a function of time (see process tomography below). In thissection we demonstrate how QuTiP and matplotlib can be used to perform a fewtypes of visualizations that often can provide additional understanding ofquantum system.
Fock-basis probability distribution¶
In quantum mechanics probability distributions plays an important role, and asin statistics, the expectation values computed from a probability distributiondoes not reveal the full story. For example, consider an quantum harmonicoscillator mode with Hamiltonian \(H = \hbar\omega a^\dagger a\), which isin a state described by its density matrix \(\rho\), and which on averageis occupied by two photons, \(\mathrm{Tr}[\rho a^\dagger a] = 2\). Giventhis information we cannot say whether the oscillator is in a Fock state,a thermal state, a coherent state, etc. By visualizing the photon distributionin the Fock state basis important clues about the underlying state can beobtained.
One convenient way to visualize a probability distribution is to use histograms.Consider the following histogram visualization of the number-basis probabilitydistribution, which can be obtained from the diagonal of the density matrix,for a few possible oscillator states with on average occupation of two photons.
First we generate the density matrices for the coherent, thermal and fock states.
In [1]: N = 20In [2]: rho_coherent = coherent_dm(N, np.sqrt(2))In [3]: rho_thermal = thermal_dm(N, 2)In [4]: rho_fock = fock_dm(N, 2)
Next, we plot histograms of the diagonals of the density matrices:
In [5]: fig, axes = plt.subplots(1, 3, figsize=(12,3))In [6]: bar0 = axes[0].bar(np.arange(0, N)-.5, rho_coherent.diag())In [7]: lbl0 = axes[0].set_title("Coherent state")In [8]: lim0 = axes[0].set_xlim([-.5, N])In [9]: bar1 = axes[1].bar(np.arange(0, N)-.5, rho_thermal.diag())In [10]: lbl1 = axes[1].set_title("Thermal state")In [11]: lim1 = axes[1].set_xlim([-.5, N])In [12]: bar2 = axes[2].bar(np.arange(0, N)-.5, rho_fock.diag())In [13]: lbl2 = axes[2].set_title("Fock state")In [14]: lim2 = axes[2].set_xlim([-.5, N])In [15]: plt.show()
All these states correspond to an average of two photons, but by visualizingthe photon distribution in Fock basis the differences between these states areeasily appreciated.
One frequently need to visualize the Fock-distribution in the way describedabove, so QuTiP provides a convenience function for doing this, seequtip.visualization.plot_fock_distribution, and the following example:
In [16]: fig, axes = plt.subplots(1, 3, figsize=(12,3))In [17]: plot_fock_distribution(rho_coherent, fig=fig, ax=axes[0], title="Coherent state");In [18]: plot_fock_distribution(rho_thermal, fig=fig, ax=axes[1], title="Thermal state");In [19]: plot_fock_distribution(rho_fock, fig=fig, ax=axes[2], title="Fock state");In [20]: fig.tight_layout()In [21]: plt.show()
Quasi-probability distributions¶
The probability distribution in the number (Fock) basis only describes theoccupation probabilities for a discrete set of states. A more completephase-space probability-distribution-like function for harmonic modes arethe Wigner and Husumi Q-functions, which are full descriptions of thequantum state (equivalent to the density matrix). These are calledquasi-distribution functions because unlike real probability distributionfunctions they can for example be negative. In addition to being more complete descriptionsof a state (compared to only the occupation probabilities plotted above),these distributions are also great for demonstrating if a quantum state isquantum mechanical, since for example a negative Wigner functionis a definite indicator that a state is distinctly nonclassical.
Wigner function¶
In QuTiP, the Wigner function for a harmonic mode can be calculated with thefunction qutip.wigner.wigner. It takes a ket or a density matrix asinput, together with arrays that define the ranges of the phase-spacecoordinates (in the x-y plane). In the following example the Wigner functionsare calculated and plotted for the same three states as in the previous section.
In [22]: xvec = np.linspace(-5,5,200)In [23]: W_coherent = wigner(rho_coherent, xvec, xvec)In [24]: W_thermal = wigner(rho_thermal, xvec, xvec)In [25]: W_fock = wigner(rho_fock, xvec, xvec)In [26]: # plot the resultsIn [27]: fig, axes = plt.subplots(1, 3, figsize=(12,3))In [28]: cont0 = axes[0].contourf(xvec, xvec, W_coherent, 100)In [29]: lbl0 = axes[0].set_title("Coherent state")In [30]: cont1 = axes[1].contourf(xvec, xvec, W_thermal, 100)In [31]: lbl1 = axes[1].set_title("Thermal state")In [32]: cont0 = axes[2].contourf(xvec, xvec, W_fock, 100)In [33]: lbl2 = axes[2].set_title("Fock state")In [34]: plt.show()
Custom Color Maps¶
The main objective when plotting a Wigner function is to demonstrate that the underlyingstate is nonclassical, as indicated by negative values in the Wigner function. Therefore,making these negative values stand out in a figure is helpful for both analysis and publicationpurposes. Unfortunately, all of the color schemes used in Matplotlib (or any other plotting software)are linear colormaps where small negative values tend to be near the same color as the zero values, andare thus hidden. To fix this dilemma, QuTiP includes a nonlinear colormap function qutip.visualization.wigner_cmapthat colors all negative values differently than positive or zero values. Below is a demonstration of how to usethis function in your Wigner figures:
In [35]: import matplotlib as mplIn [36]: from matplotlib import cmIn [37]: psi = (basis(10, 0) + basis(10, 3) + basis(10, 9)).unit()In [38]: xvec = np.linspace(-5, 5, 500)In [39]: W = wigner(psi, xvec, xvec)In [40]: wmap = wigner_cmap(W) # Generate Wigner colormapIn [41]: nrm = mpl.colors.Normalize(-W.max(), W.max())In [42]: fig, axes = plt.subplots(1, 2, figsize=(10, 4))In [43]: plt1 = axes[0].contourf(xvec, xvec, W, 100, cmap=cm.RdBu, norm=nrm)In [44]: axes[0].set_title("Standard Colormap");In [45]: cb1 = fig.colorbar(plt1, ax=axes[0])In [46]: plt2 = axes[1].contourf(xvec, xvec, W, 100, cmap=wmap) # Apply Wigner colormapIn [47]: axes[1].set_title("Wigner Colormap");In [48]: cb2 = fig.colorbar(plt2, ax=axes[1])In [49]: fig.tight_layout()In [50]: plt.show()
Husimi Q-function¶
The Husimi Q function is, like the Wigner function, a quasiprobabilitydistribution for harmonic modes. It is defined as
\[\begin{split}Q(\alpha) = \frac{1}{\pi}\left<\alpha|\rho|\alpha\right>\end{split}\]
where \(\left|\alpha\right>\) is a coherent state and\(\alpha = x + iy\). In QuTiP, the Husimi Q function can be computed givena state ket or density matrix using the function qutip.wigner.qfunc, asdemonstrated below.
In [51]: Q_coherent = qfunc(rho_coherent, xvec, xvec)In [52]: Q_thermal = qfunc(rho_thermal, xvec, xvec)In [53]: Q_fock = qfunc(rho_fock, xvec, xvec)In [54]: fig, axes = plt.subplots(1, 3, figsize=(12,3))In [55]: cont0 = axes[0].contourf(xvec, xvec, Q_coherent, 100)In [56]: lbl0 = axes[0].set_title("Coherent state")In [57]: cont1 = axes[1].contourf(xvec, xvec, Q_thermal, 100)In [58]: lbl1 = axes[1].set_title("Thermal state")In [59]: cont0 = axes[2].contourf(xvec, xvec, Q_fock, 100)In [60]: lbl2 = axes[2].set_title("Fock state")In [61]: plt.show()
Visualizing operators¶
Sometimes, it may also be useful to directly visualizing the underlying matrixrepresentation of an operator. The density matrix, for example, is an operatorwhose elements can give insights about the state it represents, but one mightalso be interesting in plotting the matrix of an Hamiltonian to inspect thestructure and relative importance of various elements.
QuTiP offers a few functions for quickly visualizing matrix data in theform of histograms, qutip.visualization.matrix_histogram andqutip.visualization.matrix_histogram_complex, and as Hinton diagram of weightedsquares, qutip.visualization.hinton. These functions takes aqutip.Qobj.Qobj as first argument, and optional arguments to, forexample, set the axis labels and figure title (see the function’s documentationfor details).
For example, to illustrate the use of qutip.visualization.matrix_histogram,let’s visualize of the Jaynes-Cummings Hamiltonian:
In [62]: N = 5In [63]: a = tensor(destroy(N), qeye(2))In [64]: b = tensor(qeye(N), destroy(2))In [65]: sx = tensor(qeye(N), sigmax())In [66]: H = a.dag() * a + sx - 0.5 * (a * b.dag() + a.dag() * b)In [67]: # visualize HIn [68]: lbls_list = [[str(d) for d in range(N)], ["u", "d"]]In [69]: xlabels = []In [70]: for inds in tomography._index_permutations([len(lbls) for lbls in lbls_list]): ....: xlabels.append("".join([lbls_list[k][inds[k]] ....: for k in range(len(lbls_list))])) ....: In [71]: fig, ax = matrix_histogram(H, xlabels, xlabels, limits=[-4,4])In [72]: ax.view_init(azim=-55, elev=45)In [73]: plt.show()
Similarly, we can use the function qutip.visualization.hinton, which isused below to visualize the corresponding steadystate density matrix:
In [74]: rho_ss = steadystate(H, [np.sqrt(0.1) * a, np.sqrt(0.4) * b.dag()])In [75]: fig, ax = hinton(rho_ss) # xlabels=xlabels, ylabels=xlabels)In [76]: plt.show()
Quantum process tomography¶
Quantum process tomography (QPT) is a useful technique for characterizing experimental implementations of quantum gates involving a small number of qubits. It can also be a useful theoretical tool that can give insight in how a process transforms states, and it can be used for example 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.
The idea is to construct a transformation matrix for a quantum process (for example a quantum gate) that describes how the density matrix of a system is transformed by the process. We can then decompose the transformation in some operator basis that represent well-defined and easily interpreted transformations of the input states.
To see how this works (see e.g. [Moh08] for more details), consider a process that is described by quantum map \(\epsilon(\rho_{\rm in}) = \rho_{\rm out}\), which can be written
(1)\[\epsilon(\rho_{\rm in}) = \rho_{\rm out} = \sum_{i}^{N^2} A_i \rho_{\rm in} A_i^\dagger,\]
where \(N\) is the number of states of the system (that is, \(\rho\) is represented by an \([N\times N]\) matrix). Given an orthogonal operator basis of our choice \(\{B_i\}_i^{N^2}\), which satisfies \({\rm Tr}[B_i^\dagger B_j] = N\delta_{ij}\), we can write the map as
(2)\[\epsilon(\rho_{\rm in}) = \rho_{\rm out} = \sum_{mn} \chi_{mn} B_m \rho_{\rm in} B_n^\dagger.\]
where \(\chi_{mn} = \sum_{ij} b_{im}b_{jn}^*\) and \(A_i = \sum_{m} b_{im}B_{m}\). Here, matrix \(\chi\) is the transformation matrix we are after, since it describes how much \(B_m \rho_{\rm in} B_n^\dagger\) contributes to \(\rho_{\rm out}\).
In a numerical simulation of a quantum process we usually do not have access to the quantum map in the form Eq. (1). Instead, what we usually can do is to calculate the propagator \(U\) for the density matrix in superoperator form, using for example the QuTiP function qutip.propagator.propagator. We can then write
\[\epsilon(\tilde{\rho}_{\rm in}) = U \tilde{\rho}_{\rm in} = \tilde{\rho}_{\rm out}\]
where \(\tilde{\rho}\) is the vector representation of the density matrix \(\rho\). If we write Eq. (2) in superoperator form as well we obtain
\[\tilde{\rho}_{\rm out} = \sum_{mn} \chi_{mn} \tilde{B}_m \tilde{B}_n^\dagger \tilde{\rho}_{\rm in} = U \tilde{\rho}_{\rm in}.\]
so we can identify
\[U = \sum_{mn} \chi_{mn} \tilde{B}_m \tilde{B}_n^\dagger.\]
Now this is a linear equation systems for the \(N^2 \times N^2\) elements in \(\chi\). We can solve it by writing \(\chi\) and the superoperator propagator as \([N^4]\) vectors, and likewise write the superoperator product \(\tilde{B}_m\tilde{B}_n^\dagger\) as a \([N^4\times N^4]\) matrix \(M\):
\[U_I = \sum_{J}^{N^4} M_{IJ} \chi_{J}\]
with the solution
\[\chi = M^{-1}U.\]
Note that to obtain \(\chi\) with this method we have to construct a matrix \(M\) with a size that is the square of the size of the superoperator for the system. Obviously, this scales very badly with increasing system size, but this method can still be a very useful for small systems (such as system comprised of a small number of coupled qubits).
Implementation in QuTiP¶
In QuTiP, the procedure described above is implemented in the function qutip.tomography.qpt, which returns the \(\chi\) matrix given a density matrix propagator. To illustrate how to use this function, let’s consider the \(i\)-SWAP gate for two qubits. In QuTiP the function qutip.gates.iswap generates the unitary transformation for the state kets:
In [77]: U_psi = iswap()
To be able to use this unitary transformation matrix as input to the function qutip.tomography.qpt, we first need to convert it to a transformation matrix for the corresponding density matrix:
In [78]: U_rho = spre(U_psi) * spost(U_psi.dag())
Next, we construct a list of operators that define the basis \(\{B_i\}\) in the form of a list of operators for each composite system. At the same time, we also construct a list of corresponding labels that will be used when plotting the \(\chi\) matrix.
In [79]: op_basis = [[qeye(2), sigmax(), sigmay(), sigmaz()]] * 2In [80]: op_label = [["i", "x", "y", "z"]] * 2
We are now ready to compute \(\chi\) using qutip.tomography.qpt, and to plot it using qutip.tomography.qpt_plot_combined.
In [81]: chi = qpt(U_rho, op_basis)In [82]: fig = qpt_plot_combined(chi, op_label, r'$i$SWAP')In [83]: plt.show()
For a slightly more advanced example, where the density matrix propagator is calculated from the dynamics of a system defined by its Hamiltonian and collapse operators using the function qutip.propagator.propagator, see notebook “Time-dependent master equation: Landau-Zener transitions” on the tutorials section on the QuTiP web site.
