Simulating a quantum circuit with decoherence and noise

Question

Based on the answers given here and here, it is pretty clear that an arbitrary quantum circuit can be simulated with matrix algebra. The difficulty is that this assume perfect fidelity. I am unsure how to generalize this method to take into account imperfect gates and decoherence. I am sure that density matrices are involved.

Answer 1

You can simulate a noisy quantum circuit using density matrix evolution (as opposed to state vector evolution,) adding some noise operators to account for imperfect gates as well.

Some good reads could be Simulating noisy quantum circuits with matrix product density operators and Qiskit textbook chapter on density matrices.

Answer 2

If the goal is to classically simulate a noisy quantum circuit using discrete noise channels, there is an approximate simulation technique that avoids using density matrices (and therefore an exponential overhead) by simulating trajectories of the noisy quantum circuit.

Suppose you are interested in simulating an $m$-qubit unitary $U$ that can be decomposed as $U = \prod_{i=1}^m U_i$. We will assume a noise model where each application of $U_i$ is followed by some channel $\Phi_i$ acting on the state, and that the Kraus decomposition for $\Phi_i$ is known and is normalized to have the form: \begin{equation}\tag{1} \Phi_i (\rho)= \sum_{j_i} p(j_i) K_{j_i} \rho K_{j_i}^\dagger \end{equation}

Then starting from the input state $|\psi\rangle$ and using the shorthand $\mathcal{U}_i(\rho) = U_i \rho U_i^\dagger$ we can compute the output as \begin{align} |\psi'\rangle \langle \psi'| &= \left(\prod_{i=1}^m \mathcal{U}_i \circ \Phi_i \right) |\psi\rangle \langle \psi| \tag{2} \\&= \mathcal{U}_m \circ \Phi_m \circ \cdots \circ\mathcal{U}_1 \circ \Phi_1(|\psi\rangle \langle \psi|) \tag{3} \\&= \sum_{j_1, \dots, j_m} p(j_1)\dots p(j_m)U_m K_{j_m} \cdots U_1 K_{j_1} |\psi\rangle \langle \psi | K_{j_1}^\dagger U_1^\dagger \cdots K_{j_m}^\dagger U_m^\dagger \tag{4} \end{align}

Then the key to trajectory simulation is to sample a vector $\mathbf{j}$ with probability $p(\mathbf{j}) = \prod_{i=1}^m p(j_i)$ given according to the Kraus operators and perform (possibly subnormalized) state vector simulation: \begin{equation} |\phi_\mathbf{j} \rangle = U_m K_{j_m} \cdots U_1 K_{j_1}|\psi\rangle \tag{5} \end{equation}

Then you repeat this $N$ times and average the state vectors together to get \begin{equation} |\tilde{\phi}\rangle \langle \tilde{\phi}| = \frac{1}{N}\sum_{\ell=1}^N |\phi_\mathbf{j}^{(\ell)} \rangle \langle \phi_\mathbf{j}^{(\ell)} | \tag{6} \end{equation} where each $\mathbf{j}$ is drawn independently at iteration $\ell$. But this is just an empirical average for the state resulting from noisy simulation, and so given large enough $N$, this empirical average will approach its expected value (by law of large numbers) which is given as \begin{align} \lim_{N\rightarrow \infty} |\tilde{\phi}\rangle \langle \tilde{\phi}| &= \underset{\mathbf{j}}{\mathbb{E}} [|\phi_\mathbf{j} \rangle \langle \phi_\mathbf{j}|] \tag{7} \\&= \sum_\mathbf{j} p(\mathbf{j})|\phi_\mathbf{j} \rangle \langle \phi_\mathbf{j}| \tag{8} \\&= |\psi'\rangle \langle \psi'| \tag{9} \end{align}

where line (8) is substituted for Equation (4). This means if we perform some large number of size-$2^n$ state vector simulations and average the results, we get an increasingly good approximation for a size $2^{2n}$ density matrix simulation. The risk is that you will not generally know what $N$ is sufficient to get good convergence, and it will depend on your given noise model. But for intermediate-size systems (say $n\geq 20$), the linear scaling in $N$ is very tempting compared to the (worst case naive estimate) $O(2^{3n})$ scaling from matrix multiplication for $n$-qubit density operators.

This is the technique implemented in the qsim circuit simulator, but with a number of optimizations that are described in the white paper (disclaimer: I am a coauthor on that preprint).

Answer 3

In addition to the rich physics answer by @user1271772 in terms of Hamiltonian, and computational approach answer by @3yakuya, I would like to add the quantum information picture here.

