Quantum simulation of quantum channels in nuclear magnetic resonance

TL;DR

This paper introduces an efficient NMR-based quantum simulation framework for quantum channels, utilizing a decomposition of non-unitary operators into unitary components, demonstrated through three key quantum channels.

Contribution

The paper presents a novel method for simulating quantum channels in NMR by decomposing non-unitary operators into unitary parts with minimal ancillary qubits, and experimentally demonstrates this approach.

Findings

Successfully simulated phase damping, amplitude damping, and depolarizing channels.

Measured fidelity and von Neumann entropy for each channel.

Experimental implementation using nuclear spins in liquid samples.

Abstract

We propose and experimentally demonstrate an efficient framework for the quantum simulation of quantum channels in NMR. Our approach relies on the suitable decomposition of non-unitary operators in a linear combination of $d$ unitary ones, which can be then experimentally implemented with the assistance of a number of ancillary qubits that grows logarithmically in $d$. As a proof-of-principle demonstration, we realize the quantum simulation of three quantum channels for a single-qubit: phase damping (PD), amplitude damping (AD), and depolarizing (DEP) channels. For these paradigmatic cases, we measure key features, such as the fidelity of the initial state and the associated von Neuman entropy for a qubit evolving through these channels. Our experiments are carried out using nuclear spins in a liquid sample and NMR control techniques.

Tables

Table 1. Table 1: The standard deviations between simulated results and theoretical predictions. The subscript in ϵ italic-ϵ \epsilon indicates the input state.

Deviations	${ϵ}_x$	${ϵ}_{\overline{y}}$	${ϵ}_z$
PD channel	0.0136	0.0158	0.0107
AD channel	0.0091	0.0098	0.0103
DEP channel	0.0994	0.0203	0.0434

Equations

$$\begin{array}[]{l}U_{0}=\mathcal{I},U_{1}=\sigma_{z},V=W=\left({\begin{array}[]{*{20}{c}}\sqrt{\frac{1+\sqrt{1-\lambda}}{2}}&\sqrt{\frac{1-\sqrt{1-\lambda}}{2}}\\ \sqrt{\frac{1-\sqrt{1-\lambda}}{2}}&{-\sqrt{\frac{1+\sqrt{1-\lambda}}{2}}}\end{array}}\right).\end{array}$$

$${\mathcal{H}}_{int}^n=\sum _{i=1}^n{\omega }_i{\sigma }_z^i+\sum _{i<j,i=1}^n\frac{\pi J_{ij}}{2}{\sigma }_z^i{\sigma }_z^j,$$

$${\rho }_{\text{thermal}}=\frac{{\mathcal{I}}^{\otimes n}}{2^n}+\sum _{i=1}^n{ϵ}_i{\sigma }_z^i,$$

$$\rho_{0}=\frac{1-\epsilon}{2^{n}}\mathcal{I}+\epsilon|0\rangle\mbox{$\langle 0|$}^{\otimes n}.$$

$$\begin{array}[]{l}U_{0}=\sigma_{x},U_{1}=i\sigma_{y},V=W=\sqrt{\frac{1}{2}}\left({\begin{array}[]{*{20}{c}}1&1\\ 1&-1\end{array}}\right).\end{array}$$

$$\left({\begin{array}[]{*{20}{c}}\sqrt{1-\frac{3p}{4}}&-\sqrt{\frac{\frac{p}{4}(1-\frac{3p}{4})}{1-\frac{p}{4}}}&-\sqrt{\frac{\frac{p}{4}(1-\frac{3p}{4})}{(1-\frac{p}{4})(1-\frac{p}{2})}}&-\sqrt{\frac{p}{4-2p}}\\ \sqrt{\frac{p}{4}}&\sqrt{1-\frac{p}{4}}&0&0\\ \sqrt{\frac{p}{4}}&-\frac{p}{4\sqrt{1-\frac{p}{4}}}&\sqrt{\frac{1-\frac{p}{2}}{1-\frac{p}{4}}}&0\\ \sqrt{\frac{p}{4}}&-\frac{p}{4\sqrt{1-\frac{p}{4}}}&-\frac{p}{4\sqrt{(1-\frac{p}{4})(1-\frac{p}{2})}}&\sqrt{\frac{4-3p}{4-2p}}\end{array}}\right).$$

$${\mathcal{M}}_{x,y,z}={\sigma }_{x,y,z}\otimes {\mathcal{I}}^{\otimes n-1}.$$

