Quantum state tomography via mutually unbiased measurements in driven cavity QED systems

We consider a cavity QED system sketched in figure 1(a), wherein a qubit (two-level atom) shown in figure 1(b) is dispersively coupled to a cavity. In order to implement single-qubit unitary transformations, we only need to add a classical driving field on the qubit, see figure 1(b). Under the rotating-wave approximation, the whole system can be described by the Hamiltonian (${\hslash }=1;$ hereafter the same)

$$\begin{eqnarray}&&H={\omega }_{r}{a}^{\dagger }a+\displaystyle \frac{{\omega }_{a}}{2}{\sigma }_{z}+g({a}^{\dagger }{\sigma }_{-}+a{\sigma }_{+})+\displaystyle \frac{{\rm{\Omega }}}{2}({\sigma }_{+}{{\rm{e}}}^{-{\rm{i}}\nu t}+{\sigma }_{-}{{\rm{e}}}^{{\rm{i}}\nu t}),\end{eqnarray} \tag{ 4 }$$

where ${\omega }_{r}$ is the cavity frequency, ${a}^{\dagger }(a)$ is the creation (annihilation) operator of the cavity photon, ${\omega }_{a}$ is the transition frequency of the qubit with the Pauli operators: ${\sigma }_{-}=| 0\rangle \langle 1|$, ${\sigma }_{+}=| 1\rangle \langle 0|$, and ${\sigma }_{z}=| 1\rangle \langle 1| -| 0\rangle \langle 0|$, g is the coupling strength between the qubit and the cavity, Ω and ν are the amplitude and the frequency of the classical driving field.

Zoom In Zoom Out Reset image size

In a frame rotating at the frequency ν of the classical driving field for both the qubit and the cavity, the Hamiltonian  (4) is changed to

$$\begin{eqnarray}&&\tilde{H}={{\rm{\Delta }}}_{r}{a}^{\dagger }a+\displaystyle \frac{{{\rm{\Delta }}}_{a}}{2}{\sigma }_{z}+g({a}^{\dagger }{\sigma }_{-}+a{\sigma }_{+})+\displaystyle \frac{{\rm{\Omega }}}{2}{\sigma }_{x},\end{eqnarray} \tag{ 5 }$$

with ${{\rm{\Delta }}}_{r}={\omega }_{r}-\nu$ being the frequency detuning of the cavity from the classical driving field, and ${{\rm{\Delta }}}_{a}={\omega }_{a}-\nu$ being the frequency detuning of the qubit from the classical driving field.

Similar to [33], in the dispersive regime (i.e., $| g/{\rm{\Delta }}| \ll 1$ with ${\rm{\Delta }}={\omega }_{a}-{\omega }_{r}$), we perform the unitary transformation $U=\mathrm{exp}[g({a}^{\dagger }{\sigma }_{-}-a{\sigma }_{+})/{\rm{\Delta }}]$ on the Hamiltonian (5) to eliminate the direct qubit-cavity coupling. Using the Hausdorff expansion to second order in the small parameter $g/{\rm{\Delta }}$, the effective Hamiltonian reads as

$$\begin{eqnarray}&&{\tilde{H}}^{\prime }={{\rm{\Delta }}}_{r}{a}^{\dagger }a+\displaystyle \frac{1}{2}{\tilde{{\rm{\Delta }}}}_{a}{\sigma }_{z}+\displaystyle \frac{{\rm{\Omega }}}{2}{\sigma }_{x},\end{eqnarray} \tag{ 6 }$$

where ${\tilde{{\rm{\Delta }}}}_{a}={{\rm{\Delta }}}_{a}+L$ with $L={g}^{2}/{\rm{\Delta }}$. As the average photon number of the cavity $\langle {a}^{\dagger }a\rangle \sim 0$, this Hamiltonian (6) can generate rotations of the qubit about any axis on the Bloch sphere by adjusting the frequency ν of the classical driving field. First, if we adjust $\nu ={\omega }_{a}+L$ so that ${\tilde{{\rm{\Delta }}}}_{a}=0$, the rotations of the qubit around x axis on the Bloch sphere can be generated, that is, the unitary transformation ${U}_{x}(\theta )={{\rm{e}}}^{{\rm{i}}\theta {\sigma }_{x}}$ with $\theta ={\rm{\Omega }}{t}_{x}/2$. Therefore, the required unitary transformation ${U}_{x}\left(\tfrac{7\pi }{4}\right)$ can be implemented with the duration time ${t}_{x}=7\pi /(2{\rm{\Omega }})$. Secondly, if the driving frequency is adjusted as $\nu ={\omega }_{a}+L-{\rm{\Omega }}$, we can realize the required Hadamard gate ${\mathbb{H}}$ with the duration time ${t}_{{\mathbb{H}}}=\pi /(\sqrt{2}{\rm{\Omega }})$. Finally, the combination of ${\mathbb{H}}$ and ${U}_{x}(\theta )$ can generate rotations of the qubit about z axis on the Bloch sphere since ${U}_{z}(\theta )={{\rm{e}}}^{{\rm{i}}\theta {\sigma }_{z}}={\mathbb{H}}{U}_{x}(\theta ){\mathbb{H}}$, and the total duration time is ${t}_{z}=(\sqrt{2}\pi +2\theta )/{\rm{\Omega }}$.