The imperfect fidelity can be modeled by two realistic scenarios:

• The noisy gates (due to realistic crosstalk or decoherence via $T_1$ [amplitude damping] and $T_2$ [phase damping])
• The false measurement apparatus (e.g. some measurement returning false outcome with some probability $p$)

The most general case incorporating both these phenomenon is a resulting density matrix which has suffered these noise types, which can be modeled as a quantum channel. It can be expressed using operator sum representation resulting into a final density matrix that can be beautifully written in terms of a superoperator as follows:

\begin{equation} \Pi_{\text{real}}(\rho)=\sum_e \bigg( a_e^ME_e\Pi^M_{\text{ideal}}(\rho)E^\dagger_e+b_e^\bar{M}E_e\Pi^\bar{M}_{\text{ideal}}(\rho)E^\dagger_e \bigg) \end{equation}

Here $\Pi_{\text{real}}(\rho)$ is the real projector incorporating the noise from the two above mentioned sources, with $\rho$ as any multiqubit state in general. $E_e$ are the error operators (usually strings of Pauli) and $M$ denotes the desired operator and $\bar{M}$ denotes the false measurement operator. $a_e^M$ and $b_e^\bar{M}$ are the probabilities of the true and false measurement outcomes. Here everything in known, in principle, given that you know what is the noise for the given circuit. Thus, you get the final density matrix with realistic noise.

You can say that after each single or two qubit gate you have an error $E_e$ in the circuit with some probability, and furthermore you can put a gate with a probability $b_e^\bar{M}$ that flips the measurement out come $M$ for any set of qubits desired.

Note that this does not contain any quantum gate required for the desired computation. Because quantum error correction is repeatedly formed and then computation, followed by desired rounds of error correction again, as the requirement maybe. You are not doing both simultaneously, they will not necessarily commute.

Answer 4

Let's call the Hamiltonian for the circuit $\hat{H}_\textrm{OQS}$, because it's an "open quantum system" which is open to decoherence from its environment, which we'll say has a Hamiltonian denoted $\hat{H}_\textrm{Bath}$. How shall we model $\hat{H}_\textrm{Bath}$? The environment is everything in the universe that interacts with $\hat{H}_\textrm{OQS}$: If it's the $\textrm{N}_2$ or $\textrm{O}_2$ molecules of the air that are interacting with the qubit, the most accurate quantum mechanical Hamiltonians for them involve their kinetic energies and the Morse/Long-range (MLR) potentials between them, which are basically harmonic oscillators until the atoms of the molecule get extremely far apart, for example at high enough temperatures for their bonds to start breaking. Likewise, if the qubits are in some solid-state material like silicon-chip qubits or qubits in GaAs quantum dots, then the atoms of the lattice can vibrate with MLR-like potentials, which are also harmonic oscillators except when the material is about to break. Whether the qubits are interacting with the vibrations of molecules in the air, or the vibrational phonons of a solid-state lattice, or even with photons in the air, $\hat{H}_\textrm{Bath}$ can be modeled as a sum of quantum harmonic oscillators:

\begin{align}\tag{1} \hat{H}_\textrm{Bath} &= \sum_{\kappa=1}^N \left( \hat{H}_{\kappa,\textrm{kinetic}} + \hat{H}_{\kappa,\textrm{potential}} \right)\\ &=\sum_{\kappa=1}^N \left( \frac{1}{2}m_\kappa \hat{\dot{X}}^2 + \frac{1}{2}m_\kappa \omega^2_\kappa \hat{X}^2\right)\tag{2}\\ &=\sum_{\kappa=1}^N \omega_\kappa\hat{a}^\dagger_\kappa \hat{a}_\kappa.\tag{3} \end{align}

The simplest coupling we can have between each state $|j\rangle$ of the mult-qubit system $\hat{H}_\textrm{OQS}$ and the oscillators of the environment, is linear coupling (and whatever the real coupling is mathematically, we can imagine that in a Taylor-type approximation of it, the first non-zero term will be linear):

\begin{align}\tag{4} \hat{H}_{\textrm{interaction},|j\rangle} &= \sum_{\kappa=1}^N c_\kappa |j\rangle\langle j | \hat{X}_\kappa \\ &= \sum_{\kappa=1}^N c_{\kappa,j} |j\rangle \langle j | \left( \hat{a_\kappa} + \hat{a_\kappa}^\dagger\right).\tag{5} \end{align}