$$\begin{array}[]{l}\mathcal{M}_{x}^{m,n}=\sigma_{x}\otimes|b(m-1,n-1)\rangle\mbox{$\langle b(m-1,n-1)|$}\ \ {\rm and}\\ \mathcal{M}_{y}^{m,n}=\sigma_{y}\otimes|b(m-1,n-1)\rangle\mbox{$\langle b(m-1,n-1)|$},\end{array}$$

$${\rho }_{\text{S}}=\frac{1}{2}\mathcal{I}+\frac{⟨{\mathcal{M}}_x⟩}{2^n}{\sigma }_x+\frac{⟨{\mathcal{M}}_y⟩}{2^n}{\sigma }_y+\frac{⟨{\mathcal{M}}_z⟩}{2^n}{\sigma }_z,$$

$$T_m=T_{m-1}+O(n).$$

Full text

I Introduction

In the last decades, significant progress has been achieved in the isolation and coherent control of quantum systems, allowing for the observation of their unitary dynamics Lloyd ; N ; QS1 ; fenggrprl ; QS2 ; fenggrsrep ; QS3 . Such a degree of controllability has resulted in the implementation of quantum machines composed of a growing number of qubits, which have been used for key tests of quantum simulations and quantum computers. As envisioned by Richard Feynman Fey82 , large-scale quantum simulators would open the door to the analysis of new quantum physical phenomena and to the study of various models that are nowadays intractable with classical computers. In opposition to quantum simulators of closed quantum systems, the simulation of open quantum systems, which has also been the subject of some research both from a theoretical SL2001 ; ED2008 ; PRB2017 ; PRA2017 ; ARXIV2017 ; wang2013solovay ; SC ; sweke2014simulation ; DiCandia15 and experimental helu ; N2011 ; NP2013 point of view, has been comparatively less explored. In this sense, both from a theoretical and experimental perspective, simulating open quantum systems pose relevant challenges. For example, understanding how quantum systems interact with their environment could potentially shed light on the physics of photosynthetic processes or transport phenomena in general Huelga13 ; Mostame12 , which in turn could help design more efficient light-harvesting devices Scully11 ; Dorfman11 ; Creatore13 . It could also help understand dissipation and thermalisation processes, or the nature of phase transitions. In the same manner, topics related to the foundations of quantum physics, as the measurement process or the quantum-to-classical transition Zurek03 , would greatly benefit from a deeper physical understanding of open quantum systems.

In this work, we consider the simulation of a general CPTP channel dynamics and provide an efficient quantum algorithm for the implementation of non-unitary quantum dynamics associated to paradigmatic quantum channels. Our approach works by decomposing the non-unitary operators into a linear combination of unitary ones. This can be physically implemented via the assistance of a number of ancillary qubits that scales logarithmically with respect to the number of the involved unitary operators. We experimentally demonstrate our proposed quantum simulation method via the implementation of a set of decoherence quantum channels on a nuclear spin-qubit with NMR control techniques. More specifically, we implement the phase damping (PD), the amplitude damping (AD), and the depolarising (DEP) channels.

II Theoretical Results

An open quantum system can be defined as a subsystem of a larger system that includes the open system and its environment and follows a unitary dynamics, as described by $\rho_{\text{se}}=U(\rho\otimes\rho_{\text{env}})U^{\dagger}$. Here, $\rho$ and $\rho_{\text{env}}$ are the initial states of the system and the environment, respectively, and are considered to be initially uncorrelated. The evolution of the principal system can be retrieved as $\rho_{\text{s}}=\text{tr}_{\text{env}}(U(\rho\otimes\rho_{\text{env}})U^{\dagger})$, where tr${}_{\text{env}}$ is the partial trace over the environment degrees of freedom N . Alternatively, the evolution of the system can also be described by a completely positive and trace-preserving map ruskai2002analysis : $\large\varepsilon(\rho)=\sum_{k}E_{k}\rho E_{k}^{\dagger}$, where $E_{k}$ are Kraus operators satisfying $\sum_{k}E_{k}^{\dagger}E_{k}=I$. Non-unitary processes of open quantum systems can also be described by master equations. While the Kraus formalism provides the description of the dynamics for a discrete time step, a master equation can provide a continuous time evolution of the density matrix that describes the open quantum system.

Our method builds upon the framework of the so-called duality quantum computing (DQC) r1 . Such a framework allows for the arbitrary sum of $d$ unitary operators acting on an $n$-qubit system by the addition of log${}_{2}(d)$ two-level ancillary systems. Considering that Kraus operators {$E_{k}$} can also be decomposed into a linear sum of $d$ unitary operators, DQC appears to be of direct applicability to the simulation of open quantum system. A schematics of our proposal follows these steps:

