Wigner State and Process Tomography on Near-Term Quantum Devices

TL;DR

This paper introduces a Wigner-based tomography method for near-term quantum devices, providing visualization tools, experimental procedures, and software to facilitate state and process characterization on quantum computers.

Contribution

It adapts a previously known ensemble-based tomography approach to pure-state quantum devices, with software tools and experimental validation on IBM quantum hardware.

Findings

Successful implementation of Wigner tomography on IBM quantum devices

Development of software package DROPStomo for practical use

Effective estimation of density and process matrices from Wigner functions

Abstract

We present an experimental scanning-based tomography approach for near-term quantum devices. The underlying method has previously been introduced in an ensemble-based NMR setting. Here we provide a tutorial-style explanation along with suitable software tools to guide experimentalists in its adaptation to near-term pure-state quantum devices. The approach is based on a Wigner-type representation of quantum states and operators. These representations provide a rich visualization of quantum operators using shapes assembled from a linear combination of spherical harmonics. These shapes (called droplets in the following) can be experimentally tomographed by measuring the expectation values of rotated axial tensor operators. We present an experimental framework for implementing the scanning-based tomography technique for circuit-based quantum computers and showcase results from IBM quantum experience. We also present a method for estimating the density and process matrices from experimentally tomographed Wigner functions (droplets). This tomography approach can be directly implemented using the Python-based software package \texttt{DROPStomo}.

Figures

Tables

Table 1. Table 1: Axial tensor operators T j ​ 0 ( ℓ ) superscript subscript 𝑇 𝑗 0 ℓ T_{j0}^{(\ell)} for one ( N = 1 𝑁 1 N=1 ) and two ( N = 2 𝑁 2 N=2 ) qubit systems.

$N$	$\ell$	$j$	$T_{j0}^{\left(\ell \right)}$
1	$\varnothing$	0	$T_{00}^{\left(\varnothing \right)}=\frac{1}{\sqrt{2}}\left(\mathrm{𝟙}\right)$
	1	1	$T_{10}^{\left(1\right)}=\frac{1}{\sqrt{2}}\left({\sigma }_z\right)$
2	$\varnothing$	0	$T_{00}^{\left(\varnothing \right)}=\frac{1}{2}\mathrm{𝟙}$
	1	1	$T_{10}^{\left(1\right)}=\frac{1}{2}\left({\sigma }_{1z}\right)$
	2	1	$T_{10}^{\left(2\right)}=\frac{1}{2}\left({\sigma }_{2z}\right)$
	12	0	$T_{00}^{\left(12\right)}=\frac{1}{2\sqrt{3}}\left({\sigma }_{1x}{\sigma }_{2x}+{\sigma }_{1y}{\sigma }_{2y}+{\sigma }_{1z}{\sigma }_{2z}\right)$
	12	1	$T_{10}^{\left(12\right)}=\frac{1}{2\sqrt{2}}\left({\sigma }_{1x}{\sigma }_{2y}-{\sigma }_{1y}{\sigma }_{2x}\right)$
	12	2	$T_{20}^{\left(12\right)}=\frac{-1}{2\sqrt{6}}\left({\sigma }_{1x}{\sigma }_{2x}+{\sigma }_{1y}{\sigma }_{2y}-2{\sigma }_{1z}{\sigma }_{2z}\right)$

Table 2. Table 2: Experimental state tomography fidelities ( ℱ s subscript ℱ 𝑠 \mathcal{F}_{s} ) corresponding to the desired single-qubit quantum states | ψ t ⟩ ket subscript 𝜓 𝑡 |\psi_{t}\rangle , see Fig.  2 .

$|{\psi }_t⟩$	${ℱ}_s$
$\frac{1}{\sqrt{2}}\left(|0⟩+|1⟩\right)$	0.9991$\pm 1\times {10}^{-3}$
$|0⟩$	0.9991$\pm 1\times {10}^{-3}$
$\frac{1}{\sqrt{2}}\left(|0⟩+i|1⟩\right)$	0.9992$\pm 1\times {10}^{-3}$
$0.885|0⟩+0.466|1⟩$	0.9990$\pm 1\times {10}^{-3}$

Table 3. Table 3: Experimental state fidelites ( ℱ s subscript ℱ 𝑠 \mathcal{F}_{s} ) corresponding to desired two-qubit quantum state | ψ t ⟩ ket subscript 𝜓 𝑡 |\psi_{t}\rangle , see Fig.  3 and 8 .

$|{\psi }_t⟩$	${ℱ}_s$
$\frac{1}{\sqrt{2}}\left(|00⟩+|11⟩\right)$	0.9989$\pm 1\times {10}^{-3}$
$\frac{1}{\sqrt{2}}\left(|00⟩+|01⟩\right)$	0.9982$\pm 1\times {10}^{-3}$

Table 4. Table 4: Experimental process tomography fidelity ( ℱ U subscript ℱ 𝑈 \mathcal{F}_{U} ) corresponding to target quantum gates U t subscript 𝑈 𝑡 U_{t} . The corresponding droplets are shown in Fig.  4 .

$U_t$	${ℱ}_U$
Hadamard (H)	$0.9506\pm 1\times {10}^{-3}$
NOT (X)	$0.9679\pm 1\times {10}^{-3}$
${\left[\frac{3\pi }{2}\right]}_y$	$0.9407\pm 1\times {10}^{-3}$

Equations

$$T_{j,\alpha \beta }^{(\ell )[N]}=R_{\alpha \beta }^{[N]}(T_{j0}^{(\ell )}{)}^{[N]}(R_{\alpha \beta }^{[N]}{)}^{†}.$$

$$R_{\alpha \beta }^{[N]}=\text{exp}(-i\alpha F_z^{[N]})\text{exp}(-i\beta F_y^{[N]})$$

$$f_j^{(\ell )}(\beta ,\alpha )=s_j⟨T_{j,\alpha \beta }^{(\ell )[N]}|A^{[N]}⟩,$$

$$f_j^{(\ell )}(\beta ,\alpha )=s_j⟨T_{j,\alpha \beta }^{(\ell )[N]}{⟩}_{A^{[N]}},$$

$$\langle{T_{j,\alpha\beta}^{(\ell)[N]}}\rangle_{A^{[N]}}=\text{tr}\big{\{}(T_{j,\alpha\beta}^{(\ell)[N]})A^{[N]}\big{\}}.$$

$$f_j^{(\ell )}(\beta ,\alpha )=s_j⟨T_{j,\alpha \beta }^{(\ell )[N]}{⟩}_{{\rho }^{[N]}}.$$

$$f_j^{(\ell )}(\beta ,\alpha )=s_j⟨T_{j0}^{(\ell )[N]}{⟩}_{{\tilde{\rho }}^{[N]}},$$

$${\tilde{\rho }}^{[N]}=(R_{\alpha \beta }^{[N]}{)}^{-1}{\rho }^{[N]}{R}_{\alpha \beta }^{[N]}.$$

$$⟨T_{10}^{(12)}{⟩}_{\tilde{\rho }(\beta ,\alpha )}$$

$$=\frac{1}{2\sqrt{2}}⟨{\sigma }_{1x}{\sigma }_{2y}{⟩}_{\tilde{\rho }(\beta ,\alpha )}-\frac{1}{2\sqrt{2}}⟨{\sigma }_{1y}{\sigma }_{2x}{⟩}_{\tilde{\rho }(\beta ,\alpha )}.$$

$$⟨{\sigma }_{1z}{\sigma }_{2z}{⟩}_{{\tilde{\rho }}_1(\beta ,\alpha )}=⟨{\sigma }_{1x}{\sigma }_{2y}{⟩}_{\tilde{\rho }(\beta ,\alpha )},$$

$${\tilde{\rho }}_1(\beta ,\alpha )=u_1\tilde{\rho }(\beta ,\alpha )u_1^{†}.$$

$$⟨{\sigma }_{1z}{\sigma }_{2z}{⟩}_{{\tilde{\rho }}_2(\beta ,\alpha )}=⟨{\sigma }_{1y}{\sigma }_{2x}{⟩}_{\tilde{\rho }(\beta ,\alpha )},$$

$${\tilde{\rho }}_2(\beta ,\alpha )=u_2\tilde{\rho }(\beta ,\alpha )u_2^{†}.$$

$$f_1^{(12)}(\beta ,\alpha )$$

$$=\frac{s_1}{2\sqrt{2}}(⟨{\sigma }_{1z}{\sigma }_{2z}{⟩}_{{\tilde{\rho }}_1(\beta ,\alpha )}+⟨{\sigma }_{1z}{\sigma }_{2z}{⟩}_{{\tilde{\rho }}_2(\beta ,\alpha )}).$$

$${\rho }^{[N]}=\sum _{a=0}^3\sum _{b=0}^3\cdots \sum _{g=0}^3{r}_{ab\dots g}({\sigma }_a\otimes {\sigma }_b\otimes \cdots \otimes {\sigma }_g)$$

$$r=⟨f_A|f_B⟩=\sum _i{\text{w}}_i{f}_A^{\ast }({\theta }_i,{\varphi }_i)f_B({\theta }_i,{\varphi }_i),$$

$${\mathcal{F}}_s=\frac{\text{tr}({\rho }^{[N]}{\rho }_t^{[N]})}{\sqrt{\text{tr}(({\rho }^{[N]}{)}^2)\text{tr}(({\rho }_t^{[N]}{)}^2)}}.$$

$$f_0^{(\varnothing )}(\beta ,\alpha )$$

$$f_1^{(1)}(\beta ,\alpha )$$

$$f_0^{(\varnothing )}(\beta ,\alpha )$$

$$f_1^{(1)}(\beta ,\alpha )$$

$$⟨\text{mathds}1⟩$$

$$⟨{\sigma }_z⟩$$

$${\rho }^{[1]}=\sum _{k=0}^3{r}_k{\sigma }_k.$$

$$r_0=⟨f_{{\sigma }_0}|f_0^{(\varnothing )}⟩.$$

$$r_k=⟨f_{{\sigma }_k}|f_1^{(1)}⟩.$$

$$\begin{array}{rl}{f}_0^{(\varnothing )}(\beta ,\alpha )& =\sqrt{\frac{1}{4\pi }}\text{expectationvalue}\ast {(T_{00}^{(\varnothing )}{)}^{[2]}}_{{\tilde{\rho }}^{[2]}}\\ & =\frac{1}{4\sqrt{\pi }}\text{expectationvalue}\ast {\text{mathds}1}_{{\tilde{\rho }}^{[2]}}\end{array}$$