For the readout of the qubit states, we need to apply another classical driving field on one side of the cavity, and then detect its SSTS on the other side, see figure 1(a). The interaction Hamiltonian between the applied driving and the cavity reads as

$$\begin{eqnarray}&&{H}_{d}=\epsilon ({a}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\omega }_{d}t}+a{{\rm{e}}}^{{\rm{i}}{\omega }_{d}t}),\end{eqnarray} \tag{ 7 }$$

where is the time-independent real amplitude and ${\omega }_{d}$ is the frequency of the applied driving. The total Hamiltonian H_T of such a driven cavity-qubit system includes H_d plus the first three terms of equation (4). In the dispersive regime and in a frame rotating at the driving frequency ${\omega }_{d}$ for the cavity, the effective Hamiltonian of the whole system is derived as

$$\begin{eqnarray}&&{H}_{T}^{{\rm{eff}}}=(-{{\rm{\Delta }}}_{{dr}}+L{\sigma }_{z}){a}^{\dagger }a+\displaystyle \frac{{\omega }_{a}+L}{2}{\sigma }_{z}+\epsilon ({a}^{\dagger }+a),\end{eqnarray} \tag{ 8 }$$

where ${{\rm{\Delta }}}_{{dr}}={\omega }_{d}-{\omega }_{r}$ is the detuning between the driving frequency and the cavity frequency.

Under the Born–Markov approximation, the dynamics of the whole system is governed by the master equation [29]

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}t}=-{\rm{i}}[{H}_{T}^{{\rm{eff}}},\rho ]+\kappa { \mathcal D }[a]\rho +{\gamma }_{1}{ \mathcal D }[{\sigma }_{-}]\rho +\displaystyle \frac{{\gamma }_{\phi }}{2}{ \mathcal D }[{\sigma }_{z}]\rho ,\end{eqnarray} \tag{ 9 }$$

where ρ is the density matrix of the system and ${ \mathcal D }[{\mathbb{A}}]\rho ={\mathbb{A}}\rho {{\mathbb{A}}}^{\dagger }-{{\mathbb{A}}}^{\dagger }{\mathbb{A}}\rho /2-\rho {{\mathbb{A}}}^{\dagger }{\mathbb{A}}/2$. Here, ${\gamma }_{1}=1/{T}_{1}$ is the qubit energy decay rate, ${\gamma }_{\phi }$ the qubit dephasing rate, and κ the photon decay rate of the cavity. From the master equation (9), the coupled equations of motion related to the desired quantity $\langle {a}^{\dagger }a\rangle$ are given by

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {a}^{\dagger }a\rangle }{{\rm{d}}t}=-2\epsilon {\rm{Im}}\langle a\rangle -\kappa \langle {a}^{\dagger }a\rangle ,\end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle a\rangle }{{\rm{d}}t}=\lambda \langle a\rangle -{\rm{i}}L\langle a{\sigma }_{z}\rangle -{\rm{i}}\epsilon ,\end{eqnarray} \tag{ 10b }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle a{\sigma }_{z}\rangle }{{\rm{d}}t}=(\lambda -{\gamma }_{1})\langle a{\sigma }_{z}\rangle -({\rm{i}}L+{\gamma }_{1})\langle a\rangle -{\rm{i}}\epsilon \langle {\sigma }_{z}\rangle ,\end{eqnarray} \tag{ 10c }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {\sigma }_{z}\rangle }{{\rm{d}}t}=-{\gamma }_{1}(\langle {\sigma }_{z}\rangle +1).\end{eqnarray} \tag{ 10d }$$

with $\lambda ={\rm{i}}{{\rm{\Delta }}}_{{dr}}-\tfrac{\kappa }{2}$. From equation (10d), we can obtain $\langle {\sigma }_{z}({t}_{{\rm{NDR}}})\rangle ={{\rm{e}}}^{-{\gamma }_{1}{t}_{{\rm{NDR}}}}[\langle {\sigma }_{z}(0)\rangle +1]-1\simeq \langle {\sigma }_{z}(0)\rangle$ since ${t}_{{\rm{NDR}}}\sim 1/\kappa$ and ${\gamma }_{1}\ll \kappa$ [34], where the subscript NDR denotes nondestructive readout. Further, the normalized SSTS ${T}_{\mathrm{ss}}=\langle {a}^{\dagger }a{\rangle }_{\mathrm{ss}}{\left(\tfrac{\kappa }{2}\right)}^{2}$ of the driven cavity is analytically derived from equations (10a)–(10c) as