($a$) A $d$-dimensional ancillary system is added to our working system (for example by the addition of $n=\log_{2}d$ qubits) and the setup is initialised in the state $|\Psi\rangle|0\rangle$, where $|\Psi\rangle$ and $|0\rangle$ are the input states of the working system and the $d$-dimensional ancilla, respectively. One additional operation $V$ is then performed on the auxiliary qudit $|0\rangle$, transforming the system to: $|\Psi\rangle|0\rangle\rightarrow\sum_{i=0}^{d-1}V_{i0}|\Psi\rangle|i\rangle$, where $V_{i0}$ are the first column elements of the unitary matrix $V$ and are determined by the target map {$E_{k}$}.

($b$) The controlled operation U_{c}=U_{0}\otimes|0\rangle\mbox{\langle 0|}+U_{1}\otimes|1\rangle\mbox{\langle 1|}+....+U_{d-1}\otimes|d-1\rangle\mbox{\langle d-1|} is implemented afterwards. Here, $U_{0}$, $U_{1}$, … , $U_{d-1}$ are the unitary basis corresponding to the decomposition of the elements {$E_{k}$}. This will result in the system evolving to the state $\sum_{i=0}^{d-1}V_{i0}U_{i}|\Psi\rangle|i\rangle$.

($c$) Operation $W$ is performed on the auxiliary system, resulting in $\sum_{i}V_{i0}U_{i}|\Psi\rangle W|i\rangle=\sum_{i}\sum_{k}W_{ki}V_{i0}U_{i}|\Psi\rangle|k\rangle$, where $W_{ki}V_{i0}$ are complex coefficients, and the sum $\sum_{i=0}^{d-1}W_{ki}V_{i0}=(WV)_{k0}$ corresponds to the $(k,0)$ element of the unitary matrix $WV$, and therefore satisfies $|\sum_{i=0}^{d-1}W_{ki}V_{i0}|\leq 1$. Thus, given a non-unitary transformation described by {$E_{k}$}, its corresponding evolution can be efficiently implemented if the unitary operations $V$, $W$, and $U_{c}$, satisfying $E_{k}=\sum_{i}W_{ki}V_{i0}U_{i}$, are found. Notice that the first column of $V$ is defined by the specific decomposition of the Kraus operators into unitary operators that is chosen, while the rest of the matrix can be arbitrarily completed, with the only requirement of it being unitary. On the other hand, matrix $W$ is uniquely determined by $V$.

($d$) Finally, measuring the corresponding final state of the working system, with the ancillary system in state |k\rangle\mbox{\langle k|}, will result in E_{k}|\Psi\rangle\mbox{\langle\Psi|}E^{\dagger}_{k}. Therefore, if we trace out from the final state of the complete system, the degrees of freedom associated to the ancillary qubits, that is, if we sum over each state |k\rangle\mbox{\langle k|}, with $\{|k\rangle\}$ a complete basis of the ancillary system, the result $\large\varepsilon(\rho)=\sum_{k}E_{k}\rho E_{k}^{\dagger}$, with \rho=|\Psi\rangle\mbox{\langle\Psi|}, will correspond to the simulation of the map {$E_{k}$}.

III The three paradigmatic channels and experiments

PD channel .– We will start the illustration of our method analysing the effect of a PD channel acting on a single-qubit walls1985effect . The effect of the PD channel is to remove the coherences of the qubit stored in the non diagonal elements of its density matrix $\rho_{\text{in}}$. In the Kraus representation, this corresponds to $E_{0}=[1\leavevmode\nobreak\ 0;0\leavevmode\nobreak\ \sqrt{1-\lambda}]$ and $E_{1}=[0\leavevmode\nobreak\ 0;0\leavevmode\nobreak\ \sqrt{\lambda}]$, where the parameter $\lambda\in[0,1]$ represents the strength of the PD channel. In Fig. 1(a) we give the quantum circuit that would realize such a noise channel according to the method introduced in this paper, which needs the addition of a single ancillary qubit. For this case, Kraus operators $E_{0}$ and $E_{1}$ can be decomposed into a linear combination of the unitary operators $\mathcal{I}$ and $\sigma_{z}$, where $\mathcal{I}$ is a $2\times 2$ identity matrix and $\sigma_{x,y,z}$ are Pauli matrices. The decomposition is given by $E_{0}=\frac{1+\sqrt{1-\lambda}}{2}\mathcal{I}+\frac{1-\sqrt{1-\lambda}}{2}\sigma_{z}$ and $E_{1}=\frac{\sqrt{\lambda}}{2}\mathcal{I}-\frac{\sqrt{\lambda}}{2}\sigma_{z}$.

