Decoherence

Introduction #

A statevector created with createQureg can only represent pure quantum states. These are states which contain no classical noise nor uncertainty. We can think of pure states as the outputs of perfect quantum computers which perform precisely the operations and gates we intend of them. But today’s quantum computer are very far from perfect. They make errors, like performing gates slightly wrong, or causing unintended interactions between qubits, or worse yet, interacting with an external environment to induce classical uncertainty upon the state. So-called mixed states cannot be described with statevectors; we must use a more descriptive representation such as a density matrix.

Density Matrices #

Creating a density matrix is just as easy as creating a statevector, using createDensityQureg.

QuESTEnv env = createQuESTEnv();

int numQubits = 10;

// begins in pure state |0><0|
Qureg qureg = createDensityQureg(numQubits, env);

Beware that a density matrix is square larger in memory than the equivalent dimension statevector. While an $n$-qubit statevector has $2^n$ amplitudes, an $n$-qubit density matrix has $2^{2n}$ amplitudes. Hence a $10$ qubit density matrix will cost the same total memory as a $20$ qubit statevector. $$|\psi\rangle = \sum\limits_n^{2^\text{qubits}} \alpha_n |n\rangle \;\;\; \rightarrow \;\;\; \{ \; \alpha_1, \;\;\; \alpha_2, \;\;\; \dots, \;\;\; \alpha_{2^{\text{qubits}}} \; \} $$ $$ \rho = \sum\limits_i^{2^{\text{qubits}}}\sum\limits_j^{2^{\text{qubits}}} \beta_{ij} \;|i\rangle\langle j | \;\;\; \rightarrow \;\;\; \begin{pmatrix} \beta_{1,1} & \beta_{1,2} & \dots & \beta_{1,2^{\text{qubits}}} \\ \beta_{2,1} \\ \vdots & & \ddots\\ \beta_{2^{\text{qubits}},1} \end{pmatrix} $$

Notice createDensityQureg returns a Qureg struct, just like createQureg. This allows our density matrix to be passed to all the initialisation, unitary, gate and operator functions that a statevector can be passed to, although behind the scenes, it will be modified in a more complicated way.

initPlusState(qureg, 0);  // qureg -> |+><+|

hadamard(qureg, 0);       // qureg ->  H qureg H

rotateY(qureg, 5, 0.5);   // qureg -> Ry(-0.5) qureg Ry(0.5)

Where a statevector $|\psi\rangle\,$ is modified by a unitary $U$ as $$ |\psi\rangle \; \rightarrow \; U \; |\psi \rangle, $$ a density matrix $\rho\,$ is being modified as $$ \rho \; \rightarrow \; U \; \rho \; U^\dagger $$

Even calculations are overloaded with their more complex density matrix forms, as we discuss later.

// computes Tr( |0><0| qureg |0><0| )
qreal prob = calcProbOfOutcome(qureg, 0, 0);

This means you can often simply replace createQureg with createDensityQureg to seamlessly switch from a statevector to a density matrix simulation, at a penalty in runtime and memory consumption. Density matrices are even destroyed with destroyQureg, the same function as for statevectors.

Mixing #

So far our register has received only unitaries and so has remained in a pure state. But the power of density matrices is to undergo mixing, and precisely describe quantum states with classical uncertainty. QuEST has many functions for simulating decoherence.

QuEST contains the canonical single-qubit channels, such as amplitude damping, dephasing and homogeneous depolarisation.

• mixDamping
• mixDephasing
• mixDepolarising

These are as trivial to simulate as unitary gates!

// decohere the least significant qubit
qreal prob = 0.1;
mixDamping(qureg, 0, prob);
mixDepolarising(qureg, 0, prob); 

QuEST also has convenient paramaterisations for inhomogeneous single-qubit Pauli channels, and two-qubit correlated dephasing and depolarising channels.

• mixPauli
• mixTwoQubitDephasing
• mixTwoQubitDepolarising

In a single function call…

mixTwoQubitDepolarising(qureg, 1, 2, prob);

we have efficiently modified the density matrix $\rho$ in-place to become $$ (1 – \text{prob}) \, \rho \; + \; \frac{\text{prob}}{15} \left( \sum\limits_{\sigma_1 \in \{ X_1, Y_1, Z_1, I_1 \} } \, \sum\limits_{ \sigma_2 \in \{X_2, Y_2, Z_2, I_,2 \} } \, \sigma_1 \, \sigma_2 \, \rho \, \sigma_1 \, \sigma_2 \right) \; – \; \frac{\text{prob}}{15} \, I_1 \, I_2 \, \rho \, I_2 \, I_1 $$

We can even create probabilistic mixtures of specific pure states by combining density matrices.

• mixDensityMatrix