$$\begin{eqnarray}&&{T}_{\mathrm{ss}}=\displaystyle \frac{\kappa }{2}{\rm{Im}}\left[\displaystyle \frac{L\langle {\sigma }_{z}(0)\rangle -{\rm{i}}(\lambda -{\gamma }_{1})}{{\lambda }^{2}+{L}^{2}-\lambda {\gamma }_{1}-{\rm{i}}L{\gamma }_{1}}\right],\end{eqnarray} \tag{ 11 }$$

where the subscript ss denotes steady state.

From equation (8), it can be seen that the interaction Hamiltonian between the qubit and the cavity is ${H}_{\mathrm{int}}=L{\sigma }_{z}{a}^{\dagger }a$, which commutes with the qubit operator ${\sigma }_{z}$, i.e., $[{H}_{\mathrm{int}},{\sigma }_{z}]=0$. Therefore, the transmission spectra readout method proposed above is of the nondestructive property [29]. The interaction Hamiltonian H_int also indicates that the cavity frequency is shifted by $-L$ (or L), if the qubit is in the computational basis state $| 0\rangle$ (or $| 1\rangle$). Hence, the shifts of the cavity frequency can mark the computational basis states of the qubit. On the other hand, only the incident photon whose frequency is equivalent to one of the state-dependent frequencies ${\omega }_{r}\pm L$ of the cavity can transmit the cavity and then be detected. Physically, the detected probability is exactly the superposed probability of the computational basis state in the qubit state. Therefore, the SSTS T_ss (11) has the feature: the relative height of each transmitted peak marked one of the computational basis states corresponds to its superposed probability in the detected qubit state. This can also be verified through numerical experiments, see section 3.3 as well. This manifest advantage allows us to directly read out all the diagonal elements of the density matrix of the detected qubit state with only one kind of such SSTS.

We now numerically demonstrate the MUBs-QST of qubit states in detail with the presented single-qubit unitary transformations in section 3.1 and the SSTS in section 3.2.

The density matrix of an arbitrary qubit state to be determined can be represented as

$$\begin{eqnarray}&&{\rho }_{2}=\displaystyle \frac{1}{2}\left(I+\displaystyle \sum _{j=1}^{3}{r}_{j}{\sigma }_{j}\right),\end{eqnarray} \tag{ 12 }$$

where r_i's are real parameters, ${\sigma }_{1}={\sigma }_{x}$, ${\sigma }_{2}={\sigma }_{y}$, and ${\sigma }_{3}={\sigma }_{z}$. Without loss of generality, if we choose r₁ = 0.5, r₂ = 0.6, and r₃ = 0.2, the density matrix (12) is specified as

$$\begin{eqnarray}{\rho }_{2}=\left(\begin{array}{cc}0.6 & 0.25-0.3{\rm{i}}\\ 0.25+0.3{\rm{i}} & 0.4\end{array}\right).\end{eqnarray} \tag{ 13 }$$

According to the MUBs-QST theory introduced in section 2, to realize the measurements in the MUBs, we first need to perform single-qubit unitary transformations V₂^k implemented in section 3.1 on the state (13). For ${V}_{2}^{1}=I$, the state ${\rho }_{2}$ remains identical. For ${V}_{2}^{2}={\mathbb{H}}$ and ${V}_{2}^{3}={U}_{x}\left(\tfrac{7\pi }{4}\right)$, the state ${\rho }_{2}$ is transformed as ${\rho }_{2}^{\prime }={V}_{2}^{2}{\rho }_{2}{{V}_{2}^{2}}^{\dagger }$ and ${\rho }_{2}^{\prime\prime }={V}_{2}^{3}{\rho }_{2}{{V}_{2}^{3}}^{\dagger }$, respectively. Then, we project the detected states ${\rho }_{2}$, ${\rho }_{2}^{\prime }$, and ${\rho }_{2}^{\prime\prime }$ onto the computational basis. The projective measurement outcomes can be determined by three kinds of SSTSs proposed in section 3.2. The numerically simulated SSTS T_ss (11) of the driven cavity as a function of the detuning ${{\rm{\Delta }}}_{{dr}}={\omega }_{d}-{\omega }_{r}$ is shown in figure 2, wherein panels (a)–(c) correspond to the detected states ${\rho }_{2}$, ${\rho }_{2}^{\prime }$, and ${\rho }_{2}^{\prime\prime }$, respectively. Here, the parameters are chosen as $({\gamma }_{1},\kappa ,L)=2\pi \times (0.19,1.69,-7.38)$ MHz [35]. As the projective measurement outcomes ${p}_{2}^{(k,l)}$ are exactly the corresponding diagonal elements $| l\rangle \langle l|$ of each detected state, we can thus directly read out all the projective measurement outcomes for each MUB from each subfigure. Specifically, in figure 2(a), the heights of the transmitted peaks marked computational basis states $| 0\rangle$ and $| 1\rangle$ are 0.602 and 0.329, respectively. Hence, the projective measurement outcomes are $({p}_{2}^{(1,0)},{p}_{2}^{(1,1)})=(0.602,0.329)$. Similarly, the projective measurement outcomes $({p}_{2}^{(2,0)},{p}_{2}^{(2,1)})=(0.751,0.207)$ and $({p}_{2}^{(3,0)},{p}_{2}^{(3,1)})=(0.801,0.166)$ can be directly read out from figure 2(b) and (c), respectively. Finally, inserting these projective measurement outcomes and the MUBs for d = 2 [18] into equation (2), we can obtain the reconstructed state normalized as