It can be easily checked that the unitary operators $V$, $W$, $U_{0}$ and $U_{1}$ that fulfil conditions $E_{k}=\sum_{i=0}^{1}W_{ki}V_{i0}U_{i}$ $(k=0,1)$ for a PD channel are given by

[TABLE]

As illustrated in Fig. 1(a), the composite system consisting of an ancillary qubit and a working qubit is initialised in state \rho^{\text{CH}}_{\text{in}}=\rho_{\text{in}}\otimes|0\rangle\mbox{\langle 0|}, with the input state of the working qubit \rho_{\text{in}}=|\phi\rangle\mbox{\langle\phi|}. In order to extract the evolution corresponding to the PD channel acting on the working qubit, we need to trace out the ancillary degrees of freedom from the final state $\rho_{\text{out}}^{\text{CH}}$. After doing so, the final state of the working qubit should correspond to $\rho_{\text{out}}=\large\varepsilon^{\text{PD}}(\rho_{\text{in}})=E_{0}\rho_{\text{in}}E_{0}^{\dagger}+E_{1}\rho_{\text{in}}E_{1}^{\dagger}$. The subspace where the ancillary qubit is in the state $|0\rangle$ will be associated with the evolution of the working system that corresponds to $E_{0}\rho_{\text{in}}E_{0}^{\dagger}$, while the subspace of the ancilla state $|1\rangle$ will be associated to $E_{1}\rho_{\text{in}}E_{1}^{\dagger}$.

In order to experimentally demonstrate our proposed quantum simulation scheme, we make use of the nuclear spins in a sample of 13C-labeled chloroform dissolved in deuterated acetone that we manipulate through techniques of NMRluxin2016 ; cory2000nmr . The nuclear spins of 13C and 1H are used to encode the two-level working qubit and the ancillary qubit, respectively. The corresponding molecule structure and parameters are illustrated in Fig. 2(a). Under the weak coupling approximation, the natural Hamiltonian of an $n$-qubit NMR system can be expressed as

[TABLE]

where $\omega_{i}$ is the chemical shift of the $i$th nucleus and $J_{ij}$ is the $J$-coupling constant between the $i$th and the $j$th nuclear spins.

In experiments, starting from an initial thermal equilibrium state, we first generate a pseudo-pure state (PPS) associated to the state $|0\rangle^{\otimes n}$, as the thermal state is not useful for quantum computation because it is a highly mixed state. For the employed liquid sample, the thermal equilibrium state can be written as

[TABLE]

where $n$ is the number of qubits, and $\epsilon_{i}$ represents the polarization of the $i$-th nucleus at room temperature. The spatial averaging technique was used to initilize our system xin2016measurement ; knill1998effective ; cory1997ensemble , taking the thermal state to the following PPS

[TABLE]

A state of this form is convenient as the term related to the identity does not evolve under any unitary propagator and cannot be observed in NMR. Therefore, we can restrict our analysis to the deviation term |0\rangle\mbox{\langle 0|}^{\otimes n} and use it to encode the behaviour of the quantum system. For our experimental analysis, we consider the following initial state for the working qubit: $|X\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$, $|-Y\rangle=(|0\rangle-i|1\rangle)/\sqrt{2}$, and $|Z\rangle=|0\rangle$. Secondly, for each given input state $\rho_{\text{in}}$, we measure the expectation values <$\sigma_{x,y,z}$>=$\text{Tr}(\rho_{\text{out}}\sigma_{x,y,z})$ on the working qubit at the output of the circuit, after it has undergone all the unitary steps. We do this for a collection of values of the parameter $\lambda$, ranging from [math] to $1$ and incremental steps of 1/20. The output state of the working qubit $\rho_{\text{out}}$ is directly obtained via single-qubit tomography.

The whole process, from the PPS |00\rangle\mbox{\langle 00|} to the end of the protocol, contains the following steps: a single rotation of the system qubit to prepare its initial state \rho_{\text{in}}\otimes|0\rangle\mbox{\langle 0|} from the PPS |00\rangle\mbox{\langle 00|}, the operations $V$, $W$, all the controlled operations, and the readout $\pi/2$ pulse. These operations are all packed up together and realised via the GRadient Ascent Pulse Engineering (GRAPE) technique Khaneja2005 ; Negrevergne2008 . The GRAPE approach provides a 5ms pulse width and over 99.5% fidelity for the whole package. Analytically, for any input state of the form $\rho_{\text{in}}=0.5\mathcal{I}+\alpha\sigma_{x}+\beta\sigma_{y}+\gamma\sigma_{z}$, the PD channel should result in a final state $\rho^{\text{th}}_{\text{out}}=0.5\mathcal{I}+\alpha\sqrt{1-\lambda}\sigma_{x}+\beta\sqrt{1-\lambda}\sigma_{y}+\gamma\sigma_{z}$. In Fig. 3(a), the expectation values of <$\sigma_{x,y,z}$> are plotted, which agree well with the theoretically expected values. These clearly show that the PD channel reduces all the magnetisation, $M_{x,y}$, in the $xy$ plane, while keeping the magnetisation, $M_{z}$, in the $z$ direction for any input state $\rho_{\text{in}}$.