$$\begin{array}{rl}{f}_1^{(1)}(\beta ,\alpha )& =\sqrt{\frac{3}{4\pi }}\text{expectationvalue}\ast {(T_{10}^{(1)}{)}^{[2]}}_{{\tilde{\rho }}^{[2]}}\\ & =\frac{1}{4}\sqrt{\frac{3}{\pi }}\text{expectationvalue}\ast {{\sigma }_{1z}}_{{\tilde{\rho }}^{[2]}},\end{array}$$

$$\begin{array}{rl}{f}_1^{(2)}(\beta ,\alpha )& =\sqrt{\frac{3}{4\pi }}\text{expectationvalue}\ast {(T_{10}^{(2)}{)}^{[2]}}_{{\tilde{\rho }}^{[2]}}\\ & =\frac{1}{4}\sqrt{\frac{3}{\pi }}\text{expectationvalue}\ast {{\sigma }_{2z}}_{{\tilde{\rho }}^{[2]}}.\end{array}$$

$$\begin{array}{rl}{f}_0^{(12)}(\beta ,\alpha )& =\sqrt{\frac{1}{4\pi }}\text{expectationvalue}\ast {(T_{00}^{(12)}{)}^{[2]}}_{{\tilde{\rho }}^{[2]}}\\ & =\frac{1}{4\sqrt{3\pi }}\text{expectationvalue}\ast {({\sigma }_{1x}{\sigma }_{2x}+{\sigma }_{1y}{\sigma }_{2y}+{\sigma }_{1z}{\sigma }_{2z})}_{{\tilde{\rho }}^{[2]}}\end{array}$$

$$\begin{array}{rl}{f}_1^{(12)}(\beta ,\alpha )& =\sqrt{\frac{3}{4\pi }}\text{expectationvalue}\ast {(T_{10}^{(12)}{)}^{[2]}}_{{\tilde{\rho }}^{[2]}}\\ & =\frac{1}{4}\sqrt{\frac{3}{2\pi }}\text{expectationvalue}\ast {({\sigma }_{1x}{\sigma }_{2y}-{\sigma }_{1y}{\sigma }_{2x})}_{{\tilde{\rho }}^{[2]}}\end{array}$$

Full text

1 Introduction

Quantum tomography is an essential tool in quantum information processing to characterize quantum systems. Here, we use a phase-space tomography approach to recover finite-dimensional Wigner representations, with a particular focus on the DROPS (Discrete Representation of OPeratorS) representation [Garon et al.(2015)Garon, Zeier, and Glaser]. The DROPS representation follows the general strategy of Stratonovich [Stratonovich(1957)], which specifies criteria for the definition of continuous Wigner functions for finite-dimensional quantum systems. This representation is based on the mapping of an arbitrary operator to a set of spherical functions denoted as droplets. It provides an intuitive visualization approach to better understand quantum systems and quantum dynamics. An example of visualization of a two-qubit state is shown in Fig. 1 and a detailed summary of the visualization approach is provided in supplementary Sec. S-I. This interactive DROPS visualization is implemented in the free SpinDrops [Tesch et al.()Tesch, Glaser, and Glaser] software.

The characteristic shapes of the droplets arising from the DROPS representation is an outcome of an abstract mapping. In our recent studies [Leiner et al.(2017)Leiner, Zeier, and Glaser, Leiner and Glaser(2018)], these shapes were related to experimentally measurable quantities, and a scanning-based tomography approach was developed to measure the droplets corresponding to quantum states and unitary processes. The procedures were experimentally implemented on an NMR quantum information processor [Cory et al.(1997)Cory, Fahmy, and Havel], an ensemble quantum computer, where expectation values of observables can be directly measured [Jones(2011)]. This paper adapts the formalism of scanning-based tomography and presents an approach to experimentally implement state and process tomography on a pure-state quantum computer. In contrast to an NMR quantum information processor, on a pure-state quantum computer, expectation values are measured by many repetitions of projective measurements on individual quantum systems [Nielsen and Chuang(2010)]. This study shows the flexibility of the scanning-based Wigner tomography approach and provides a tutorial-style approach for its implementation on current state-of-the-art pure-state quantum computing devices.

This particular tomography technique directly provides visual three-dimensional droplets based on experimental data, which helps to identify different experimental errors, such as gate imperfections, etc. We discuss this in more detail in Sec. 7, showcasing the physical intuitiveness of the DROPS representation in comparison with skyscraper visualizations [Nielsen and Chuang(2010)]. We also illustrate how to estimate the matrix representation of density and process operators based on the experimentally measured droplet functions. The theory and experimental approaches presented here can be applied on any near-term quantum device. Here, we focus on the superconducting qubit-based IBM quantum platform for performing experiments. We also provide the Python-based software package DROPStomo [Devra et al.()Devra, Glaser, Huber, and Glaser] related to this work. DROPStomo allows the use of Wigner state and process tomography on a simulator as well as on a quantum computer. This software package is discussed in Sec. 9.

2 Illustration of main results

In this section, we highlight the main results of our study. The scanning-based tomography approach estimates the expectation values of observable operators at different well-defined points on the Bloch sphere and combines them to form three-dimensional droplet functions. For example, in our experimental demonstration of state and process tomography, we use a simple equiangular sampling scheme with a combination of eight polar $\beta\in\{0,\frac{\pi}{7},\cdots\pi\}$ and fifteen azimuthal angles $\alpha\in\{0,\frac{2\pi}{14},\cdots 2\pi\}$. This corresponds to a number of sample points $N_{p}=8\cdot 15=120$. More sophisticated sampling schemes are discussed in Sec. 8. Experimental results for state tomography are shown in Fig. 2 for a single-qubit system and in Fig. 3 for a two-qubit Bell state.

For process tomography, the mapping of a unitary process matrix (U) onto a density matrix is a key step. This can be achieved by using an ancilla qubit and implementing a controlled process operation (cU) using the ancilla as a control qubit. Examples of process tomography results are showcased in Fig. 4. All the experiments in this study were performed on ibm$\_$lagos device with $N_{s}=8192$ shots per sample point, i.e., for every combination of angles $\beta$, and $\alpha$. Overall, a total number of shots $N_{tot}=N_{s}\cdot N_{p}$ were acquired. The simulated droplets shown in the results are plotted with high resolution, whereas the experimental droplets are interpolated between the experimentally determined sampling points using the Matlab surf function [Mat(2021)]. However, in the supplementary Sec. S-VIII, we also show plots of the experimental and simulated droplets with the same resolution, both of which show comparable plotting artifacts of the Matlab display function.

3 Theory of Wigner quantum state tomography

We are interested in experimentally scanning the Wigner representations of a density operator. This is a special case of quantum state tomography (QST), a vital tool in quantum computing and quantum information processing. Since the beginning of the field, there has been a lot of work performed in this direction [Vogel and Risken(1989), Leonhardt(1995), White et al.(1999)White, James, Eberhard, and Kwiat]. Recent studies use neural networks [Torlai et al.(2018)Torlai, Mazzola, Carrasquilla, Troyer, Melko, and Carleo, Neugebauer et al.(2020)Neugebauer, Fischer, Jäger, Czischek, Jochim, Weidemüller, and Gärttner] and compressed sensing [Gross et al.(2010)Gross, Liu, Flammia, Becker, and Eisert] to access the information about an unknown experimental state. Here we use a phase space [Koczor et al.(2020)Koczor, Zeier, and Glaser] tomography approach, which is helpful in experimentally visualizing quantum operators in finite-dimensional quantum systems. These, in general, can be any quantum operators such as density operators, quantum processes (propagators), etc. This section describes the scanning-based tomography approach for these operators.

3.1 Summary of scanning tomography approach

A general procedure for performing tomography in the context of Wigner representations using a scanning approach is described in the study [Leiner et al.(2017)Leiner, Zeier, and Glaser] (see results 1 and 2). Here we summarize this approach. In the following, we focus without loss of generality on a system consisting of $N$ qubits. Consider a multi-qubit quantum operator $A$, which is also represented by a set of rank $j$ and label $\ell$ spherical droplet functions $f^{(\ell)}=\sum_{j\in J(\ell)}f^{(\ell)}_{j}(\beta,\alpha)$ as described in Sec. S-I. To distinguish the size of different matrices, in the following, we refer to operators as $A^{[N]}$, where $N$ is the number of qubits and therefore, the operator $A^{[N]}$ is represented by a matrix of size $2^{N}\times 2^{N}$. The main aim is to experimentally measure spherical droplet functions $f^{(\ell)}_{j}$ representing a quantum operator $A$ and this can be done by experimentally estimating the scalar products of rotated axial tensor operators $T_{j,\alpha\beta}^{(\ell)[N]}$ with operator $A^{[N]}$, where

[TABLE]

The term $(T_{j,\alpha\beta}^{(\ell)})^{[N]}$ is the rotated version of axial tensor operators $(T_{j0}^{(\ell)})^{[N]}$ of rank $j$ and order $m=0$. The rotation operator is given by

[TABLE]

where $F_{z}=\frac{1}{2}\sum_{k=1}^{N}\sigma_{kz}^{[N]}$, and $F_{y}=\frac{1}{2}\sum_{k=1}^{N}\sigma_{ky}^{[N]}$. $R_{\alpha\beta}^{[N]}$ corresponds to a rotation around the $y$ axis by polar angle $\beta\in[0,\pi]$ followed by rotation around the $z$ axis by azimuthal angle $\alpha\in[0,2\pi)$. Here we use the shorthand notation $\sigma_{ka}=\mathds{1}\otimes\dots\otimes\mathds{1}\otimes\sigma_{a}\otimes\mathds{1}\otimes\dots\otimes\mathds{1}$, where $\sigma_{a}$ is located on the $k^{th}$ position and $a\in\{x,y,z\}$. For given angles $\beta$ and $\alpha$, droplet function $f_{j}^{(\ell)}$ can be calculated by

[TABLE]

which can be equivalently written using a shorthand notation as

[TABLE]

where $s_{j}=\sqrt{(2j+1)/(4\pi)}$ and the scalar product $\langle{T_{j,\alpha\beta}^{(\ell)[N]}}|{A^{[N]}}\rangle$ is expressed as the expectation value of $T_{j,\alpha\beta}^{(\ell)[N]}$ for the density operator ${A^{[N]}}$:

[TABLE]