$$\begin{eqnarray}{\tilde{\rho }}_{2}=\left(\begin{array}{cc}0.6593 & 0.3174-0.3711{\rm{i}}\\ 0.3174+0.3711{\rm{i}} & 0.3407\end{array}\right).\end{eqnarray} \tag{ 14 }$$

However, it is noted that the reconstructed density matrix (14) is unphysical because it has a negative eigenvalue violating the property of positive semidefiniteness of all physical density matrices. To avoid this problem, we employ the commonly used maximum likelihood estimation (MLE) technique [2] to derive a physical density matrix most likely to have returned the numerically simulated results (14). In this manner, a physical density matrix is obtained as

$$\begin{eqnarray}\ {\tilde{\rho }}_{2}^{\prime }=\left(\begin{array}{cc}0.6567 & 0.3087-0.3607{\rm{i}}\\ 0.3087+0.3607{\rm{i}} & 0.3433\end{array}\right).\end{eqnarray} \tag{ 15 }$$

And its fidelity [36] is calculated as

$$\begin{eqnarray}&&F({\rho }_{2},{\tilde{\rho }}_{2}^{\prime })={\rm{Tr}}(\sqrt{\sqrt{{\rho }_{2}}{\tilde{\rho }}_{2}^{\prime }\sqrt{{\rho }_{2}}})=0.952.\end{eqnarray} \tag{ 16 }$$

Zoom In Zoom Out Reset image size

We consider a setup shown in figure 1(a), wherein a qutrit (three-level atom with the cascade configuration, without loss of generality) displayed in figure 1(c) is dispersively coupled to a cavity. The Hamiltonian of the qutrit-cavity system reads as

$$\begin{eqnarray}&&{ \mathcal H }={\omega }_{r}{a}^{\dagger }a+\displaystyle \sum _{j=0}^{2}{\omega }_{j}{{\rm{\Pi }}}_{{jj}}+\displaystyle \sum _{j=0}^{1}{g}_{j}({a}^{\dagger }{{\rm{\Pi }}}_{j,j+1}+a{{\rm{\Pi }}}_{j+1,j}),\end{eqnarray} \tag{ 17 }$$

where ${\omega }_{j}$ is the frequency of the energy level $| j\rangle$, g_j is the coupling strength between the cavity and the transition $| j\rangle \leftrightarrow | j+1\rangle$, and ${{\rm{\Pi }}}_{{{jj}}^{\prime }}=| j\rangle \langle {j}^{\prime }|$.

In the dispersive regime, i.e., $| {g}_{0}/{{\rm{\Delta }}}_{0}^{2},{g}_{1}/{{\rm{\Delta }}}_{1}^{2}| \ll 1$ with the frequency detunings ${{\rm{\Delta }}}_{0}={\omega }_{1}-{\omega }_{0}-{\omega }_{r}$ and ${{\rm{\Delta }}}_{1}={\omega }_{2}-{\omega }_{1}-{\omega }_{r}$, and after adiabatically eliminating the cavity mode, the effective Hamiltonian of the system can be obtained as

$$\begin{eqnarray}&&{{ \mathcal H }}_{{\rm{eff}}}={L}_{1}| 1\rangle \langle 1| +{L}_{2}| 2\rangle \langle 2| \end{eqnarray} \tag{ 18 }$$

with ${L}_{1}={g}_{0}^{2}/{{\rm{\Delta }}}_{0}$ and ${L}_{2}={g}_{1}^{2}/{{\rm{\Delta }}}_{1}$. This Hamiltonian can generate the required phase operation R with the possible duration time ${t}_{R}=-10\pi /(3{L}_{1})=-4\pi /(3{L}_{2})$, and its inverse operator ${R}^{-1}$ with the duration time $2{t}_{R}$.