AD channel .– We move now to analyse the case of the AD channel fanoc2009 , which is characterised by taking every input to a specific state. The AD channel is described in the Kraus representation via the operators $M_{0}=[1\leavevmode\nobreak\ 0;0\leavevmode\nobreak\ \sqrt{1-\lambda}]$ and $M_{1}=[0\leavevmode\nobreak\ \sqrt{\lambda};0\leavevmode\nobreak\ 0]$. Alternatively, the AD process can be represented as $M_{0}\rho_{\text{in}}M_{0}^{\dagger}+\lambda S_{0}\rho_{\text{in}}S_{0}^{\dagger}$, where $S_{0}$ is a Kraus operator corresponding to the completely positive and trace preserving process described by the set of Kraus operators $\{S_{0},S_{1}\}$, with $S_{0}=[0\leavevmode\nobreak\ 1;0\leavevmode\nobreak\ 0]$ and $S_{1}=[0\leavevmode\nobreak\ 0;1\leavevmode\nobreak\ 0]$. For experimental convenience, we choose to implement this second decomposition in terms of $M_{0}$ and $S_{0}$. We do this, on the one hand, because the $M_{0}\rho_{\text{in}}M_{0}^{\dagger}$ part can be directly obtained from the simulation of the PD channel, and on the other hand, because the simulation of $S_{0}$ is specially convenient as it does not depend on parameter $\lambda$, and therefore a single experimental run serves to compute the effect of any value of $\lambda$, clearly reducing the experimental requirements. The evolution associated to Kraus operators $S_{0}$ and $S_{1}$ can easily be given by the operators $V$, $W$, $U_{0}$, and $U_{1}$ taking values

[TABLE]

Thus, the experiment is performed in two steps, corresponding to the quantum circuit shown in Fig. 1(a) with two different settings of the operators $V$, $W$, $U_{0}$, and $U_{1}$. The first setting is chosen according to Eq. (1), and only subspace $|0\rangle$ of the ancillary qubit is measured, which is associated to the transformation $M_{0}\rho_{\text{in}}M_{0}^{\dagger}$. The second setting is that shown in Eq. (5), and we only measure the subspace of the ancillary qubit corresponding to state $|0\rangle$, which leads to the term $S_{0}\rho_{\text{in}}S_{0}^{\dagger}$.

We use the same sample as that of the previous experiment in order to experimentally simulate the dynamics of the AD channel. We follow the same experimental steps as those in the previous section, performing the experiment twice, for two different setting of the quantum gates in the circuit. As an example, for the case in which the expectation value <$\sigma_{y}$> of the final state is measured for the initial state $|X\rangle$, we firstly prepare the initial state \rho^{\text{CH}}_{\text{in}}=|X\rangle\mbox{\langle X|}\otimes|0\rangle\mbox{\langle 0|} from the PPS and drive it following quantum circuit shown in Fig. 1(a), as described by Eq. (1). Then, the observable \sigma_{y}\otimes|0\rangle\mbox{\langle 0|} is measured to provide the $y$-element associated to the evolution $M_{0}\rho_{\text{in}}M_{0}^{\dagger}$. Next, the same preparation and measurement are performed but this time utilizing the setting of unitary operators in Eq. (5), the results corresponding now to the $y$-element associated to $S_{0}\rho_{\text{in}}S_{0}^{\dagger}$. Combining these two results, one obtains the desired value <$\sigma_{y}$> for a qubit undergoing an AD channel. As in the previous experiment, the GRAPE technique is employed to generate the evolution corresponding to the quantum circuits.

For the input state $\rho_{\text{in}}$, the AD channel can be shown to result in $\rho^{\text{th}}_{\text{out}}=0.5\mathcal{I}+\alpha\sqrt{1-\lambda}\sigma_{x}+\beta\sqrt{1-\lambda}\sigma_{y}+(\gamma(1-\lambda)+0.5\lambda)\sigma_{z}$. In Fig. 3(b), we show the experimental measurement of $\langle\sigma_{x,y,z}\rangle$, necessary for the reconstruction of the system qubit, and how these measurement compare to the analytically computed values. The experimental results show a good agreement with the theoretical predictions, which clearly show that the AD channel damps the system towards the ground state |0\rangle\mbox{\langle 0|}, reducing the magnetisation in the $xy$ plane, while increasing it in the $z$ direction. This could be of interest in the initialisation of a system that is in an arbitrary state.