Note that axial tensor operators are Hermitian, i.e., $({T_{j,\alpha\beta}^{(\ell)[N]}})^{\dagger}=({T_{j,\alpha\beta}^{(\ell)[N]}})$ [Leiner et al.(2017)Leiner, Zeier, and Glaser].

For the experimental tomography of droplet functions $f_{j}^{(\ell)}(\beta,\alpha)$ using the scanning approach, a multitude of choices for the set of sampling angles $\beta$ and $\alpha$ can be used. Fig. 5 illustrates two such sampling techniques: equiangular and REPULSION [Bak and Nielsen(1997)]. For simplicity, in the demonstration experiments shown here, we use an equiangular grid, but a numerical study using more sophisticated sampling techniques is presented and discussed in Sec. 8.

3.2 Wigner quantum state tomography

In this section, we first present an algorithm for Wigner quantum state tomography and then elaborate each step individually. A general step-wise procedure (see Fig. 6) to experimentally measure a spherical droplet function $f^{(\ell)}_{j}$ representing a density matrix $\rho^{[N]}=|\psi^{[N]}\rangle\langle\psi^{[N]}|$ is the following:

1. (i)

Preparation ($\mathcal{P}$): Prepare the desired quantum state $\rho^{[N]}$ from a defined initial state $\rho_{i}^{[N]}$. 2. (ii)

Rotation ($\mathcal{R}$): Rotate the density operator $\rho^{[N]}$ inversely for scanning. 3. (iii)

Detection-associated rotations ($\mathcal{D}$): Apply local unitary operations to measure expectation values of Pauli operator components of axial tensor operators $T_{j0}^{(\ell)[N]}$ (see Table 1) that are not directly measurable.

These steps are repeated for the set of angles $\beta\in[0,\pi]$, $\alpha\in[0,2\pi]$ and for different local unitary operators $u_{n}$ (vide infra), rank $j$, and label $\ell$ to experimentally scan the droplet functions $f_{j}^{(\ell)}$. In the rest of this section, we elaborate on each step of the presented algorithm.

Step 1: The first step of the algorithm, ‘Preparation’, can be achieved by applying unitary operations depending on the initial and desired state.

Step 2: Since our operator of interest is a density matrix $\rho^{[N]}$, Eq. 4 takes the following form:

[TABLE]

Instead of rotating the axial tensor operators $T_{j0}^{(\ell)[N]}$ as shown in Eq. 1, it is equivalent (and experimentally more convenient) to rotate the density matrix $\rho^{[N]}$ inversely, such that:

[TABLE]

where

[TABLE]

The axial tensor operators $(T_{j0}^{(\ell)})$ are explicitly given in [Garon et al.(2015)Garon, Zeier, and Glaser, Leiner et al.(2020)Leiner, Zeier, and Glaser] for systems consisting of up to six qubits. In Table 1, we summarize the axial tensor operators for one and two-qubit systems.

Step 3: Depending on the number of qubits $N$ and rank $j$, the axial tensors $T_{j0}^{(\ell)}$ consist of different Pauli operators (see Table 1), but it might not be possible to measure these components directly depending on the specific quantum computing hardware. In a typical pure-state quantum computing device, the measurement is done along the $z$ axis, which implies that the directly measurable operators are: $\mathds{1}$, $\sigma_{1z}$, $\sigma_{2z}$, and $\sigma_{1z}\sigma_{2z}$. In this case, measuring expectation values of non-directly measurable operators can be achieved with the help of local unitary operations $u_{n}$. For example, consider the expectation value $\langle T_{10}^{(12)}\rangle_{\tilde{\rho}(\beta,\alpha)}$ given by:

[TABLE]

In the first term ($n=1$) of Eq. 9, the expectation value $\langle{\sigma_{1x}\sigma_{2y}}\rangle_{\tilde{\rho}(\beta,\alpha)}$ needs to be determined. This can be achieved by measuring instead the expectation value

[TABLE]

where

[TABLE]

The density operator $\tilde{\rho}_{1}(\beta,\alpha)$ is obtained from $\tilde{\rho}(\beta,\alpha)$ by applying a $-\pi/2$ rotation around the $y$ axis (for bringing the $x$ axis to the $z$ axis) to the first qubit and a $\pi/2$ rotation around the $x$ axis (for bringing the $y$ axis to the $z$ axis) to the second qubit. This corresponds to the local unitary transformation u_{1}=\big{(}(R_{0,\frac{\pi}{2}})^{-1}\otimes\mathds{1}\big{)}\cdot(\mathds{1}\otimes(R_{\frac{\pi}{2},\frac{\pi}{2}})^{-1}).

Similarly, in the second term ($n=2$) of Eq. 9, the expectation value $\langle{\sigma_{1y}\sigma_{2x}}\rangle_{\tilde{\rho}(\beta,\alpha)}$ needs to be determined. This can be achieved by measuring instead the expectation value

[TABLE]

where

[TABLE]

The density operator $\tilde{\rho}_{2}(\beta,\alpha)$ is obtained from $\tilde{\rho}(\beta,\alpha)$ by applying a $\pi/2$ rotation around the $x$ axis (for bringing the $y$ axis to the $z$ axis) to the first qubit and a $-\pi/2$ rotation around the $y$ axis (for bringing the $x$ axis to the $z$ axis) to the second qubit. This corresponds to the local unitary transformation u_{2}=\big{(}(R_{{\frac{\pi}{2}},\frac{\pi}{2}})^{-1}\otimes\mathds{1}\big{)}\cdot(\mathds{1}\otimes(R_{0,\frac{\pi}{2}})^{-1}).

Overall, for $\ell=\{12\}$ and $j=1$, the droplet function from Eq. 7 can be expressed as

[TABLE]

Hence, only projective measurements along the $z$ axis are required, c.f. S-III.

3.3 Estimation of a density matrix from the droplet functions

In this section, we show how the matrix form of a density operator can be estimated based on its experimentally measured DROPS representation. A general $N$-qubit density matrix [Nielsen and Chuang(2010)] can be expressed in terms of Pauli operators as:

[TABLE]

where $\sigma_{0}$ is $\mathds{1}$ (the 2$\times$2 identity matrix), while $\sigma_{1},\sigma_{2}$, and $\sigma_{3}$ are standard Pauli matrices $\sigma_{x},\sigma_{y}$, and $\sigma_{z}$ respectively. The terms $r_{ab\dots g}$ are real coefficients. Given the DROPS representation of a density operator, these coefficients can be computed by calculating the scalar product between basis droplets (ideally simulated without noise) and experimental droplets [Leiner et al.(2017)Leiner, Zeier, and Glaser]. The basis droplets can be generated using the definitions provided in supplementary Sec. S-VI for one and two qubits.

The scalar product between two tensor operators can be approximated by the discretized scalar product between their droplet functions, see supplementary Sec. S-II. In the general case of a droplet $f_{A}$ with complex values $f_{A}(\theta_{i},\phi_{i})$ and another droplet $f_{B}$ with complex values $f_{B}(\theta_{i},\phi_{i})$ at the grid points, the scalar product or overlap between the two droplets is given by

[TABLE]

where $f_{A}^{*}(\theta_{i},\phi_{i})$ is the complex conjugate of $f_{A}(\theta_{i},\phi_{i})$. The sampling weights $\text{w}_{i}$ corresponding to an equiangular grid for calculating the scalar product are provided in supplementary Sec. S-II. The sampling weights for other sampling techniques, such as REPULSION and Lebedev, are available in [Hogben et al.(2011)Hogben, Krzystyniak, Charnock, Hore, and Kuprov], and [University()], respectively.

Based on the estimated experimental density operator $\rho^{[N]}$, the state fidelity ($\mathcal{F}_{s}$) [Liang et al.(2019)Liang, Yeh, Mendonça, Teh, Reid, and Drummond] with which a desired state $\rho_{t}^{[N]}$ has been reached can be calculated using the normalized scalar product:

[TABLE]

In the next section, we focus on the experimental implementation of the presented algorithm for one and two-qubit systems and showcase the experimental results performed on the ibm_lagos device.

4 Experimental implementation of Wigner state tomography

In this section we discuss how the Wigner quantum state tomography algorithm can be implemented on an experimental near-term quantum device. For concreteness, we will present the implementation using Qiskit [et al.(2019)], an open-source quantum development platform for simulations and experiments. Wigner state tomography can be used directly for one- and two-qubit systems using the Python-based software package DROPStomo [Devra et al.()Devra, Glaser, Huber, and Glaser].

4.1 One qubit

For a system consisting of one qubit ($N=1$), there are only two possible values for the rank $j$: $j=0$ for $\ell=\{\emptyset\}$ and $j=1$ for $\ell=\{1\}$. Hence, the single-qubit density matrix $\rho^{[1]}$ represented by spherical functions $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$ can be tomographed by measuring the expectation values from Eq. 7 as

[TABLE]

Substituting the explicit form of the tensor operators $T_{00}^{(\emptyset)}$ and $T_{10}^{(1)}$ from Table 1 gives

[TABLE]

where $\mathds{1}$ is the $2\times 2$ identity matrix. The required expectation values can be computed directly on a pure-state quantum computer based on the outcome probabilities $p_{0}$ and $p_{1}$ corresponding to state $|0\rangle$ and state $|1\rangle$ respectively:

[TABLE]

The probabilities $p_{0}$ and $p_{1}$ can be experimentally measured by repeating an experiment multiple times, where the number of repetitions is also referred to as number of shots $N_{s}$. We discuss this in more detail in supplementary Sec. S-III. In the considered case, the detection-associated rotation ($\mathcal{D}$) step is not required.

The probabilities $p_{0}$ and $p_{1}$ required for computing the expectation values for droplets $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$ can be measured in the same experiment; hence the algorithm does not require to be repeated for rank $j$ and label $\ell$. In supplementary Sec. S-V, we explicitly provide the quantum circuit (Fig. S3) for performing tomography of a qubit in the state $|\psi\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle$ along with the corresponding plot of simulated and experimental expectation values (Fig. S3) of droplet function $f_{1}^{(1)}(\beta,\alpha)$.

Based on the experimentally tomographed droplet functions $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$, the corresponding density matrix can be estimated as discussed in Sec. 3.3. For a single qubit, the density matrix can be expressed in terms of Pauli operators as

[TABLE]

To compute the coefficient $r_{0}$, a scalar product is calculated using Eq. 16 between the ideal basis droplet ($f_{\sigma_{0}}$) and the experimentally tomographed $f_{0}^{(\emptyset)}$ droplet:

[TABLE]