On the other hand, similar to section 3.1, a classical driving field with the amplitude ${{\rm{\Omega }}}_{1}$ and the frequency ${\nu }_{1}$ is applied between $| 0\rangle$ and $| 1\rangle$. By adjusting the frequency ${\nu }_{1}$, an arbitrary unitary transformation between $| 0\rangle$ and $| 1\rangle$ can be produced, i.e., ${U}^{01}({\theta }_{1},{\phi }_{1})={U}_{z}^{01}\left(\tfrac{{\phi }_{1}}{2}\right){U}_{x}^{01}({\theta }_{1}){U}_{z}^{01}\left(-\tfrac{{\phi }_{1}}{2}\right)$ with the total duration time ${t}_{01}=2({\theta }_{1}+\sqrt{2}\pi +2\pi )/{{\rm{\Omega }}}_{1}$. Likewise, we add another classical driving field with the amplitude ${{\rm{\Omega }}}_{2}$ and the frequency ${\nu }_{2}$ between $| 1\rangle$ and $| 2\rangle$. We can generate an arbitrary unitary transformation between $| 1\rangle$ and $| 2\rangle$, ${U}^{12}({\theta }_{2},{\phi }_{2})={U}_{z}^{12}\left(\tfrac{{\phi }_{2}}{2}\right){U}_{x}^{12}({\theta }_{2}){U}_{z}^{12}(-\tfrac{{\phi }_{2}}{2})$, with the total duration time ${t}_{12}=2({\theta }_{2}+\sqrt{2}\pi +2\pi )/{{\rm{\Omega }}}_{2}$ by adjusting the frequency ${\nu }_{2}$. Since the transition $| 0\rangle \leftrightarrow | 2\rangle$ is forbidden, any unitary transformation between $| 0\rangle$ and $| 2\rangle$ can be constructed by a sequence of ${U}^{01}({\theta }_{1},{\phi }_{1})$ and ${U}^{12}({\theta }_{2},{\phi }_{2})$, e.g., ${U}^{02}({\theta }_{3},{\phi }_{3})={U}_{x}^{01}\left(\tfrac{\pi }{2}\right){U}^{12}({\theta }_{3},{\phi }_{3}){U}_{x}^{01}(\tfrac{\pi }{2})$ with the total duration time ${t}_{02}=2\pi /{{\rm{\Omega }}}_{1}+2({\theta }_{3}+\sqrt{2}\pi +2\pi )/{{\rm{\Omega }}}_{2}$. Furthermore, an arbitrary phase operation ${ \mathcal R }({\theta }_{4},{\phi }_{4})$ can be constructed with a sequence of ${U}^{01}({\theta }_{1},{\phi }_{1})$ and ${U}^{02}({\theta }_{3},{\phi }_{3})$, such as ${ \mathcal R }({\theta }_{4},{\phi }_{4})={{ \mathcal R }}^{01}({\theta }_{4}){{ \mathcal R }}^{02}({\phi }_{4})$ with ${{ \mathcal R }}^{01}({\theta }_{4})={U}^{01}\left(\tfrac{3\pi }{4},\tfrac{\pi }{2}\right){U}_{x}^{01}({\theta }_{4}){U}^{01}(\tfrac{\pi }{4},\tfrac{\pi }{2})$ and ${{ \mathcal R }}^{02}({\phi }_{4})={U}^{02}\left(\tfrac{3\pi }{4},0\right){U}^{02}({\phi }_{4},\tfrac{\pi }{2}){U}^{02}(\tfrac{3\pi }{4},\pi )$. The total duration time is ${t}_{{R}}=2(8\pi +2\sqrt{2}\pi +{\theta }_{4})/{{\rm{\Omega }}}_{1}+(15\pi +6\sqrt{2}\pi +2{\phi }_{4})/{{\rm{\Omega }}}_{2}$. Finally, the required inverse operator of the Fourier transformation can be constructed with the combination of the above unitary transformations, e.g., ${F}^{-1}=-{\rm{i}}{U}^{12}\left(\tfrac{7\pi }{4},-\tfrac{\pi }{6}\right){U}^{01}(\pi +\mathrm{arccos}\tfrac{\sqrt{2}}{\sqrt{3}},-\tfrac{\pi }{6}){U}^{02}(\tfrac{\pi }{4},-\tfrac{5\pi }{6}){ \mathcal R }(\tfrac{11\pi }{6},\tfrac{5\pi }{3})$. The total duration time is ${t}_{{F}^{-1}}=\left(83\pi /3+6\sqrt{2}\pi +2\mathrm{arccos}\tfrac{\sqrt{2}}{\sqrt{3}}\right)/{{\rm{\Omega }}}_{1}+(91\pi /3+10\sqrt{2}\pi )/{{\rm{\Omega }}}_{2}$.