DEP channel .– To complete our study of decoherence channels, we consider the DEP channel $\large\varepsilon^{\text{DEP}}$ Bennett1999 . For this case, the Kraus representation is given by $E_{0}=\sqrt{1-\frac{3p}{4}}\mathcal{I}$, $E_{1}=\sqrt{\frac{p}{4}}\sigma_{x}$, $E_{2}=\sqrt{\frac{p}{4}}\sigma_{y}$, and $E_{3}=\sqrt{\frac{p}{4}}\sigma_{z}$. It can be trivially shown that the effect of the DEP channel on an initial state $\rho_{\text{in}}$ is to evolve towards the maximally mixed state $\large\varepsilon^{\text{DEP}}(\rho_{\text{in}})=p\mathcal{I}/2+(1-p)\rho_{\text{in}}$, with some probability $p$.

In Fig. 1(b), the quantum circuit to realise the quantum simulation of the DEP channel is depicted, following our proposed protocol. In this case, we set $U_{0}=\mathcal{I}$, $U_{1}=\sigma_{x}$, $U_{2}=\sigma_{y}$, and $U_{3}=\sigma_{z}$, while operator $V$ is given by

[TABLE]

On the other hand, operation $W$ is fixed to a $4\times 4$ identity matrix. The quantum circuit for the DEP channel is then implemented by combining two ancillary qubits initially prepared in the state |00\rangle\mbox{\langle 00|} and a system qubit in input state $\rho_{\text{in}}$. At the end of the protocol, the two ancillary qubits are traced out to acquire the output state of the DEP channel $\rho_{\text{out}}=\large\varepsilon^{\text{DEP}}(\rho_{\text{in}})$.

Experimentally, we need a three-qubit quantum-information processor, which is implemented via diethyl fluoromalonate dissolved in d6 acetone in NMR, where the nuclear spins of 13C, 1H, and 19F in the diethyl fluoromalonate molecule act as the system qubit and the two ancillary qubits, respectively. Figure 2(b) shows the corresponding structure and parameters. The spatial averaging technique is again used to prepare the PPS |000\rangle\mbox{\langle 000|} xine2015 . The unitary operators are implemented via the GRAPE technique that provides a 10ms pulse width. For the DEP channel, we only carry out $xx$-, $\bar{y}y$-, and $zz$-experiments, which are enough to demonstrate the properties of the DEP channel. For instance, for the case of the $\bar{y}y$-experiment, we prepare the system in the state |-Y\rangle\mbox{\langle-Y|}\otimes|00\rangle\mbox{\langle 00|} by applying a $\pi/2$ pulse around the $x$ axis to the system qubit 13C. Then, we measure the observable $\sigma_{y}\otimes\mathcal{I}\otimes\mathcal{I}$ which provides us with the expectation value <$\sigma_{y}$>=Tr$(\rho_{\text{out}}\sigma_{y})$. $\rho_{\text{out}}$ is the output state of the system qubit 13C after tracing out the ancillary qubits 1H and 19F at the end of quantum circuit. Two other experiments are performed in a similar fashion similarly corresponding to $\bar{y}y$-experiments. Fig. (4) illustrates the corresponding results of three experiments $xx$, $\bar{y}y$, and $zz$ for different values of $p$, which presents a good agreement between the theoretical predictions and the experiments.

IV Discussion

In order to evaluate the accuracy of our simulations, we have computed the fidelity $F(\rho^{\text{exp}}_{\text{out}},\rho^{\text{th}}_{\text{out}})$ between the reconstructed single-qubit density matrix $\rho^{\text{exp}}_{\text{out}}$ and the ideal state $\rho^{\text{th}}_{\text{out}}$ for each of the PD and AD channels using the following procedure. These results are illustrated in Fig. 5. The adopted fidelity definition through the whole work is $F=\text{Tr}(\rho_{a}\rho_{b})/\sqrt{\text{Tr}(\rho_{a}^{2})\text{Tr}(\rho_{b}^{2})}$. In our experiments, the average fidelities between the reconstructed single-qubit state $\rho^{\text{exp}}_{\text{out}}$ and the ideal output state $\rho^{\text{th}}_{\text{out}}$ are around 99.52% and 99.87% for the PD and AD channels, respectively.