To compute the remaining coefficients $r_{k}$ for $k\in\{1,2,3\}$, we calculate the scalar product between all ideal basis droplets with label $\ell=1$, rank $j=1$, and the experimentally tomographed rank $j=1$ droplet $f_{1}^{(1)}$:

[TABLE]

The ideal basis droplet functions $f_{\sigma_{k}}$ are provided in supplementary Sec. S-VI.I. Hence, by calculating these overlap coefficients $r_{k}$, an estimate of the density matrix can be obtained using Eq. 21. In Table 2, we provide experimental state tomography fidelities computed using Eq. 17 for different states considered in Fig. 2. The standard deviations given in Table 2 were estimated by conducting the experiment three times with the state $|\psi_{t}\rangle=0.885|0\rangle+0.466|1\rangle$, for a reference. To avoid redundancy, in Fig. 2, we only plot the rank $j=1$ droplets, as the rank $j=0$ droplet is a sphere of radius $\sqrt{1/(8\pi)}$ in the single-qubit case.

4.2 Two qubits

For a two-qubit system ($N=2$), there are four possible labels $\ell$, and based on these labels, there are different ranks $j$ as shown in Table 1. Therefore, for a two-qubit system, the Wigner quantum state tomography requires measuring the following spherical droplets $f_{j}^{(\ell)}$:

[TABLE]

for the identity droplet ($j$=0, and $\ell$=$\emptyset$), here $\mathds{1}$ is the 4$\times$4 identity matrix. For droplets of rank $j=1$, for each qubit (with labels $\ell=1$ and $\ell=2$) we have:

[TABLE]

[TABLE]

For bilinear terms with label $\ell=12$, we have

[TABLE]

for rank $j=0$,

[TABLE]

for rank $j=1$, and

[TABLE]

for rank $j=2$. Different Pauli operator expectation values are required in Eq. 25 to Eq. 29, and some of them are not directly measurable. In this case, we use the last block of the algorithm called detection-associated rotations ($\mathcal{D}$) as explained with an example in Sec. 3.2. The rotations required for step $\mathcal{D}$ can be implemented in terms of local $\mathrm{U}_{3}$ gates as described in supplementary Sec. S-IV. Similar to the computation of expectation values of linear terms in Eq. 20, the expectation values of bilinear terms can be computed by combining the outcome probabilities:

[TABLE]

where $p_{ab}$ for $a,b\in\{0,1\}$ is the probability corresponding to state $|ab\rangle$. We refer to supplementary Sec. S-III for more information. In Fig. 7, we explicitly show the quantum circuits for state tomography of a Bell state $|\psi\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$. The different circuits in the figure are used for calculating the different expectation values, which are then combined to form a particular droplet function $f_{j}^{(\ell)}$. Hence, for a two-qubit Wigner quantum state tomography, a maximum of five quantum circuits are required, which are repeated for all combinations of values of the angles $\beta$ and $\alpha$.

Similar to the one-qubit system, we can estimate the density matrix based on the experimentally tomographed droplets. The density matrix for a two-qubit system can be expressed in terms of Pauli operators as

[TABLE]

The terms $r_{kl}$ with $k,l\in\{0,1,2,3\}$ are real coefficients and can be calculated by computing the scalar product between the droplet functions as shown in Eq. 16. To compute $r_{00}$, the scalar product is calculated between the simulated ideal basis droplet ($f_{\sigma_{0}}$) with label $\ell=\emptyset$, rank $j=0$ and the experimentally tomographed droplet function $f_{0}^{(\emptyset)}$ as,

[TABLE]

To calculate the coefficients $r_{k0}$ for $k\in\{1,2,3\}$ the scalar product is computed between all ideal basis droplets with label $\ell=1$, rank $j=1$ and the experimentally tomographed droplet $f_{1}^{(1)}$:

[TABLE]

Similarly, the coefficients $r_{0l}$ for $l\in\{1,2,3\}$ can be computed by calculating the scalar product between the ideal basis droplets of label $\ell=2$, rank $j=1$ and experimentally tomographed droplet $f_{1}^{(2)}$:

[TABLE]

The remaining bilinear coefficients $r_{kl}$ for $k,l\in\{1,2,3\}$ can be calculated by computing the scalar product between the ideal bilinear basis droplets $f_{\sigma_{1k}\sigma_{2l}}$ with the sum of the experimentally tomographed droplets $f_{0}^{(12)}$, $f_{1}^{(12)}$, and $f_{2}^{(12)}$:

[TABLE]

The ideal basis droplets are provided in the supplementary Sec. S-VI.II. Hence, by calculating these coefficients for every value of $k$ and $l$, a density matrix can be estimated using Eq. 31 and the state fidelity ($\mathcal{F}_{s}$) can be computed using Eq. 17. In Table 3, we present the experimental state fidelities of two-qubit examples. The standard deviation given in Table 3 were estimated by conducting the experiment three times with the state $|\psi_{t}\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$, for a reference. In Fig. 3 and Fig. 8, we show the experimentally tomographed and theoretical droplets for a maximally entangled Bell state and a separable quantum state, respectively. The bilinear droplets are combined to $f^{(12)}=f_{0}^{(12)}+f_{1}^{(12)}+f_{2}^{(12)}$. In Sec. S-VIII we also provide plots of droplets corresponding to individual ranks. The separable state used here as an example for tomography has also been used as an example for visualization in Fig. 1.

In the next section, we take this scanning tomography approach a step further and explain how this can be used to experimentally tomograph shapes representing unitary quantum processes.

5 Theory of Wigner quantum process tomography

In general, process tomography [Nielsen and Chuang(2010), Chuang and Nielsen(1997), Childs et al.(2001)Childs, Chuang, and Leung] is a method to characterize a quantum process based on experimentally measured data. Here we focus on the tomography of unitary processes and translate the previously introduced Wigner quantum process tomography of known propagators [Leiner and Glaser(2018)] in ensemble quantum devices to the setting of pure state near-term quantum devices. The considered unitary processes could refer to quantum gates, time evolution operators, or pulse sequences.

As shown in Sec. 3.2, if the operator of interest is a quantum state or a density operator represented by $\rho^{[N]}$, the spherical droplet function $f_{j}^{(\ell)}(\beta,\alpha)$ can be measured experimentally. In the case of process tomography, our operator of interest is an $N$ qubit quantum process $U^{[N]}$. It is possible to scan the Wigner representation of an arbitrary operator $A^{[N]}$ if it can be experimentally mapped onto the density operator. In the next section, we present the algorithm for Wigner process tomography based on a method to map a unitary process matrix onto a density matrix [Leiner and Glaser(2018), Fahmy et al.(2008)Fahmy, Marx, Bermel, and Glaser, Myers et al.(2001)Myers, Fahmy, Glaser, and Marx, Marx et al.(2010)Marx, Fahmy, Kauffman, Lomonaco, Spörl, Pomplun, Schulte-Herbrüggen, Myers, and Glaser].

5.1 Mapping of a unitary process matrix onto a density matrix

Mapping a general unitary matrix $U^{[N]}$ onto a Hermitian density matrix of the same dimension is not possible. Here, we double the dimension of the density matrix by using an additional ancilla qubit $q_{0}$ and by a controlled process $cU^{[N+1]}$ operation we inscribe the unitary $U^{[N]}$ (and its adjoint $(U^{[N]})^{\dagger}$) in an off-diagonal block of the density matrix $\rho^{[N+1]}$ as shown below. Under $cU^{[N+1]}$, the unitary $U^{[N]}$ acts only on the target qubits $\textit{q}_{1},\dots,\textit{q}_{N}$ if the control qubit $\textit{q}_{0}$ is in state $\ket{1}$. The corresponding matrix representation of the controlled process $cU^{[N+1]}$ is

[TABLE]

where the top diagonal block corresponds to a $2^{N}\times 2^{N}$ dimensional identity matrix $\mathds{1}^{[N]}$ and the lower diagonal block is the unitary $U^{[N]}$. The off-diagonal blocks are $2^{N}\times 2^{N}$ dimensional zero matrices.

As shown in [Leiner and Glaser(2018)] for ensemble quantum processors, $U^{[N]}$ can be mapped onto the density operator $\rho^{[N+1]}$ by preparing the ancilla (control) qubit $\textit{q}_{0}$ in the superposition state $\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ and the remaining system qubits $\textit{q}_{1},\dots,\textit{q}_{N}$ in the maximally mixed state. Hence, the prepared density operator is

[TABLE]

and the density operator after application of $cU^{[N+1]}$ is

[TABLE]

which can be rewritten in block matrix form as

[TABLE]

Using this approach, the unitary $U^{[N]}$ is now imprinted onto the density operator $\rho^{[N+1]}$ of the augmented system. Since the experimental implementation of a controlled process $cU^{[N+1]}$ requires the knowledge of $U^{[N]}$, this version of Wigner process tomography described here is in general only applicable for known processes [Leiner and Glaser(2018)].

5.2 Wigner quantum process tomography

Here, we first present the algorithm for process tomography and then explain each step individually for pure-state quantum processors. A droplet function $f_{j}^{(\ell)}$ representing a quantum process $U^{[N]}$ can be experimentally measured using the following steps (see Fig. 9):

1. (i)

Preparation ($\mathcal{P}$): Prepare ancilla qubit $\textit{q}_{0}$ in the superposition state $\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ and effectively create the maximally mixed state of the system qubits $\textit{q}_{1},\dots,\textit{q}_{N}$ by temporal averaging. 2. (ii)

Mapping ($\mathcal{M}$): Implement the $cU^{[N+1]}$ operation to map the process $U^{[N]}$ onto $\rho_{U}^{[N+1]}$. 3. (iii)

Rotation ($\mathcal{R}$): Rotate the system qubits $\textit{q}_{1},\dots,\textit{q}_{N}$ inversely for scanning. 4. (iv)

Detection-associated rotations ($\mathcal{D}$): Apply local unitary operations to measure required expectation values of Pauli operator components of axial tensor operators $T_{j0}^{(\ell)[N]}$ (see Table 1) that are not directly measurable.

These four steps are repeated for a set of angles $\beta\in[0,\pi]$ and $\alpha\in[0,2\pi]$ and for different $n$, rank $j$ and labels $\ell$ to calculate the droplet function $f_{j}^{(\ell)}(\beta,\alpha)$. Now we elaborate each step individually.