Similar to section 3.2, a classical driving field with the same interaction Hamiltonian (7) is employed to achieve the readout of the qutrit states. The total Hamiltonian of such a driven cavity-qutrit system is ${{ \mathcal H }}_{T}={ \mathcal H }+{H}_{d}$. In the dispersive regime and in a frame rotating at the driving frequency ${\omega }_{d}$ for the cavity, the efficient Hamiltonian of the whole system is derived as

$$\begin{eqnarray}&&{{ \mathcal H }}_{T}^{{\rm{eff}}}=\left(-{{\rm{\Delta }}}_{{dr}}+\displaystyle \sum _{j=0}^{2}{S}_{i}{{\rm{\Pi }}}_{{jj}}\right){a}^{\dagger }a+\displaystyle \sum _{j=0}^{2}({\omega }_{j}+{L}_{j}){{\rm{\Pi }}}_{{jj}}+\epsilon ({a}^{\dagger }+a),\end{eqnarray} \tag{ 19 }$$

where ${L}_{j}=\tfrac{{g}_{j-1}^{2}}{{{\rm{\Delta }}}_{j-1}}$ and ${S}_{j}=\tfrac{{g}_{j-1}^{2}}{{{\rm{\Delta }}}_{j-1}}-\tfrac{{g}_{j}^{2}}{{{\rm{\Delta }}}_{j}}$ are Lamb shift and Stark shift of energy level $| j\rangle$, respectively.

The master equation to describe the dynamics of the whole system is given by [29]

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}t}=-{\rm{i}}[{{ \mathcal H }}_{T}^{{\rm{eff}}},\rho ]+\kappa { \mathcal D }[a]\rho +\displaystyle \sum _{j=1,2}{\gamma }_{1}^{j}{ \mathcal D }[{{\rm{\Pi }}}_{j-1,j}]\rho +\displaystyle \sum _{j=1,2}{\gamma }_{\phi }^{j}{ \mathcal D }[{{\rm{\Pi }}}_{{jj}}]\rho ,\end{eqnarray} \tag{ 20 }$$

where ${\gamma }_{1}^{j}=1/{T}_{1}^{j}$ and ${\gamma }_{\phi }^{j}$ are the energy decay rate and the dephasing rate for energy level $| j\rangle$, respectively. From the master equation (20), we can derive a set of coupled equations of emotion associated with the desired quantity $\langle {a}^{\dagger }a\rangle$, which includes equation (10a) and

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle a\rangle }{{\rm{d}}t}=\lambda \langle a\rangle -{\rm{i}}{S}_{0}\langle {{\rm{\Pi }}}_{00}a\rangle -{\rm{i}}{S}_{1}\langle {{\rm{\Pi }}}_{11}a\rangle -{\rm{i}}{S}_{2}\langle {{\rm{\Pi }}}_{22}a\rangle -{\rm{i}}\epsilon ,\end{eqnarray} \tag{ 21a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {{\rm{\Pi }}}_{00}a\rangle }{{\rm{d}}t}=\lambda \langle {{\rm{\Pi }}}_{00}a\rangle -{\rm{i}}{S}_{0}\langle {{\rm{\Pi }}}_{00}a\rangle +{\gamma }_{1}^{1}\langle {{\rm{\Pi }}}_{11}a\rangle -{\rm{i}}\epsilon \langle {{\rm{\Pi }}}_{00}\rangle ,\end{eqnarray} \tag{ 21b }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {{\rm{\Pi }}}_{11}a\rangle }{{\rm{d}}t}=\lambda \langle {{\rm{\Pi }}}_{11}a\rangle -{\rm{i}}{S}_{1}\langle {{\rm{\Pi }}}_{11}a\rangle -{\gamma }_{1}^{1}\langle {{\rm{\Pi }}}_{11}a\rangle +{\gamma }_{1}^{2}\langle {{\rm{\Pi }}}_{22}a\rangle -{\rm{i}}\epsilon \langle {{\rm{\Pi }}}_{11}\rangle ,\end{eqnarray} \tag{ 21c }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {{\rm{\Pi }}}_{22}a\rangle }{{\rm{d}}t}=\lambda \langle {{\rm{\Pi }}}_{22}a\rangle -{\rm{i}}{S}_{2}\langle {{\rm{\Pi }}}_{22}a\rangle -{\gamma }_{1}^{2}\langle {{\rm{\Pi }}}_{22}a\rangle -{\rm{i}}\epsilon \langle {{\rm{\Pi }}}_{22}\rangle ,\end{eqnarray} \tag{ 21d }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {{\rm{\Pi }}}_{00}\rangle }{{\rm{d}}t}={\gamma }_{1}^{1}\langle {{\rm{\Pi }}}_{11}\rangle ,\end{eqnarray} \tag{ 21e }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {{\rm{\Pi }}}_{11}\rangle }{{\rm{d}}t}=-{\gamma }_{1}^{1}\langle {{\rm{\Pi }}}_{11}\rangle +{\gamma }_{1}^{2}\langle {{\rm{\Pi }}}_{22}\rangle ,\end{eqnarray} \tag{ 21f }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {{\rm{\Pi }}}_{22}\rangle }{{\rm{d}}t}=-{\gamma }_{1}^{2}\langle {{\rm{\Pi }}}_{22}\rangle .\end{eqnarray} \tag{ 21g }$$