For an $n$-qubit quantum system $\rho_{\text{SA}}$, with one system qubit and $(n-1)$ ancillary qubits, the operation of tracing out the ancillary qubits, in order to obtain the state of our system qubit $\rho_{\text{S}}=\text{Tr}_{\text{A}}(\rho_{\text{SA}})$, can be realised by measuring the following operators of the output state $\rho_{\text{SA}}$,

[TABLE]

In an NMR platform, $2^{n-1}$ peaks will be observed, with the $m$-th peak providing the expectation values of operators

[TABLE]

where $b(m-1,n-1)$ is the binary representation of number $m-1$ in $n-1$ bits.

Summing the following results over $m$ from 1 to $2^{n-1}$ leads to $\mathcal{M}_{x,y}=\sum_{m=1}^{2^{n-1}}\mathcal{M}_{x,y}^{m,n}$. To measure the observable $\mathcal{M}_{z}$, we apply an additional readout pulse ($\pi/2$ pulse around $y$ axes) on the system qubit at the end, which transfers the magnetization in the $z$ direction to the $x$ direction. In this manner, the expectation value of $\mathcal{M}_{x}$ corresponds the value of the desired observable $\mathcal{M}_{z}$. Moreover, single-qubit tomography of the system qubit can easily be realised using the following rule,

[TABLE]

where the coefficient $2^{n}$ is a normalisation constant, and $\left\langle\mathcal{M}_{x,y,z}\right\rangle$ is the expectation value of observable $\mathcal{M}_{x,y,z}$, $\text{Tr}(\rho_{\text{SA}}\mathcal{M}_{x,y,z})$.

We further complete our analysis by the study of the behaviour of some additional properties under these quantum channels. More precisely, we look at the fidelity $F(\rho_{\text{out}},\rho_{\text{in}})$ and the von Neumann entropy $S(\rho_{\text{out}})=-\text{Tr}(\rho_{\text{out}}\text{log}_{2}\rho_{\text{out}})$, for an input state $|X\rangle$ as it undergoes the PD and AD channels. These results are illustrated in Fig. (6). $F(\rho_{\text{out}},\rho_{\text{in}})$ reflects the strength of the quantum channel acting on a qubit, which decreases for higher strength $\lambda$. $S(\rho_{\text{out}})$ quantifies the strength of the entanglement between the system qubit and the ancillary system, such that $S$ will increase together with the strength $\lambda$ of the PD channel. On the opposite, for the AD channel the entropy clearly shows a maximum for $\lambda=0.5$, while it vanishes for minimal ($\lambda=0$) and maximal ($\lambda=1$) values of $\lambda$. This happens because under the PD channel the initial state $|X\rangle$ will tend towards a maximally mixed state while, for an AD channel, it will gradually tend towards the ground state $|0\rangle$ through intermediate mixed states, respectively.

Conclusion

We have proposed a method for the quantum simulation of open quantum dynamics and experimentally implemented the proposed simulation procedure, realizing proof-of-principle experiments in an NMR setup. Our experiment is a small-scale demonstration of the working principles of the proposed techniques, which can be considered as building blocks for more involved protocols. The experimental results show a high degree of correspondence with the theoretical predictions, showing the capacity of our method to simulate paradigmatic decoherence channels. A natural extension of this work is the development of methods to construct algorithms to simulate the dynamics of open quantum systems in higher dimensions. For example, using the Weyl operator basis, any three-dimensional channel in the Kraus representation can be decomposed into a linear combination of Weyl operators: $M=\sum_{n,m=0}^{2}b_{nm}U_{nm}$ welyl , where $U_{nm}$ is the Weyl operator and $b_{nm}$ are coefficients. Namely, we can perform any Kraus operator in the form of a linear combination of Weyl operators with the proposed method.

Acknowledgements.

T. X. , S. W. and G. L. are grateful to the following funding sources: National Natural Science Foundation of China under Grants No. 11175094 and No. 91221205; National Basic Research Program of China under Grant No. 2015CB921002. J. S. P. and E. S. acknowledge financial support from grants: Spanish MINECO/FEDER FIS2015-69983-P and Basque Government IT986-16.

Appendix A: Scaling of the Protocol

The size of the ancillary system in our protocol is given by the greatest of these two: the number of Kraus operators $d_{1}$, and the number of unitary operators $d_{2}$ onto which the Kraus operators are decomposed. For an $n$-qubit system, with a Hilbert space dimension $d_{S}=2^{n}$, any operator can be decomposed in the Weyl basis welyl as the complex superposition of maximally $d_{S}^{2}$ unitary operators, which are also traceless, and trace-wise orthogonal. Therefore, all the simulated Kraus operators, which act on a system of $n$ qubits, can be decomposed into a basis of not more than $N=2^{2n}$ unitary operators. As a consequence, the total number of ancillary qubits is upper bounded by $\log_{2}(N)=2n$. This is similar to other simulation approaches, like for example the Stinespring dilatation method, which also takes a maximum of $2n$ ancillary qubits.