Step 1: The algorithm starts with the initial state $\rho_{i}=|0\dots 0\rangle\langle 0\dots 0|$. The state $\rho_{0}$ (Eq. 37) is prepared by applying a Hadamard gate to qubit $q_{0}$ to achieve an equal superposition. The maximally mixed state of qubits $q_{1},\dots,q_{N}$ is created by temporally averaging experiments for all the possible computational basis states by applying appropriate local NOT gates (see supplementary Sec. S-VII). We discuss this in detail for a single-qubit system in Sec. 6.

Step 2: The operation $cU^{[N+1]}$ used for mapping can be experimentally implemented by decomposing it into elementary gates [Barenco et al.(1995)Barenco, Bennett, Cleve, DiVincenzo, Margolus, Shor, Sleator, Smolin, and Weinfurter] or using pulse-level control methods [Khaneja et al.(2005)Khaneja, Reiss, Kehlet, Schulte-Herbrüggen, and Glaser, Devra et al.(2018)Devra, Prabhu, Singh, Arvind, and Dorai].

Step 3: Since our operator of interest is an $N$-qubit unitary process matrix $U^{[N]}$, Eq. 3 takes the form

[TABLE]

As shown in supplementary Sec. S-IX, Eq. 40 can be rewritten as

[TABLE]

Similarly as in state tomography (Sec. 3.2), instead of rotating the axial tensor operator $T_{j0}^{(\ell)[N]}$, we equivalently rotate the density matrix of the system qubits $\textit{q}_{1},\dots,\textit{q}_{N}$ inversely, such that:

[TABLE]

where

[TABLE]

and $R_{\alpha\beta}^{[N+1]}=\mathds{1}^{[1]}\otimes R_{\alpha\beta}^{[N]}$ which corresponds to the rotation of only the system qubits $\textit{q}_{1},\dots,\textit{q}_{N}$ for scanning. Using the relation $\sigma^{+}=\frac{1}{2}(\sigma_{x}+i\sigma_{y})$, Eq. 42 can be rewritten in terms of Pauli operators as:

[TABLE]

Step 4: In analogy to the case of Wigner state tomography, the expectation values of Pauli operators, which are not directly observable, can be measured with the help of local unitary operations $u_{n}$ (detection-associated rotations).

5.3 Estimation of unitary process matrices from droplet functions

Similar to the estimation of density matrices in the case of Wigner state tomography, unitary process matrices can also be estimated from the experimentally tomographed droplets. A general $N$-qubit unitary process matrix can be expressed in terms of Pauli operators as:

[TABLE]

where $\sigma_{0}$ is a 2$\times$2 identity ($\mathds{1}$) matrix, while $\sigma_{1}$, $\sigma_{2}$ and $\sigma_{3}$ are the standard Pauli matrices $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$, respectively. The complex coefficients $c_{ab\dots g}$ can be computed by calculating the scalar product between basis droplets (ideally simulated without noise) and experimental droplets, as shown in Eq. 16. The basis droplets can be generated using the definitions provided in supplementary Sec. S-VI. Based on the estimated process matrix $U^{[N]}$, the process tomography fidelity $\mathcal{F}_{U}$ [Glaser et al.(1998)Glaser, Schulte-Herbrüggen, Sieveking, Schedletzky, Nielsen, Sørensen, and Griesinger] can be calculated using the relation:

[TABLE]

where $U^{[N]}_{t}$ is a target unitary process matrix.

6 Experimental implementation of Wigner process tomography

This section describes the approach to implementing the above-mentioned Wigner process tomography on experimental quantum devices. Here, we present the simulated and experimental process tomography results performed on IBM quantum devices for a pure state of an individual quantum system. The quantum circuits provided here are general and can be directly adapted to other near-term quantum devices. The Wigner process tomography can be directly implemented using the Python-based software package DROPStomo [Devra et al.()Devra, Glaser, Huber, and Glaser] for a single-qubit system.

6.1 Single qubit system

For the Wigner process tomography of a single-qubit ($N=1$) system, the total number of qubits required is two ($q_{0}$ and $q_{1}$), where $q_{0}$ is an ancilla qubit and $q_{1}$ is the system qubit. For a single-qubit system, the possible values of rank $j$ are (c.f. Table 1): $j=0$ for label $\ell=\emptyset$, and $j=1$ for label $\ell=1$. Hence, a single-qubit unitary process is represented by the spherical functions $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$, which can be calculated based on the measured expectation values of Eq. 44 as:

[TABLE]

Substituting the explicit form of the tensor operators $T_{00}^{(\emptyset)}$ and $T_{10}^{(1)}$ from Table 1 gives

[TABLE]

We first focus on the preparation step ($\mathcal{{P}}$) of the algorithm, i.e., preparing qubits $\textit{q}_{0}$ and $\textit{q}_{1}$ in a state whose density matrix corresponds to Eq. 37 for $N=1$. The preparation of the ancilla qubit $q_{0}$ in the superposition state can be achieved straightforwardly by applying the Hadamard (H) gate on $q_{0}$. A relatively straightforward approach to preparing the system qubit in the maximally mixed state would be to prepare it in the state $|0\rangle$ and to repeat each experiment multiple times, where, in each repetition, it is randomly decided (with 50% probability) whether to flip the qubit using a NOT gate. However, the standard deviation of the average population of state $|1\rangle$ from the expected value of 0.5 decreases only with the inverse of the square root of the number of repetitions in this probabilistic approach. Hence, a large number of repetitions would be required to minimize this additional noise source due to the imperfect realization of the completely mixed state. For example, to achieve a standard deviation of less than 0.005, about 8000 repetitions would be necessary, even in the case of perfect gate operation. In contrast, the completely mixed state can be created exactly using a temporal averaging approach [Knill et al.(1998)Knill, Chuang, and Laflamme, Preskill(1998)], where the experiment is only repeated for the set of computational basis states. In the case of a single system qubit, this requires only two experiments: one experiment without a flip of initial state $|0\rangle$ and one experiment with a flip, see supplementary Sec. S-VII.

Interestingly, temporal averaging was initially introduced in quantum information processing to mimic experiments of a pure state by averaging expectation values obtained by measuring a set of experiments on (partially) mixed states of an ensemble quantum processor [Knill et al.(1998)Knill, Chuang, and Laflamme]. In contrast, here experiments on the maximally mixed state of the system qubits $\textit{q}_{1},\dots,\textit{q}_{N}$ are mimicked by averaging expectation values obtained by measuring a set of experiments with pure states.

As an example, in Fig. 10 we present the set of quantum circuits for process tomography where the unitary operator of interest is the Hadamard gate (H). Hence, a controlled Hadamard gate is applied in the mapping step ($\mathcal{M}$) of each experiment in (a)-(d). After the preparation step $\mathcal{P}$, the system qubit $q_{1}$ is in state $|0\rangle$ for circuits (a) and (c), whereas it is switched to $|1\rangle$ by applying a NOT (X) gate on $q_{1}$ in circuits (b) and (d). The (temporally) averaged expectation values of experiments (a) and (b) provide $\langle{\sigma_{0x}}\rangle$ (and also $\langle{\sigma_{0x}\sigma_{1z}}\rangle$). Similarly, the (temporally) averaged expectation values of experiments (c) and (d) provide $\langle{\sigma_{0y}}\rangle$ (and also $\langle{\sigma_{0y}\sigma_{1z}}\rangle$). Using Eq. 48, these expectation values can be combined to obtain the droplet functions $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$, see Fig 4 (first column) for the combined ($f=f_{0}^{(\emptyset)}+f_{1}^{(1)}$) droplets [Leiner and Glaser(2018)] and Fig. S6 (first row) for the individual droplets $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$. Fig. 4 also shows experimental droplets of the NOT (X) gate and the process corresponding to a rotation of $\frac{3\pi}{2}$ around the $y$ axis. The separate droplets $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$ for rank $j=0$ and $j=1$ are provided in the supplementary Fig. S6.

Based on the experimentally measured droplet functions $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$, a process matrix can be estimated (as shown in Sec. 5.3): any single-qubit unitary process can be expressed in terms of Pauli operators [Nielsen and Chuang(2010)] as

[TABLE]

with in general complex coefficients $c_{k}$ for $k\in\{0,1,2,3\}$. Using Eq. 16, the coefficients $c_{k}$ are obtained by calculating the scalar product between the basis droplets $f_{\sigma_{k}}$ (refer to Sec. S-VI.I) and the sum of the experimentally tomographed droplets $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$:

[TABLE]

Table 4 summarizes the experimental process fidelities of the gates considered above. The standard deviation given in Table 4 were estimated by conducting the experiment three times for the NOT gate for a reference.

7 Understanding errors using Wigner state and process tomography

Quantum devices are prone to different kinds of errors both in the implementation of desired states and of quantum gates. Here, we focus on the example of rotation errors, e.g., due to errors in pulse calibrations, etc. Visualizing or identifying these errors directly is useful in quantum information processing. Here, we show how the DROPS representation is helpful to achieve this. As described in the caption of Fig. 1, the radius (distance from the origin to a point on the sphere) of a droplet represents the absolute value of a droplet function $f^{(\ell)}$, and color represents the phase $\varphi=\text{arg}[f^{(\ell)}]$. In addition, the direction of a qubit droplet reflects the direction of the Bloch vector for quantum states (see: Fig. 1 and 2) and of the rotation axis for single-qubit quantum gates (see Fig. 4).

As an example, we show the experimental tomography result of the desired quantum state $|\psi\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|01\rangle)$ with some error in the state preparation. We deliberately introduce an additional rotation of $\mathrm{U}_{3}(\pi/12,0,0)$ on qubit $q_{1}$, and $\mathrm{U}_{3}(\pi/9,\pi/12,0)$ on qubit $q_{2}$ in the preparation step. Fig. 11 shows that these kinds of errors are directly visible in the DROPS representations (on the right) as misalignment of the linear droplet functions $f^{(1)}$ and $f^{(2)}$ compared to the experimental tomography results of the case without rotation errors as shown in Fig. 8 and to the ideal case shown in Fig 1. Note that such a direct physical interpretation of the error terms is not possible using the standard skyscraper visualization [Nielsen and Chuang(2010)] of the density matrix (on the left). In Fig. 11, only the skyscraper visualization of the real part of the density matrix is shown. The imaginary part is plotted in Fig. S9 along with decomposed bilinear droplet functions for the state with and without rotation errors. In Fig. S10, we also show the droplet plots from a different perspective to emphasize the misalignment errors.

8 Scanning using different sampling schemes: a numerical study