Similar to the qubit case, from equations (21e)–(21g), we can obtain $\langle {{\rm{\Pi }}}_{00}\rangle \simeq \langle {{\rm{\Pi }}}_{00}(0)\rangle$, $\langle {{\rm{\Pi }}}_{11}\rangle \simeq \langle {{\rm{\Pi }}}_{11}(0)\rangle$, and $\langle {{\rm{\Pi }}}_{22}\rangle \simeq \langle {{\rm{\Pi }}}_{22}(0)\rangle$. Further, the normalized SSTS of the driven cavity can be analytically derived from equations (10a) and (21a)–(21d) as

$$\begin{eqnarray}{T}_{\mathrm{ss}} & = & -\displaystyle \frac{\kappa }{2}{\rm{Im}}\left[\displaystyle \frac{{\rm{i}}}{\lambda }+\displaystyle \frac{{S}_{0}}{\lambda }\left(\displaystyle \frac{\langle {{\rm{\Pi }}}_{00}(0)\rangle }{A}+\displaystyle \frac{{\gamma }_{1}^{1}\langle {{\rm{\Pi }}}_{11}(0)\rangle }{{AB}}+\displaystyle \frac{{\gamma }_{1}^{1}{\gamma }_{1}^{2}\langle {{\rm{\Pi }}}_{22}(0)\rangle }{{ABC}}\right)\right.\\ & & +\left.\displaystyle \frac{{S}_{1}}{\lambda }\left(\displaystyle \frac{\langle {{\rm{\Pi }}}_{11}(0)\rangle }{B}+\displaystyle \frac{{\gamma }_{1}^{2}\langle {{\rm{\Pi }}}_{22}(0)\rangle }{{BC}}\right)+\displaystyle \frac{{S}_{2}}{\lambda }\displaystyle \frac{\langle {{\rm{\Pi }}}_{22}(0)\rangle }{C}\right],\end{eqnarray} \tag{ 22 }$$

with $A={\rm{i}}{S}_{0}-\lambda$, $B={\rm{i}}{S}_{1}-\lambda +{\gamma }_{1}^{1}$, and $C={\rm{i}}{S}_{2}-\lambda +{\gamma }_{1}^{2}$.

The physical mechanism of the SSTS readout method is similar to the qubit case explained in section 3.2. This is also a kind of nondestructive measurement due to the fact that the interaction Hamiltonian between the qutrit and the cavity (i.e., Stark shift term in equation (19)) ${{ \mathcal H }}_{\mathrm{int}}={S}_{j}{{\rm{\Pi }}}_{{jj}}{a}^{\dagger }a$ commutes with the qutrit operator ${{\rm{\Pi }}}_{{jj}}$, that is, $[{{ \mathcal H }}_{\mathrm{int}},{{\rm{\Pi }}}_{{jj}}]=0$. Theoretical analysis and numerical experiments also indicate that multiple peaks emerge in the SSTS (22) with the same feature as the qubit case. This feature is also verified in section 4.3.

With the implemented single-qutrit unitary transformations in section 4.1 and the presented SSTS in section 4.2, we numerically show in detail how to realize the MUBs-QST of qutrit states.

An arbitrary qutrit state can be represented as

$$\begin{eqnarray}&&{\rho }_{3}=\displaystyle \frac{1}{3}\left(I+\displaystyle \sum _{j=1}^{8}{r}_{j}{\alpha }_{j}\right),\end{eqnarray} \tag{ 23 }$$

where r_j's are real parameters, and ${\alpha }_{j}$'s are SU(3) generators [37]. Without loss of generality, if we choose r₁ = 0.3, r₂ = 0.24, r₃ = 0.3, r₄ = 0.36, r₅ = 0.3, r₆ = 0.42, r₇ = 0.36, r₈ = 0.1, the qutrit state (23) is specified as

$$\begin{eqnarray}{\rho }_{3}=\left(\begin{array}{ccc}0.453 & 0.1-0.08{\rm{i}} & 0.12-0.1{\rm{i}}\\ 0.1+0.08{\rm{i}} & 0.253 & 0.14-0.12{\rm{i}}\\ 0.12+0.1{\rm{i}} & 0.14+0.12{\rm{i}} & 0.294\end{array}\right).\end{eqnarray} \tag{ 24 }$$