Equation 5 is linear in $|j\rangle\langle j |$ and linear in $\hat{X}_\kappa$ and any more complicated Hamiltonian would appear exactly like Eq. 5 in a first-order expansion, which is already much better than simply ignoring decoherence (also note that I've absorbed the mass into the $\hat{X}$ operators to form scaled coordinates $\hat{Q}$ before using second quantization to represent the harmonic oscillators in terms of creation and annihilation operators, with the zero-point energy term ignored since it's proportional to the identity matrix and therefore has no non-trivial effect on the circuit's dynamics).

We now have the Feynman-Vernon Hamiltonian (which was actually studied earlier in a single-author paper about masers by Feynman's student Wells circa 1959-1961, and likely studied with different notation by Bloch and Wangsness before their 1948 paper, and apparently was used earliest in a Russian-language paper by Bogoliubov, also see how the Froehlich Hamiltonian dating back to at least 1950 is presented here):

\begin{align}\tag{6} \hat{H} & = \hat{H}_\textrm{OQS} + \hat{H}_\textrm{Bath} + \hat{H}_\textrm{interaction}\\ &= \hat{H}_\textrm{OQS} + \sum_{\kappa=1}^N \omega_\kappa\hat{a}^\dagger_\kappa \hat{a}_\kappa+ \sum_j^M\sum_{\kappa=1}^N c_{\kappa,j} |j\rangle \langle j | \left( \hat{a_\kappa} + \hat{a_\kappa}^\dagger\right),\tag{7} \end{align}

where $M=2^n$ for an n-qubit circuit.

You seem to already know that a closed quantum system (with no decoherence) with (constant-in-time) Hamiltonian $\hat{H}$ and initial density matrix $\rho$ evolves according to:

$$ \rho(t) = e^{-\frac{\textrm{i}}{\hbar}Ht}\rho(0)e^{\frac{\textrm{i}}{\hbar}Ht}\tag{8}. $$

This is exactly the formula you would use to get the dynamics of the combined system comprised of the circuit plus its environment, so you could do the same matrix algebra to get the dynamics, this time including decoherence/noise. There might be many vibrational modes $N$ though, and many levels of the harmonic oscillators that are important, so the matrices might get huge quickly. Fortunately Feynman figured out how to solve this problem. If there's so many vibrational modes interacting with the system that we can say that the system's degrees of freedom interact with a mode of frequency $\omega$ with strength proportional to $J(\omega)$ (this is called the "spectral distribution function" or "spectral density function") for all non-negative real numbers $\omega$, and we have the initial condition:

\begin{align} \rho(0) &= \rho_{\textrm{OQS}}(0) \otimes \rho_\textrm{Bath}(0)\tag{9} \\ &= \rho_{\textrm{OQS}}(0) \otimes e^{-\beta H_\textrm{Bath}},\tag{10} \end{align}

with $\beta \equiv \frac{1}{k_B T}$ at temperature T and with $k_B$ being the Boltzmann constant, then the dynamics of the density matrix of the circuit is given by the following double Feynman integral (it would be a single Feynman integral if we only had to evolve the wavefunction, but since circuit will get entangled with the bath, we need to propagate a density matrix, which is like propagating two wavefunctions):

\begin{align} \rho_{\textrm{OQS}}(t) &= \textrm{Tr}_\textrm{Bath} \rho(t) \tag{11}\\ &= \int \int \hat{O} \rho(0) I\mathcal{D}j(t) \mathcal{D}j^\prime(t),\tag{12} \end{align}

where $I$ is known as the Feynman-Vernon influence functional (describing the influence of noise/decoherence) and $\hat{O}$ represents the dynamics of the quantum circuit without noise/decoherence.

If there's enough interest, I can give exact formulas for $I$ and $\hat{O}$ but it will involve a lot of typing.

The main point from a practical point of view, is that open source software in MATLAB and Python are available for calculating this double Feynman integral numerically for at least up to 16 qubits (maybe more) for if the bath is mathematically friendly enough for the calculation.

If you need to simulate more qubits or a more mathematically unfriendly bath, or need a faster calculation done but don't mind losing accuracy, Markovian (or even non-Markovian) master equations can be used to simulate $\rho_\textrm{Bath}(t)$ for the same Hamiltonian as above. Some of these master equations can give results quite close to the exact Feynman integral described above, especially if you use a unitary transform to look at the Hamiltonian in a frame which depends on parameters that are variationally optimized to minimize the Feynman-Bogoliubov bound on the free energy (essentially making the circuit-bath coupling weak enough in that frame to use weak-coupling master equation techniques without significant loss of accuracy), then transform back to the laboratory frame.