As we have seen, scanning is a key step in the Wigner tomography approach, and hence choosing a suitable sampling scheme on a sphere is important. This is a topic of interest in the general field of signal processing [Kennedy and Sadeghi(2013)], and a number of different sampling schemes have been proposed in the literature.

In the ideal case of negligible experimental noise, only a small number of sampling points would be necessary to determine the correct expansion coefficients of spherical harmonics as each droplet function is band-limited [Khalid et al.(2014)Khalid, Kennedy, and McEwen, Leiner and Glaser(2018)]. An advantage of using a large number of sampling points is to obtain a direct view of the experimentally measured droplet shapes without additional signal processing steps, such as the extraction of expansion coefficients of spherical harmonics or the estimation of the matrix representation of an operator. Note that a larger number of sampling points $N_{p}$ does not necessarily entail an increase in the total number of experiments $N_{tot}(=N_{p}\cdot N_{s})$ because the number $N_{s}$ of shots per sampling point can also be adapted to each sampling scheme. In the following, we will compare the performance of different sampling schemes for (approximately) the same total number of shots ($N_{tot}$) given by the product of the number of sampling points ($N_{p}$) and the number of shots per sampling point ($N_{s}$).

In Fig. 12 the mean fidelity ($\bar{\mathcal{F}_{s}}$) of the tomographed state is shown as a function of the total number of experiments for the Lebedev [Lebedev(1976)], REPULSION [Bak and Nielsen(1997)], and SHREWD [Edén and Levitt(1998)] sampling schemes along with the simple equiangular grid and the standard tomography method [Nielsen and Chuang(2010), James et al.(2001)James, Kwiat, Munro, and White]. For more detailed information, including standard deviations, see supplementary Sec. S-X. We only consider the noise due to a limited number of shots. In the simple case of an equiangular grid [Driscoll and Healy(1994), McEwen and Wiaux(2011)] of eight polar angles $\beta\in\{0,\frac{\pi}{7},\dots\pi\}$ and fifteen azimuthal angles $\alpha\in\{0,\frac{2\pi}{14},\dots 2\pi\}$ as shown in Fig. 5a, the total number of grid points is 120. For Lebedev, REPULSION, and SHREWD 110 grid points were used. Since the number of sampling points in both cases is similar, for simplicity, the same number of shots per sampling point was chosen. In contrast, for the standard tomography method, only three measurement settings are required for the case of a single-qubit. The forty-fold decrease in the number of sampling points was compensated by correspondingly increasing the number of shots per measurement setting by a factor of 40 to arrive at the same total number of shots $N_{tot}$ as in the previous cases.

In the supplementary Sec. S-X, we also provide results for the state $|\psi\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$. The plots indicate that the mean fidelity differs depending on sampling techniques and tomography methods and also show that an equiangular grid is not an optimal choice, as expected because the density of grid points is higher near the poles compared to the equator of the sphere. This numerical study is expected to help make an informed choice of the sampling scheme for quantum computing applications in which an estimate of the quantum state with high precision is required.

We used the Spinach [Hogben et al.(2011)Hogben, Krzystyniak, Charnock, Hore, and Kuprov] software to generate angles and weights for the REPULSION and the SHREWD sampling techniques. For standard state tomography, the maximum-likelihood estimation (MLE) method [Smolin et al.(2012)Smolin, Gambetta, and Smith, Singh et al.(2016)Singh, Arvind, and Dorai] was used on the numerical data to estimate a valid quantum state. We used the corresponding qiskit [et al.(2019)] classes to perform the standard state tomography based on the MLE method.

9 DROPStomo: A Python-based software package for Wigner state and process tomography

DROPStomo [Devra et al.()Devra, Glaser, Huber, and Glaser] is a Python-based software package for performing Wigner state tomography for one- and two-qubit systems and process tomography for a single-qubit system. With DROPStomo, users can simulate (on a simulator or on quantum hardware) and analyze the tomographed droplets interactively. The package is based on the Qiskit framework [et al.(2019)]. However, it is straightforward to adapt it to other frameworks. DROPStomo can be installed and imported using the following command:

Here, we give a code snippet for performing Wigner state tomography for the one-qubit state $|\psi\rangle=\frac{|0\rangle+|1\rangle}{\sqrt{2}}$ for eight polar angles $\beta\in[0,\pi]$ and fifteen azimuthal angles $\alpha\in[0,2\pi]$.

We provide the extended tutorial codes for two-qubit Wigner state tomography and one-qubit Wigner process tomography in our repository [Devra et al.()Devra, Glaser, Huber, and Glaser].

10 Discussion

In this work, we developed a general approach for Wigner tomography of quantum states and quantum processes for pure-state quantum devices with projective measurements by adapting the methodology described in [Leiner et al.(2017)Leiner, Zeier, and Glaser] and [Leiner and Glaser(2018)]. We demonstrated the experimental implementation of these tomography approaches on IBM quantum devices. The experimentally measured droplet shapes provide a unique, vivid visual representation of abstract quantum operators, which reflects their essential features. For example, for a single-qubit, the droplet orientation provides essentially the same information as the Bloch vector representation. However, the DROPS representation is not limited to single-qubit systems but can also be used to visualize states and processes in multi-qubit systems, where different droplets provide information about the expectation values of multi-qubit operators, such as $\sigma_{1x}\sigma_{2x}$ or $\sigma_{1y}\sigma_{2z}\sigma_{3y}$ etc. The presented approach has similar limitations with respect to the number of qubits due to the exponential growth of the Hilbert space dimension as conventional tomography methods. However, for a small number of qubits, Wigner-type DROPS tomography forms an easily implementable alternative approach with additional benefits and without any additional cost in terms of experimental time and signal-to-noise ratio. In particular, the DROPS visualization allows one to directly see the kind of errors present in a given realization of quantum states and processes.

Originally, Wigner state and process tomography was developed for ensemble quantum processors, such as NMR [Leiner et al.(2017)Leiner, Zeier, and Glaser, Leiner and Glaser(2018)]. The main purpose of this paper was to show that it is indeed possible to also apply it to standard quantum devices based on pure states. To achieve this, the following points have been explicitly addressed and discussed:

1. (a)

The description of DROPS tomography was rewritten using the language of quantum information processing instead of the language of NMR. A simple example is the consistent use of Pauli operators $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$ instead of the corresponding spin operators $I_{x}$, $I_{y}$, and $I_{z}$, which are not only unfamiliar outside of the NMR community but also differ by a factor of two. Another example is the description of quantum operations in terms of elementary quantum gates using the QASM [Cross et al.(2022)Cross, Javadi-Abhari, Alexander, De Beaudrap, Bishop, Heidel, Ryan, Sivarajah, Smolin, Gambetta, and Johnson] nomenclature instead of their descriptions in terms of rotations and pulse sequences. 2. (b)

Whereas it is natural to measure expectation values directly on ensemble quantum processors, in pure-state quantum processors, expectation values are typically estimated by averaging the outcomes of projective measurements for many repetitions of the experiment. The measurement of expectation values of single-qubit and multi-qubit operators necessary for DROPS tomography is explicitly discussed. 3. (c)

The fact that DROPS tomography of unitary processes requires an ancilla qubit to be prepared in the completely mixed state could create the false impression that it cannot be applied to pure-state quantum processors. We removed this hurdle by explaining and demonstrating how the concept of temporal averaging can be used to circumvent this problem. 4. (d)

We also showed how to implement discretized scalar products between droplet functions defined on a finite number of sampling points and how to use them to extract the standard matrix representation and the fidelity of states and processes based on experimentally measured droplets. Furthermore, we presented the results of a numerical study of the effect of different sampling schemes on the fidelity with which states can be experimentally tomographed. 5. (e)

Finally, for a convenient adaption of the presented approaches, we provided the Python package DROPStomo [Devra et al.()Devra, Glaser, Huber, and Glaser] for a direct implementation using Qiskit, which can also be adapted to other frameworks in a straightforward way.

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie-Sklodowska-Curie grant agreement No 765267 (QuSCo). S.G. acknowledges funding by the German Research Foundation (DFG) under Germany’s Excellence Strategy – EXC-2111 – 390814868. D.H. acknowledges support from the Verband der chemischen Industrie e.V (VCI). The project is part of the Munich Quantum Valley (MQV) initiative, which is supported by the Bavarian state government with funds from the Hightech Agenda Bayern Plus. We thank Frederik vom Ende for his useful comments on the manuscript. We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team. The quantum circuits presented in this paper were prepared using the latex package Q-circuit [Eastin and Flammia(2004)].

S-I Visualization of operators

A method for representing and visualizing arbitrary quantum operators such as density matrices, Hamiltonians, unitary processes, etc., was proposed in [Garon et al.(2015)Garon, Zeier, and Glaser], which is also known as DROPS representation (Discrete Representation of OPeratorS). Although the DROPS representation can be applied to arbitrary finite-dimensional quantum systems [Garon et al.(2015)Garon, Zeier, and Glaser, Leiner et al.(2020)Leiner, Zeier, and Glaser], here we focus on systems consisting of one or more qubits which are of particular interest in quantum information processing.

The method maps an arbitrary operator $A$ on a set of spherical functions $f^{(\ell)}$ which can easily be visualized, e.g. in the polar representation used here (or alternatively as colored spherical surfaces) [Leiner et al.(2017)Leiner, Zeier, and Glaser]. Panel (a) of Fig. S1 shows an example where the density operator $A$ of a two-qubit state is mapped to a set of four individual droplet functions $f^{(\ell)}$.

For systems consisting of up to two qubits, each of the individual droplet functions ($f^{\{\emptyset\}}$, $f^{\{1\}}$, $f^{\{2\}}$, and $f^{\{12\}}$) simply corresponds to one of the four possible (sub) sets of qubits: $\{\emptyset\}$ labels the empty set, $\{1\}$ labels the subset consisting only of the first qubit, $\{2\}$ labels the subset consisting only of the second qubit, and $\{12\}$ labels the set consisting of both qubits. Each subset is specified by a corresponding superscript label $\ell$. (For systems consisting of more than two qubits, specifying the (sub) systems of qubits is not sufficient and additional selection criteria, such as permutation symmetry, are necessary to specify each droplet [Garon et al.(2015)Garon, Zeier, and Glaser].) The method first decomposes any given operator $A$ as a sum of operators $A^{(\ell)}$, which are defined based on the criteria of the labels $\ell$ specified above. For example, the operator $A^{\{1\}}$ acts only on the first qubit, whereas, e.g., the operator $A^{\{12\}}$ acts on both the first and the second qubit etc.