In order to realize the measurements in the MUBs, we first perform single-qutrit unitary transformations V₃^k implemented in section 4.1 on the qutrit state (24). The qutrit state (24) remains the same for ${V}_{3}^{1}=I$. For the other unitary transformations ${V}_{3}^{2}$, ${V}_{3}^{3}$, and ${V}_{3}^{4}$, the transformed states are given by ${\rho }_{3}^{\prime }={V}_{3}^{2}{\rho }_{3}{{V}_{3}^{2}}^{\dagger }$, ${\rho }_{3}^{\prime\prime }={V}_{3}^{3}{\rho }_{3}{{V}_{3}^{3}}^{\dagger }$, and ${\rho }_{3}^{\prime\prime\prime }={V}_{3}^{4}{\rho }_{3}{{V}_{3}^{4}}^{\dagger }$, respectively. Then, these detected states ${\rho }_{3}$, ${\rho }_{3}^{\prime }$, ${\rho }_{3}^{\prime\prime }$, and ${\rho }_{3}^{\prime\prime\prime }$ are projected onto the computational basis. The projective measurement outcomes can be read out by four kinds of SSTSs presented in section 4.2. In figure 3, we numerically show the SSTS T_ss (22) of the driven cavity as a function of the detuning ${{\rm{\Delta }}}_{{dr}}={\omega }_{d}-{\omega }_{r}$. Panels (a)–(d) correspond to the detected states ${\rho }_{3}$, ${\rho }_{3}^{\prime }$, ${\rho }_{3}^{\prime\prime }$, and ${\rho }_{3}^{\prime\prime\prime }$, respectively. The parameters are choose as $({\gamma }_{1}^{1},{\gamma }_{1}^{2},\kappa )=2\pi \times (0.199,0.227,1.69)$ MHz [38] and ${S}_{(\mathrm{0,1,2})}=4\pi \times (10.0,5.9,3.4)$ MHz (twice of the Stark shifts in [38] to increase the distinguishability of the transmitted peaks). From the position and the height of the transmitted peaks, we can directly read out all the diagonal elements of the detected states, which are exactly the projective measurement outcomes. Specifically, the projective measurement outcomes $({p}_{3}^{(1,0)},{p}_{3}^{(1,1)},{p}_{3}^{(1,2)})=(0.457,0.222,0.244)$ can be directly read out from figure 3(a). Likewise, we can directly read out the projective measurement outcomes $({p}_{3}^{(2,0)},{p}_{3}^{(2,1)},{p}_{3}^{(2,2)})=(0.578,0.232,0.135)$, $({p}_{3}^{(3,0)},{p}_{3}^{(3,1)},{p}_{3}^{(3,2)})=(0.463,0.302,0.161)$, and $({p}_{3}^{(4,0)},{p}_{3}^{(4,1)},{p}_{3}^{(4,2)})=(0.257,0.373,0.262)$ from figure 3(b), (c), and (d), respectively. Finally, substituting these projective measurement outcomes and the MUBs for d = 3 [22] into equation (2), the normalized reconstructed state is obtained as

$$\begin{eqnarray}{\tilde{\rho }}_{3}=\left(\begin{array}{ccc}0.551 & 0.1399-0.1108{\rm{i}} & 0.1618-0.1356{\rm{i}}\\ 0.1399+0.1108{\rm{i}} & 0.2085 & 0.2755-0.1472{\rm{i}}\\ 0.1618+0.1356{\rm{i}} & 0.2755+0.1472{\rm{i}} & 0.2405\end{array}\right).\end{eqnarray} \tag{ 25 }$$

Note that the reconstructed density matrix (25) is also unphysical since it contains a negative eigenvalue. Likewise, this problem is resolved with the commonly used MLE technique [2], which allows us to derive a physical density matrix most likely to produce the numerically simulated results (25). In this way, we can obtain a physical density matrix

$$\begin{eqnarray}{\tilde{\rho }}_{3}^{\prime }=\left(\begin{array}{ccc}0.5267 & 0.1575-0.0913{\rm{i}} & 0.117-0.146{\rm{i}}\\ 0.1575+0.0913{\rm{i}} & 0.2224 & 0.2153-0.0923{\rm{i}}\\ 0.117+0.146{\rm{i}} & 0.2153+0.0923{\rm{i}} & 0.2509\end{array}\right)\end{eqnarray} \tag{ 26 }$$

with the fidelity [36]

$$\begin{eqnarray}&&F({\rho }_{3},{\tilde{\rho }}_{3}^{\prime })={\rm{Tr}}(\sqrt{\sqrt{{\rho }_{3}}{\tilde{\rho }}_{3}^{\prime }\sqrt{{\rho }_{3}}})=0.961.\end{eqnarray} \tag{ 27 }$$

Zoom In Zoom Out Reset image size