In order to count the number of required gates, we split our protocol in two parts. On the one hand, we have the initial and final operations $V$ and $W$, which act on the ancillary system and are, in general, arbitrary matrices. It is known that an arbitrary unitary operation acting on an $M$-qubit system can always be implemented with a circuit containing a total of $O(M^{3}2^{2M})$ single qubit and CNOT gates N ; complex1 ; complex2 . Therefore, in the most unfavourable case, where a total of $2n$ ancillary qubits are required, our method would employ up to $O(8n^{3}2^{4n})$ single-qubit and two-qubit gates to implement the $V$ and $W$ operations.

On the other hand, we have the controlled unitary operations acting on the target system of $n$ qubits. This operations are not arbitrary, but they correspond to a specific basis of unitary operators. One can, for example, choose a basis consisting of the tensor product of Pauli operators.

In this case, it can be shown that the gate complexity for each of the controlled operations goes like $O(n^{2})$ complex1 . Let $C_{m}(U)$ denote a controlled gate where the number of control qubits is $m$ and $U$ acts on a target system of $n$ qubits. We will use $T_{m}$ to denote the gate cost of decomposing $C_{m}(U)$. The circuit in Fig. 7 shows a suitable decomposition of $C_{2n}(U)$. Moreover, the Toffolli gate over $2n$ qubits can be decomposed into $O(n)$ single-qubit and CNOT gates. On the other hand, $M^{\dagger}$ and $M$, which fulfil $M^{2}=U$, can also be decomposed into $n$ single-qubit gates, as $U$ is a tensor product of Pauli matrices. The cost of decomposing the gates $C_{1}(M)$ and $C_{1}(M^{\dagger})$ is therefore $O(n)$. From such a decomposition, the following recurrence relation can be inferred,

[TABLE]

Hence, the total gate complexity to implement each controlled unitary operator $C_{2n}(U)$ is proportional to $O(n^{2})$.

According to the discussion above, a total of $2^{2n}$ unitary operations form a complete basis of the $n$-qubit system, and therefore the implementation of these basis operators controlled with respect to the ancillary system takes a total of $O(n^{2}2^{2n})$ single qubit and CNOT gates. In total our algorithm in the most general case can be associated to a gate complexity of $O(8n^{3}2^{4n}+n^{2}2^{2n})$. If we consider other simulation methods, like the Stinespring dilatation, were the system is enlarged to accommodate a $2n$-qubit environment in the most general case, we find that to perform an arbitrary unitary operation in the enlarged Hilbert space, we need up to $O(27n^{3}2^{6n})$ single qubit and CNOT gates, which is an exponential factor more gates than we need. This is because, while our method keeps the size of the operations either to the size of the system or to that of the ancillary space, the Stinespring dilatation method needs to perform operations on the complete system-plus-ancilla space. However, for specific cases where the number of Kraus operators is small and their decomposition requires a large number of unitary operators, it can be the case that the Stinespring dilatation method is more convenient. Ultimately, the comparison should be done case by case.

Appendix B: Measurements in NMR

While NMR spectroscopy is a so-called ensemble weak measurement, which does not collapse the total wave function, expectation values of arbitrary global spin observables can be measured, and with these one can reproduce the outcome of projective measurements, which can be distinguished by the spectra of the NMR ensamble and individually operated on with selective pulses in NMR. In this manner, one can imitate the outcomes of projective measurements and their associated probabilities lee ; auccaise ; laflamme . Besides, in NMR, the measurement of the expectation value of an observable corresponds to the spectroscopy of macroscopic ensembles of quantum spins, which results in a usually significantly precise and stable measurement. Indeed, the precision of the measured data is such that the error bars are typically smaller than the plotted dots, as it is the case for the experimental data presented throughout this paper.

Finally, the minor deviations of the measured data can be associated to imperfections of the PPS initialization, imprecisions of the GRAPE pulses, and dephasing effects caused by decoherence, which are the leading sources of error in our setup. We have numerically simulated the GRAPE pulses including a contrasted decoherence model for our qubits, in order to estimate an error bar for each simulated channel. We compute the standard deviation of our simulated data as $\epsilon=\sqrt{\sum_{i=1}^{M}(x^{i}_{\text{sim}}-x^{i}_{\text{th}})^{2}/(M-1)}$, with $M$ the number of sampling points. In Tab. 1 we give the results for different input states in the PD, AD, and DEP channels.