Qureg rho = createDensityQureg(numQubits, env);
Qureg tmp = createDensityQureg(numQubits, env);
initPlusState(tmp);

tGate(rho, 0);
sGate(tmp, 0);
mixDensityMatrix(rho, 0.7, tmp);

This code left qureg rho in state $$ \rho \;\; = \;\; 0.3 \; T \, |0\rangle\langle 0| \, T^\dagger \;\; + \;\; 0.7 \; S \, |+\rangle\langle +| \, S^\dagger $$

But that’s not all; we can specify and simulate general Kraus maps with any dimension!

• mixKrausMap
• mixTwoQubitKrausMap
• mixMultiQubitKrausMap

These allow us to specify any discrete trace-preserving quantum channel on our register, using arrays of ComplexMatrix2, ComplexMatrix4 or ComplexMatrixN.

int qubit = 0;
int numOps = 2;
ComplexMatrix2 ops[2];

// manually specify a damping channel of 0.1 prob
// (C99 will fill non-specified matrix elements with 0)

ops[0] = {
    .real = {{1, 0},
             {0, sqrt(0.9)},
    .imag = {{0}}}; 

ops[1] = {
    .real = {{0, sqrt(0.1)}},
             {0, 0}},
    .imag = {{0}}};

mixKrausMap(qureg, qubit, ops, numOps);

To effect channels of 3 or more qubits, you must use the ComplexMatrixN type, which can be stored either in dynamic memory, or kept in the local stack (C only), and initialised either element-by-element or in batch (C only).

int numTargs = 3;
int targets[3] = {0, 1, 2};

ComplexMatrixN a = createComplexMatrixN(numTargs);
a.real[0][0] = .1;
a.imag[0][2] = .5;
...

ComplexMatrixN b = createComplexMatrixN(numTargs);
initComplexMatrixN(b, 
    (qreal[8][8]) { {0, 0.1, 0.5, ...}, ...},
    (qreal[8][8]) { ... });

ComplexMatrixN c = getStaticComplexMatrixN(numTargs,
    ( { {0.1, 0.2, ...}, ...} ),
    ( { ... } ) );

int numOps = 3;
ComplexMatrixN ops[] = {a, b, c};

mixMultiQubitKrausMap(qureg, targets, numTargs, ops, numOps);

destroyComplexMatrixN(a);
destroyComplexMatrixN(b);

Calculations #

Now that we have density matrices in mixed states, what calculations can we do with them?

We can firstly ask all the same questions we could of statevectors, like probabilities

• calcProbOfOutcome
• calcProbOfAllOutcomes
• calcTotalProb

For example, the code

qreal a = calcProbOfOutcome(psi, q, x);
qreal b = calcProbOfOutcome(rho, q, x);

respectively calculates $$ a = \Big| \; |{x}\rangle\langle {x}|_{q} \; |\psi\rangle \; \Big| ^2 $$ $$ b = \text{trace}\Big( \; |{x}\rangle\langle {x}|_{q} \;\; \rho \;\; |{x}\rangle\langle {x}|_{q} \; \Big) $$

and both in time $O\left(2^{ \text{numQubits} }\right)$. Depending on the noise in our state, we may find that these quantities differ; that decoherence has effected the outcome statistics of our qubits.

We can compute expectation values in much the same way

• calcExpecDiagonalOp
• calcExpecPauliProd
• calcExpecPauliHamil
• calcExpecPauliSum

which behind the scenes, could compute $$ \langle \psi | \, H \, | \psi \rangle \;\;\;\;\; \text{or} \;\;\;\;\; \text{trace}\left( H \, \rho \right) $$

Density matrices can answer more questions than a statevector, such as through general measures of noise impact.

• calcPurity
• calcFidelity

Qureg psi = createQureg(numQubits, env);
Qureg rho = createDensityQureg(numQubits, env);

for (int r=0; r<10; r++) {

    // evolve both registers via 2nd order Trotter
    applyTrotterCircuit(psi, hamil, dt, 2, 1);
    applyTrotterCircuit(rho, hamil, dt, 2, 1);

    // decohere random qubit of rho
    mixDepolarising(rho, rand()%numQubits, 0.1);
}

// study impact of noise
qreal pur = calcPurity(rho);
qreal fid = calcFidelity(rho, psi);

Denisty matrices can furthermore undergo more advanced measures

• calcDensityInnerProduct
• calcHilbertSchmidtDistance

$$ \text{trace} \left( \; {\rho_1}^\dagger \; \rho_2 \; \right) \;\;\;\;\;\;\text{and}\;\;\;\;\;\; \sqrt{ \; \text{trace}\left[ \; (\rho_1 – \rho_2) \; (\rho_1 – \rho_2 )^\dagger \; \right] } $$