As shown on the left side of the double-headed arrow in Eq. S1 (vide infra), each droplet operator $A^{(\ell)}$ can be expanded in an operator basis consisting of irreducible spherical tensor operator components $T^{(\ell)}_{jm}$ with expansion coefficients $c_{jm}^{(\ell)}$ [Garon et al.(2015)Garon, Zeier, and Glaser, Leiner et al.(2020)Leiner, Zeier, and Glaser], where $j$ is the rank and $m$ is the order of the spherical tensor operators $T^{(\ell)}_{jm}$. Based on the well-known correspondence [Silver(1976)] between irreducible tensor operator components $T_{jm}$ and spherical harmonics $Y_{jm}$, a bijective mapping between the droplet operators $A^{(\ell)}$ and corresponding spherical droplet functions $f^{(\ell)}$ can be defined as

[TABLE]

where identical expansion coefficients $c_{jm}^{(\ell)}$ are used on both sides of the double-headed arrow. Eq. S1 can be rewritten in the more compact form

[TABLE]

where we defined the rank-$j$ droplet operators $A^{(\ell)}_{j}$ and droplet functions $f^{(\ell)}_{j}$ as

[TABLE]

As shown for the example of the droplet $f^{\{12\}}$ in panel (b) of Fig. S1, each spherical droplet function $f^{(\ell)}$ can be expressed as a sum of spherical functions $f^{(\ell)}_{j}$ with different ranks $j$. For the case of $f^{\{12\}}$, the rank $j$ can be 0, 1, or 2. (Note that in the special case of the operator A represented in Fig. S1, the droplet function with rank 0 happens to vanish, i.e., $f^{\{12\}}_{0}=0$.) This graphically illustrates the decomposition given on the right side of the double-headed arrow in Eq. S2.

The decomposition of $f^{(\ell)}_{j}$ in terms of the spherical harmonics $Y_{jm}$ with the expansion coefficients $c_{jm}^{(\ell)}$ (see right hand side of Eq. S3) is illustrated in panel (c) of Fig. S1 for the rank $j=1$ spherical function $f^{\{12\}}_{1}$, which is decomposed in terms of the spherical harmonics $Y_{1,-1}$, $Y_{1,0}$, and $Y_{1,1}$ with the expansion coefficients $c_{1,-1}=\dfrac{i}{2}$, $c_{1,0}=0$, and $Y_{1,1}=\dfrac{i}{2}$.

S-II Scalar product for tensor operators and spherical functions

The scalar product between two tensor operators $T_{A}$ and $T_{B}$ is defined as

[TABLE]

where $T_{A}^{\dagger}$ is the adjoint (conjugate transpose) of the operator $T_{A}$. If the tensor operators $T_{A}$ and $T_{B}$ are mapped to spherical droplet functions $f_{A}(\theta,\phi)$ and $f_{B}(\theta,\phi)$ using Eq. S1, the scalar product of Eq. S4 is by construction identical to the scalar product of the droplet functions defined as [Leiner et al.(2017)Leiner, Zeier, and Glaser]

[TABLE]

This definition corresponds to the following simple procedure: for each point on the surface of a sphere, the complex conjugate value of the spherical function $f_{A}$ is multiplied by the value of the spherical function $f_{B}$ and the integration of the resulting product $f_{A}^{*}f_{B}$ over the surface of the sphere is the value of the desired scalar product $\langle{f_{A}|f_{B}}\rangle$.

If the spherical functions are only known at a finite number of sample points ($\theta_{i}$,$\phi_{i}$) the scalar product can be approximated by a corresponding discretized scalar product of the form

[TABLE]

Here, the weights $\text{w}_{i}$ reflect the relative surface area represented by each sample point, which depends on the distribution of the sample points on the sphere for a chosen sampling scheme.

As an illustrative example, let us consider the simple case of an equiangular grid, where the sample points are equally distributed along polar ($\theta\in[0,\pi]$) angles as:

[TABLE]

for $k=1,2,\dots,M+1$. The points are distributed for azimuthal ($\phi\in[0,2\pi)$) angles as:

[TABLE]

for $l=1,2,\dots,2M$, where $M$ is a constant and $d=\dfrac{\pi}{M}$ is the angle increment. In this case Eq. S6 can be written in the form

[TABLE]

At the north pole (where $k=1$, corresponding to $\theta_{1}=0$) and at the south pole (where $k=M+1$, corresponding to $\theta_{M+1}=\pi$), the weights are

[TABLE]

For all other sampling points, the weights are given by

[TABLE]

where $\theta_{k}$ is given in Eq. S7. In Eq. S10 and Eq. S11, $l$ runs from 1 to 2$M$ and in Eq. S11, $k$ runs from 2 to $M$.

In the experiments, for redundancy, for each polar angle $\theta_{k}$, we measured not $2M$ but $2M+1$ phase angles $\phi_{l}$ for sampling points, i.e., according to Eq. S8 the phase angle $\phi_{2M+1}=2\pi$ and hence corresponds to $\phi_{1}=0$. This can be simply taken into account by scaling the weights by half for these specific points, i.e., $\text{w}_{k,1}=\text{w}_{k,2M+1}=\frac{1}{2}\text{w}_{k,l}$, where $\text{w}_{k,l}$ is given in Eq. S10 and Eq. S11.

S-III Estimation of expectation values

A general single-qubit state is given by,

[TABLE]

where the probability $p_{0}=|c_{0}|^{2}$ to find the state $|0\rangle$ is given by the expectation value $\langle P_{0}\rangle=\langle\psi|P_{0}|\psi\rangle=|c_{0}|^{2}$ of the projection operator $P_{0}=|0\rangle\langle 0|$ and the probability $p_{1}=|c_{1}|^{2}$ to find the state in $|1\rangle$ is given by the expectation value $\langle P_{1}\rangle=\langle\psi|P_{1}|\psi\rangle=|c_{1}|^{2}$ of the projection operator $P_{1}=|1\rangle\langle 1|$, where the matrix forms of the projection operators $P_{0}$ and $P_{1}$ are

[TABLE]

As the Pauli operator $\sigma_{z}$ can be written as

[TABLE]

the expectation value of $\sigma_{z}$ is given by

[TABLE]

Similarly, the expectation value of the identity matrix $\mathds{1}$ is

[TABLE]

Experimentally, the probabilities $p_{0}$ and $p_{1}$ can be estimated based on the number of shots $N_{0}$ and $N_{1}$ for which the states $|0\rangle$ and $|1\rangle$ are found: $p_{0}\approx N_{0}/N_{s}$ and $p_{1}\approx N_{1}/N_{s}$, where $N_{s}=N_{0}+N_{1}$ is the total number of shots per sampling point.

For a two-qubit system ($N=2$), a general quantum state has the form

[TABLE]

and the expectation value of the operators $\mathds{1}^{[2]}=\mathds{1}^{[1]}\otimes\mathds{1}^{[1]}$, $\sigma_{1z}^{[2]}=\sigma_{z}^{[1]}\otimes\mathds{1}^{[1]}$, $\sigma_{2z}^{[2]}=\mathds{1}^{[1]}\otimes\sigma_{z}^{[1]}$, and $(\sigma_{1z}\sigma_{2z})^{[2]}$ are given by

[TABLE]

where $P_{ab}=|ab\rangle\langle ab|$ and $p_{ab}$ for $a,b\in\{0,1\}$ are the projection operators and probabilities corresponding to the computational basis states $|ab\rangle$. Hence, the required expectation values for the presented Wigner tomography scheme can be calculated by combining the estimated probabilities $p_{ab}$ to find the system in the computational basis states. Experimentally, the probabilities $p_{00}$, $p_{01}$, $p_{10}$, and $p_{11}$ can be estimated based on the number of shots $N_{00}$, $N_{01}$, $N_{10}$, and $N_{11}$ for which the states $|00\rangle$, $|01\rangle$, $|10\rangle$, and $|11\rangle$ are found: $p_{00}\approx N_{00}/N_{s}$, $p_{01}\approx N_{01}/N_{s}$, $p_{10}\approx N_{10}/N_{s}$, and $p_{11}\approx N_{11}/N_{s}$, where $N_{s}=N_{00}+N_{01}+N_{10}+N_{11}$ is the total number of shots.

Note that other expectation values, such as $\langle\sigma_{1x}\rangle$, $\langle\sigma_{1x}\sigma_{2x}\rangle$, etc. can be obtained by applying appropriate local unitary operations $u_{n}$ (so-called detection-associated rotation, see Sec. 3.2) before the projective measurement of the computational basis states.

S-IV Generic rotation matrix

The state of a qubit can be transformed using unitary rotation operations. The matrix form of the most-general single-qubit rotation, which is also known as $\mathrm{U}_{3}$ gate, is (up to a global phase) given by

[TABLE]

This can also be written in the form of the following Euler angle decomposition

[TABLE]

This corresponds to a rotation of angle $\lambda$ around the $z$ axis, followed by a rotation $\theta$ around the $y$ axis, and followed by a rotation of $\phi$ around the $z$ axis. The values of the three Euler angles $\theta,\phi,\lambda$ can be adjusted to implement any desired unitary operation or a gate [Goldstein et al.(2002)Goldstein, Poole, and Safko]. For example, $\mathrm{U}_{3}(\theta,-\pi/2,\pi/2)$ corresponds to a rotation around the $x$ axis by an angle $\theta$ and $\mathrm{U}_{3}(\theta,0,0)$ corresponds to a rotation around the $y$ axis by an angle $\theta$. The rotation operation $R_{\alpha\beta}$ used for scanning in the state tomography algorithm corresponds in this nomenclature to $\mathrm{U}_{3}(\beta,\alpha,0)$ for $\beta\in[0,\pi]$ and $\alpha\in[0,2\pi]$. Similarly, different steps of the presented Wigner tomography algorithm can be designed with the help of the general rotation gate $\mathrm{U}_{3}$.

S-V Quantum circuit and plots of corresponding expectation values for one-qubit Wigner state tomography

In Fig. S3, we show the quantum circuit for performing state tomography of a single-qubit considering the example shown in Fig. 2b in the main text, where the state of the qubit is the superposition state $|\psi\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$. As shown in Fig. S3, this state is prepared from the initial state $|\psi_{i}\rangle=|0\rangle$ using $\mathrm{U}_{3}(\frac{\pi}{2},0,0)$ (alternately, a Hadamard gate could be used). This preparation step is followed by the rotation step $\mathcal{R}$ which is implemented by $\mathrm{U}_{3}^{-1}(\beta,\alpha,0)=\mathrm{U}_{3}(-\beta,0,-\alpha)$.

In Fig. S3, we plot corresponding simulated and experimental expectation values $\langle\sigma_{z}\rangle_{\tilde{\rho}^{[1]}}$ for all combinations of scanning angles $\beta_{k}=(k-1)\dfrac{\pi}{7}$ and $\alpha_{l}=(l-1)\dfrac{\pi}{7}$, where $k=1,2,\dots,8$ and $l=1,2,\dots,15$. Fig. S3a shows the ideal expectation values with no noise, Fig. S3b shows the simulated expectation values with shot noise for number of shots $N_{s}=8192$, and Fig. S3c shows the experimental expectation values for the same number of shots $N_{s}$. Based on these expectation values, the rank $j=1$ droplet function $f_{1}^{(1)}$ can be calculated using Eq. 19 and Fig. S3 shows also plots of the corresponding tomographed droplets.

S-VI Generating basis droplets

Here, we summarize the explicit form of the ideal spherical droplet functions $f_{\sigma_{k}}(\beta,\alpha)$ in terms of linear combinations of complex spherical harmonics $Y_{jm}(\beta,\alpha)$. This decomposition is based on the general basis transformations from Pauli to tensor operators, c.f. Ref. [Garon et al.(2015)Garon, Zeier, and Glaser]. As shown in Sec. 4.1 and Sec. 4.2, these basis droplet functions $f_{\sigma_{k}}(\beta,\alpha)$ can be used to estimate the density matrix corresponding to an experimentally measured droplet.

S-VI.I One qubit

The basis droplet functions ($f_{\sigma_{k}}$) for $k\in\{0,1,2,3\}$ in terms of complex spherical harmonics ($Y_{jm}$) for one qubit are:

[TABLE]

where $\sigma_{0}$ corresponds to the $2\times 2$ identity operator.

S-VI.II Two qubits

The basis droplet functions ($f_{\sigma_{k}}$) for $k\in\{0,1,\dots,15\}$ in terms of complex spherical harmonics ($Y_{jm}$) for two qubit are:

[TABLE]

where $\sigma_{0}$ corresponds to the $4\times 4$ identity operator.

S-VII Temporal averaging to create maximally mixed states

For Wigner quantum process tomography, the preparation of the system qubit in a maximally mixed state is required (see preparation ($\mathcal{P}$) step in Sec. 5.2). As stated in Sec 6.1, in our experimental implementation, a temporal averaging approach [Knill et al.(1998)Knill, Chuang, and Laflamme, Preskill(1998)] was used to prepare the maximally mixed state. A detailed explanation of this approach and alternative methods are provided in the following.

S-VII.I Creating maximally mixed state from pure states

In general, there is an infinite number of different ways to prepare the “maximally mixed state” $\rho_{mm}$ of a single-qubit ($N=1$) based on an ensemble of pure single-qubit states [Preskill(1998)]. For example, the maximally mixed state can be realized using an ensemble of two extremal pure states $|{\uparrow_{\hat{n}}}\rangle$ and $|{\downarrow_{\hat{n}}}\rangle$ corresponding to two antipodal points on the Bloch sphere:

[TABLE]

Hence, for the single-qubit case, the preparation of a maximally mixed state $\rho_{mm}^{[1]}$ (see Eq. S23) can be realized by simply repeating the experiment twice and averaging the measurement results (which we refer to as temporal averaging here): once by starting with the pure-state $|0\rangle$ (corresponding to the density operator |0\rangle\langle 0|=\bigl{(}\begin{smallmatrix}1&0\\ 0&0\end{smallmatrix}\bigr{)}) and once by starting with the pure-state $|1\rangle$ (corresponding to the density operator |1\rangle\langle 1|=\bigl{(}\begin{smallmatrix}0&0\\ 0&1\end{smallmatrix}\bigr{)}). Note that the average density operator of the two experiments is identical to the desired maximally mixed state $\rho_{mm}^{[1]}$.

An alternative way to create a maximally mixed state is by creating a maximal entanglement between the system and an ancilla qubit and partially tracing out of the ancilla qubit. However, this approach requires additional resources (both in a number of qubits and gates) compared to the scheme presented here.

S-VII.II Creating the maximally mixed state of two-qubit by temporal averaging

For the case of a two-qubit ($N=2$) system, the maximally mixed state $\rho_{mm}^{[2]}$ is of the form

[TABLE]

In this case, the preparation of $\rho_{mm}^{[2]}$ can be realized by simply repeating the experiment four times (for all the computational basis states) and averaging the measurement results.

S-VII.III Creating the state $|+\rangle\langle+|\otimes\frac{1}{2^{N}}\mathds{1}^{[N]}$ required for process tomography of unitary gates

Here we give a detailed calculation of the preparation step $\mathcal{{P}}$ (see Sec. 5) for one and two-qubit Wigner process tomography. For single-qubit process tomography, two qubits $q_{0}$ and $q_{1}$ are required and the preparation step requires preparing qubit $q_{0}$ in an equal superposition state $|+\rangle=\dfrac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ and $q_{1}$ in a maximally mixed state $\dfrac{1}{2}\mathds{1}^{[1]}$. The corresponding density matrix is given by

[TABLE]

This maximally mixed state of qubit $q_{1}$ is achieved using deterministic temporal averaging (see Eq. S23), i.e., by averaging measurement result of two experiments initialized in states $|\psi_{1}\rangle=|+\rangle\otimes|0\rangle=|+0\rangle$ and $|\psi_{2}\rangle=|+\rangle\otimes|1\rangle=|+1\rangle$. The two corresponding density matrices are

[TABLE]

and

[TABLE]

Therefore, the density matrix $\rho_{0}^{[2]}$ can be prepared by averaging the density operators $\rho_{1}^{[2]}$ and $\rho_{2}^{[2]}$:

[TABLE]

Similarly, for a two-qubit process tomography, three qubits $q_{0}$, $q_{1}$, and $q_{2}$ are required. The preparation step requires preparing qubit $q_{0}$ in an equal superposition state $|+\rangle$ and qubits $q_{1}$ and $q_{2}$ in a maximally mixed state $\dfrac{1}{4}\mathds{1}^{[2]}$. This resultant density matrix $\rho_{0}^{[3]}$ can be achieved using Eq. S24 as follows:

[TABLE]

This corresponds to performing four different experiments and averaging the output density operators.

S-VII.IV Expectation values for mixed states

Here, we explicitly show that the expectation values for a mixed state can be expressed as an average of the expectation values of pure states. The expectation value of an observable $A$ is given by

[TABLE]

The density matrix $\rho$ can be written in its spectral decomposition as

[TABLE]

where $p_{k}$ are non-negative coefficients which add up to one and $|\beta_{k}\rangle$ are the basis states. Eq. S30 can be rewritten as:

[TABLE]

Hence, for the special case where $\rho$ is the completely mixed state with identical coefficients $p_{k}=\frac{1}{2^{N}}$, the expectation value $\langle{A}\rangle_{\rho}$ is simply the average of the expectation values for the pure states $|\beta_{k}\rangle$. For example, for a single qubit the maximally mixed state $\mathds{1}^{[1]}/2$ can be written as: where $|\beta_{1}\rangle=|0\rangle$, $|\beta_{2}\rangle=|1\rangle$, and $p_{1}=p_{2}=\frac{1}{2}$. This implies that experimentally estimating the expectation value of an observable $A$ for the mixed state $\frac{\mathds{1}^{[1]}}{2}$ can be realized by averaging the expectation values from two different experiments, one in which the qubit is in the pure state $|0\rangle$, and another in which the qubit is in the pure state $|1\rangle$.

S-VIII Additional figures for Wigner state and process tomography results

Whereas in Figs. 3 and 8 of the main text all bilinear terms $f_{j}^{\{12\}}$ with $j\in\{0,1,2\}$ are merged into a single droplet $f^{\{12\}}$, for completeness, here we also plot the individual droplets $f_{j}^{\{12\}}$ separately in Figs. S4 and S5. Furthermore, in contrast to Figs. 3 and 8, where the simulated droplets were calculated and plotted with high resolution using a large number of sampling points, in Figs. S4 and S5 we plot the simulated droplet with the same resolution used in the experiments, i.e., eight polar $\beta\in\{0,\frac{\pi}{7},\cdots\pi\}$ and fifteen azimuthal $\alpha\in\{0,\frac{2\pi}{14},\cdots 2\pi\}$ angles. For the same resolution, simulated and experimental droplet functions show comparable plotting artifacts of the Matlab display function. Similarly, Fig. S6 shows the individual droplets $f_{0}^{(\emptyset)}$ and $f_{1}^{(1)}$ for the process tomography results shown in Fig. 4 of the main text.

S-IX Proof of equivalence of Eq. 40 and Eq. 41

After the mapping of the process matrix onto the density matrix, the Eq. 39 can be written as:

[TABLE]

where $\sigma^{+}=\frac{1}{2}(\sigma_{x}+i\sigma_{y})$, and $\sigma^{-}=\frac{1}{2}(\sigma_{x}-i\sigma_{y})$. The equivalence of Eq. 41 and Eq. 40 can be shown by the following steps which transform Eq. 41 to Eq. 40:

[TABLE]

Note that axial tensor operators are Hermitian, i.e., $({T_{j,\alpha\beta}^{(\ell)[N]}})^{\dagger}=({T_{j,\alpha\beta}^{(\ell)[N]}})$ [Leiner et al.(2017)Leiner, Zeier, and Glaser].

S-X Plot of the mean fidelity as a function of the total number of shots

Here, we present the plot of the mean fidelities $\bar{\mathcal{F}_{s}}$ for four different sampling schemes and the standard state tomography method as a function of the total number of shots $N_{tot}$. Fig. S7 shows the same data as in Fig. 12 and in addition provides the standard deviations. Fig. S8 shows the same numerical study for an equal superposition state $|\psi\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$.

S-XI Additional figures for visualizing rotation errors

For visualizing the result of rotation errors using the DROPS representation, we consider an example in Fig. 11 where for conciseness only the skyscraper visualization of the real part of the density matrix is shown. For completeness, in Fig. S9 we also show the imaginary part of the skyscraper plot along with decomposed bilinear droplets for the state with and without rotation errors. Fig. S10 shows the droplet plots from a different perspective to emphasize the misalignment errors.