1 Introduction

Quantum information processing may potentially revolutionize classical computation based on ordinary bits. However, current constructions of universal quantum computer still cannot compete with classical computational machines. One of the main difficulties lies in the fact that systems on a quantum scale are extremely susceptible to any kind of external noise, as well as to erroneous action of quantum gates in a circuit. Therefore, to handle qubits effectively, there is a need for methods protecting quantum information against all possible disturbances. Two general approaches to this problem have been developed. The first one is based on the so-called decoherence-free subspaces and exploits particular states of the Hilbert space that are immune to certain errors—a readable review of this methods can be found in Refs. [1, 2]. An alternative technique is based on quantum error-correcting codes (QECC), which are quantum counterparts of the classical error-correcting codes. Quantum error-correcting codes are constructs which protect quantum information against some specified errors. This method of error correction has been extensively studied in the case of unitary noise operations—see Ref. [3] for a comprehensive introduction to this field. However, any real quantum operation can be only approximately unitary, and thus, one has to consider non-unitary noise operations. Here, the progress is slower than in the unitary case, mainly because of the increased complexity of the problem. For instance, in the unitary case, any product of two Kraus operators is normal, so its numerical range is determined by the spectrum, which is not longer true in the general case. The need for a constructive method of finding quantum error correction codes for non-unitary noise models provides a motivation toward this work.

The main aim of this paper is to propose a method of finding quantum error correction codes for a class of non-unitary noise operators. The method, in the form that is presented here, can be effectively applied to noise models with short Kraus decomposition (consisting of two operators). The main advantage of this method is that it allows to solve an algebraic problem (often untractable with other approaches) using an elementary geometrical construction. The paper is organized as follows. In Sect. 2, we review some basic notions related to quantum error correction and specify the general form of a non-unitary quantum channel. In Sect. 3, we recall definitions of generalized numerical range, introduce the concept of nuclear numerical range and present some of its properties. In Sect. 4, we describe a geometric method of obtaining quantum error correction codes for block-diagonal Kraus operators. Precisely, using compression formalism based on the Knill–Laflamme conditions [4, 5] and the notion of nuclear numerical range, we obtain projectors on the subspaces of the quantum error correction code. This is the main contribution of this work. In Sect. 5, we provide two examples of block-diagonal channels and obtain the corresponding quantum error correction codes using method described in Sect. 4. We conclude this paper with a summary of results obtained and discuss possible generalizations of the method to higher-dimensional problems.

2 Quantum error correction

In this section, we recall the definition of Kraus operators and their role in description of noise in quantum systems. We invoke the Knill–Laflamme conditions [4] for an error correction code of a particular noise model. In the end of this section, we present a block-diagonal model of quantum noise that will be explored further in this work.

2.1 Kraus representation of a quantum noise

Consider a quantum system \(\rho \) described in an n-dimensional Hilbert space \(\mathcal {H}^{n}\). Assume that the system evolves according to some given error process (noise channel) represented by a superoperator \({\varPhi }\) acting on \(\mathcal {H}^n\). According to the Kraus Representation Theorem [6], any such superoperator can be written as the sum of m matrix operators \(A_i \in \mathcal {M}_n(\mathbb {C})\), where \(\mathcal {M}_n(\mathbb {C})\) is the space of all complex-valued matrices of order n:

$$\begin{aligned} {\varPhi }[\cdot ] = \sum _{i=1}^m A^{\dagger }_i[\cdot ]A_i, \end{aligned}$$
(1)

where \(\dagger \) denotes the adjoint operator. Matrices \(A_i\) are called Kraus operators, and they satisfy the trace-preserving condition: \(\sum _{i=1}^{m} A_i^{\dagger } A_i = \mathbb {I}\). In this paper, we consider (unless otherwise stated) models of noise acting on two-qubit systems described by two Kraus operators. The dimension of matrices \(\{A_i\}\) is thus \(n = 4\), and there are \(m = 2\) of them. To solve the quantum error correction problem for a map \({\varPhi }\), one has to look for subspaces \(\mathcal {H}^n_C \subset \mathcal {H}^n\), which satisfy Knill–Laflamme Conditions [4, 5]. Error-correcting code, labeled by , is itself a quantum state in subspace \(\mathcal {H}^n_C\) of dimension n. We denote by a particular basis of this subspace (correction code basis), so that and by the projection operator on . According to [4], the following conditions are sufficient to reconstruct information about the system \(\rho \) subjected to errors described by the set of m Kraus operators \(\{A_i\}\):

$$\begin{aligned} PA^{\dagger }_iA_jP = \lambda _{ij} P \quad \text {for} \quad i,j = 1, 2, \ldots m. \end{aligned}$$
(2)

where \(\lambda _{ij}\) are called compression values of the error-correcting code. The problem of determining the projectors P is related to an algebraic compression problem [7]. In our case, the invariant subspace of code, \(\mathcal {H}^n_C\), is two-dimensional, so P is a projector on a two-dimensional subspace (\(P = P_2\)). Denoting \(T_{ij} \equiv A^{\dagger }_iA_j\) we can write:

$$\begin{aligned} P_2 T_{ij} P_2 = \lambda _{ij} P_2, \quad i,j = 1, 2, \ldots , m \end{aligned}$$
(3)

Determining quantum error correction code for the error model described by the set of Kraus operators \(\{A_i\}\) is equivalent to finding a subspace \(P_2\) that satisfies the above set of equations.

2.2 Noise models with block-diagonal Kraus operators

The problem we address in this paper involves finding correction subspace \(\mathcal {H}^n_C\), which amounts to determining projections \(P_2\) from the last equation. The main difficulty lies in the fact that \(P_2\) has to satisfy the Knill–Laflamme conditions (3) for all operators \(T_{ij}\) simultaneously. In what follows, we consider the Kraus operators with a block-diagonal structure. It will be convenient to introduce a short-hand notation:

$$\begin{aligned} T_{11}&= A_1^{\dagger } A_1 = E_{11} \oplus F_{11} = \begin{bmatrix} E_{11}&0 \nonumber \\ 0&F_{11} \end{bmatrix}, \\ T_{12}&= A_1^{\dagger } A_2 = E_{12} \oplus F_{12} = \begin{bmatrix} E_{12}&0 \\ 0&F_{12}, \end{bmatrix}, \nonumber \\ T_{22}&= A_2^{\dagger } A_2 = \mathbb {I} - T_{11}, \qquad T_{21} = A_2^{\dagger } A_1 = T_{12}^{\dagger }, \end{aligned}$$
(4)

where \(E_{11}, E_{12}, F_{11}, F_{12} \in \mathcal {M}_2(\mathbb {C})\). Note that if Eq. (3) is fulfilled for \(T_{12}\), it is also fulfilled for the adjoint \(T_{12}^{\dagger } = T_{21}\). Moreover, due to the normalization condition, we can write \(T_{22} = \mathbb {I}-T_{11}\). Thus, we effectively look for simultaneous solutions of the compression problem for two operators \(T_{11}\) and \(T_{12}\). There are several approaches of solving problems of this kind if matrices \(T_{ij}\) are normal (see for example the eigenvector-pairing method introduced in [8]). Here, we consider a broader range of models of non-unitary Kraus operators, for which matrices \(T_{ij}\) need not be normal. Thus, the techniques developed in literature cannot be applied here in a straightforward manner. In Sect. 4, we will develop a new method for solving this type of problems.

3 Mathematical tools

In this section, we review the notion of higher order numerical range and introduce the concept of nuclear numerical range. We state several its properties that will be further explored to determine the correction code subspaces for models considered in Sects. 4 and 5.

3.1 Higher order numerical range

The algebraic form of Eq. (3) suggests that the problem can be approached using the so-called rank- k numerical range \({\varLambda }_k(T)\) of matrix T, introduced in [8],

$$\begin{aligned} {\varLambda }_k(T) = \{ \lambda \in \mathbb {C}: P T P = \lambda P \text { for some } P \in \mathcal {P}_k \}, \end{aligned}$$
(5)

where \(\mathcal {P}_k\) is the set of all rank-k projections on space \(\mathcal {H}^n\). Unit vectors which yield compression value \(\lambda = \lambda _0\) are called generating vectors (or generators) of \(\lambda _0\). The two special cases (\(k = 1\) and \(k = 2\)) are of particular interest in this paper. Notice that setting \(k = 1\) yields the standard numerical range \({\varLambda }_1(T)\), often denoted by W(T) [9]:

(6)

On the other hand, the case \(k = 2\) gives the numerical range of rank two:

$$\begin{aligned} {\varLambda }_2(T) = \left\{ \lambda \in \mathbb {C}: P_2 T P_2 = \lambda P_2,\, P_2 = P_2^{\dagger } \text { and } P_2P_2=P_2 \right\} , \end{aligned}$$

Based on the definition of higher rank numerical range \({\varLambda }_k\), one can derive the following properties [8] :

  1. (P1)

    For any \(a,\, b \in \mathbb {C}\), \({\varLambda }_k(aA+b I) = a {\varLambda }_k(A) + b.\)

  2. (P2)

    For any unitary \(U \in \mathcal {M}_n\), \({\varLambda }_k(U^{\dagger } A U) = {\varLambda }_k(A).\)

  3. (P3)

    If A is Hermitian and \(\{ \lambda _i\}\) is an ordered set of eigenvalues of a matrix A, then \({\varLambda }_k(A) = [\lambda _{n-k+1}(A), \lambda _k(A)]\), where interval is an empty set if \(\lambda _{n-k+1}(A) > \lambda _k(A)\) when \(k > n/2\).

  4. (P4)

    For any \(k_1,\, k_2 \in \mathbb {N}\), \({\varLambda }_{k_1}(A) \cap {\varLambda }_{k_2}(B) \subseteq {\varLambda }_{k_1+k_2}(A\oplus B)\).

Having defined higher order numerical range, we are ready to introduce a new structure, the nuclear numerical range.

3.2 Nuclear numerical range

In this section, we define a structure that will prove useful to determine error correction subspaces.

Definition 1

Given two square, complex-valued matrices of the same size: A, \(Z \in \mathcal {M}_n(\mathbb {C})\), the Z-nuclear numerical range of A, labeled by W(A|Z), is defined as the set of expectation values of A among pure states which belong to the kernel of Z,

(7)

By definition, the nuclear numerical range is a set in the complex plane contained in W(A), that is, \(W(A|Z) \subseteq W(A)\) for all matrices Z. Note that the notion of nuclear numerical range W(A|Z) belongs to the class of restricted numerical ranges described in [10], since the set of states used in Definition 1 is restricted to the set containing the kernel (nucleus) of Z. The standard numerical range W(A) is convex, while the set formed by Z-nuclear numerical range is not convex in general. In this paper, we mostly consider W(A|Z) for matrices of size \(n = 2\). Let us now present some properties of nuclear numerical range.

Let \(X,\,Y\) be two nonzero matrices in \(\mathcal {M}_n(\mathbb {C})\) and U be some unitary matrix in \(\mathcal {M}_n(\mathbb {C})\). Then, the following properties follow directly from the definition of W(A|Z):

  1. (N1)

    \(W(A|0) = W(A)\) is the standard numerical range.

  2. (N2)

    \(W(A|X^{\dagger } X) = \varnothing \) for any non-singular X and all matrices A.

  3. (N3)

    \(W(A|A-\alpha \mathbb {I}) = \{ \alpha \}\), \(\alpha \in W(A)\).

  4. (N4)

    \(W(A|Z) = W(A|\,(-1) \cdot Z)\).

  5. (N5)

    \(W(X|X-Y) = W(Y|Y-X)\).

  6. (N6)

    \(W(A + \alpha Z|Z) = W(A|Z)\).

  7. (N7)

    If Z is traceless, then W(A|Z) is not empty.

Moreover, the following properties hold for two-dimensional matrices (\(n = 2\)). Let us assume that Z is a hermitian matrix, and let \(\nu _1\) and \(\nu _2 > \nu _1\) denote its eigenvalues. Then:

  1. (N8)

    W(A|Z) forms a (possibly degenerate) elliptic curve.

  2. (N9)

    If \({{\mathrm{sgn}}}(\nu _1) \ne {{\mathrm{sgn}}}(\nu _2)\), then W(A|Z) is not empty.

  3. (N10)

    If Z is normal and \(\lambda \in [\nu _1, \nu _2]\), then \(W(A|Z-\lambda \mathbb {I}) = W(A)\).

The proof of (N\(1)-(\)N6) is straightforward from the definition of W(A|Z) and properties (P\(1)-(\)P2) from previous subsection. To prove property (N7) note that W(A|Z) always contains the point \(\frac{1}{2}\hbox {Tr }A\). Thus, \(\frac{1}{2}\hbox {Tr }A = 0 \in W(Z)\) and W(A|Z) is not empty. In order to prove property (N8), we use the following theorem, proved in “Appendix”.

Theorem 1

Let Z be a normal matrix of order two with real-valued entries. Furthermore, let A be an arbitrary complex-valued matrix of order two. Consider:

$$\begin{aligned} Z&= \begin{bmatrix} 2 a&b \\ b&2 c \end{bmatrix} \in \mathcal {M}_2^{\mathbb {R}}, \quad A = \begin{bmatrix} d&f \\ g&h \end{bmatrix} \in \mathcal {M}_2^{\mathbb {C}}. \end{aligned}$$

Then, there exists a set of normalized states , , parametrized by a phase \(\varphi \in [0, 2\pi )\) and a real number \(\lambda \in W(Z)\), which satisfies the following set of simultaneous equations:

(8)

where \(z(\varphi , \lambda )\) forms an elliptic disk in the complex plane parametrized by \(\varphi \) and \(\lambda \):

$$\begin{aligned} z(\varphi , \lambda ) = z_0 + w \lambda + p(\lambda ) [q \cos \varphi + r \sin \varphi ]. \end{aligned}$$
(9)

Variables \(z_0\), w, \(p(\lambda )\), q and r are defined in Eq. (44). The family of states is given by:

(10)

where \(\tan 2\alpha = b / (a-c)\) and \(\cos \theta = (\lambda - a - c) / \sqrt{b^2 + (a-c)^2}\).

From above theorem we can deduce the following simple corollary:

Corollary 1

\(W(A|Z - \lambda \mathbb {I}) = z(\varphi , \lambda )\), where \(\varphi \in [0, 2\pi )\) and \(z(\varphi , \lambda )\) is given by Eq. (9).

The proof of property (N8) follows directly from Theorem 1 if we take \(\lambda = 0\). In order to prove property (N9) note that by the elliptic range theorem [11], the set W(Z) forms an elliptic disk with foci at eigenvalues of Z, that is, \(\nu _1\) and \(\nu _2\). By the convexity property, W(Z) must also contain a line with endpoints \((\nu _1, \nu _2)\). Since \({{\mathrm{sgn}}}(\nu _1) = -{{\mathrm{sgn}}}(\nu _2)\), the line either passes through the origin or is a singular point at the origin. Thus, there exists such that and W(A|Z) is not an empty set. Property (N10) can be proven after noticing that if one allows \(\lambda \) to take all possible real values between the two eigenvalues of Z, then there is no restriction on vector (cosine of azimuthal angle given by Eq. (42) in the Appendix takes all real values between \(-1\) and 1). This means that \(W(A|Z-\lambda \mathbb {I})\) is the full numerical range W(A).

It is worth to emphasize that the unitary invariance of standard numerical range does not hold in the case of nuclear numerical range, that is \(W(A|UZU^{\dagger }) \ne W(A|Z)\). It can be easly seen by considering:

$$\begin{aligned} A = \begin{bmatrix} a&b \\ c&d \end{bmatrix}, \quad Z = \begin{bmatrix} 0&2 \\ 0&0 \end{bmatrix}, \quad U = \frac{1}{\sqrt{2}}\begin{bmatrix} 1&1 \\ 1&-1 \end{bmatrix}. \end{aligned}$$
(11)

Computing W(A|Z) yields a single point \(\lambda = a \in \mathbb {C}\), whereas \(W(A| UZU^{\dagger }) = \lambda ' = (a+b+c+d)/2 \in \mathbb {C}\). Thus, for a general matrix A, \(\lambda \ne \lambda '\). Similary, we can find \(W(UAU^{\dagger }|Z) \ne W(A|Z)\).

4 Determination of code subspaces for Kraus operators with block-diagonal structure

In this section, we describe a method of constructing projectors \(P_2\) onto the code subspace \(\mathcal {H}^n_C\) for not-normal matrices. Let us now return to the basic problem and recall Eq. (3). Because of the block-diagonal structure of \(T_{ij}\), we expect that \(P_2\) will have a similar block-diagonal structure. We emphasize that this choice of \(P_2\) is not the most general possible and reduces the total set of possible correction codes we can obtain. Let us call the upper and lower of the blocks of \(P_2\) by \(P_E\) and \(P_F\), respectively. This allows us to write:

$$\begin{aligned}&P_{2} = P_E \oplus P_F = \begin{bmatrix} P_{E}&0 \\ 0&P_{F} \end{bmatrix}. \end{aligned}$$
(12)

Using property (P4) and setting \(k_1 = k_2 = 2\), we can reduce the initial problem of determining one 4-by-4 projection matrix \(P_2\) to a problem of finding two 2-by-2 projection matrices \(P_E\) and \(P_F\). Writing explicitly:

$$\begin{aligned} {\varLambda }_{2}(T_{ij}) = {\varLambda }_{2} (E_{ij} \oplus F_{ij}) \supseteq {\varLambda }_{1}(E_{ij}) \cap {\varLambda }_{1}(F_{ij}). \end{aligned}$$
(13)

We can rewrite Eq. (3) in the matrix notation:

$$\begin{aligned} P_2T_{ij}P_2 = \begin{bmatrix} P_{E}E_{ij}P_E&0 \\ 0&P_FF_{ij}P_F \end{bmatrix} = \begin{bmatrix} \lambda _{ij}^E P_E&0 \\ 0&\lambda _{ij}^F P_F \end{bmatrix}. \end{aligned}$$
(14)

By definition, \(\lambda _{ij}^E \in {\varLambda }_1(E_{ij})\) and \(\lambda _{ij}^F \in {\varLambda }_1(F_{ij})\) are points in respective standard numerical ranges. The above equality is satisfied only if the condition \(\lambda _{ij}^E = \lambda _{ij}^F = \lambda _{ij}\) holds for all values of i and j. This means that both points, \(\lambda _{ij}^E\) and \(\lambda _{ij}^F\), lie in the intersection \({\varLambda }_{1}(E_{ij}) \cap {\varLambda }_{1}(F_{ij})\) for every choice of ij. In order to determine \(P_E\) and \(P_F\) (which are projections onto points in that intersection), note that we can rewrite them in respective correction code bases as:

(15)

where both vectors are normalized . Note that the states and are both generators of the set of compression values \(\lambda _{ij}\), for \(i, j = 1, 2\). The intersections that we are interested in, \({\varLambda }_{1}(E_{ij}) \cap {\varLambda }_{1}(F_{ij})\), are sets of points for which the following statement holds:

(16)

Due to the structure assumed in Eq. (4) not all of the above equations are independent. Having this in mind, we are left with the following set of two equations:

(17)

To solve this set of equations, we use the concept of nuclear numerical range discussed in Sect. 4. Consider the following two sets:

$$\begin{aligned} W(E_{12}|E_{11} - \lambda _{11} \mathbb {I}) \qquad \text {and} \qquad W(F_{12}|F_{11} - \lambda _{11} \mathbb {I}). \end{aligned}$$

Notice that points in the complex plane which satisfy Eq. (17) are exactly the ones that constitute the intersection of above nuclear numerical ranges. Let us label this intersection by \({\varGamma }(\lambda _{11})\):

$$\begin{aligned} {\varGamma }(\lambda _{11}) = W(E_{12}|E_{11} - \lambda _{11} \mathbb {I}) \cap W(F_{12}|F_{11} - \lambda _{11} \mathbb {I}) \subset \mathbb {C}. \end{aligned}$$
(18)

Clearly, \({\varGamma }(\lambda _{11}) \subset W(E_{12}) \cap W(F_{12})\). Since both matrices \(E_{11}\) and \(F_{11}\) are normal, by property (N8) we conclude that for a given value of \(\lambda _{11}\) these two sets are elliptic curves in the complex plane. If we treat \(\lambda _{11}\) as a parameter whose range is the appropriate line segment \(W(E_{11})\), then by property (N10) we have \(W(E_{12}|E_{11} - \lambda _{11} \mathbb {I}) = W(E_{12})\). Similar statement holds for a pair \(F_{12}\) and \(F_{11}\). In our problem, in order to satisfy Eq. (17), the number \(\lambda _{11}\) must be contained in both \(W(E_{11})\) and \(W(F_{11})\). Thus \(\lambda _{11}\) must be contained in the intersection \(W(E_{11}) \cap W(F_{11})\), which is simply a line segment. If we denote the sets of eigenvalues of matrices \(E_{11}\) and \(F_{11}\) by \(\{\nu _{E}\}\) and \(\{\nu _{F}\}\), respectively, then in order to fulfill the first condition from Eq. (17), one has to satisfy:

$$\begin{aligned} \max \limits _{ \{\nu _{E}\} }> \lambda _{11} > \min \limits _{ \{\nu _{F}\} } \quad \text {or} \quad \min \limits _{ \{\nu _{E}\} }< \lambda _{11} < \max \limits _{ \{\nu _{F}\} }. \end{aligned}$$
(19)

Let us label by \({\varOmega }\) the set of all allowable values of \(\lambda _{11}\). Then, the set \({\varOmega }\) is not empty if and only if:

$$\begin{aligned} \hbox {Tr }E_{11} \pm \left[ \left( \hbox {Tr }E_{11}\right) ^2-4\det E_{11} \right] ^{1/2} \gtrless \hbox {Tr }F_{11} \mp \left[ \left( \hbox {Tr }F_{11}\right) ^2-4\det F_{11} \right] ^{1/2}, \end{aligned}$$
(20)

where the symbol \(\gtrless \) corresponds to \(+\) and − signs respectively, so Eq. (20) contains two inequalities. The above equation follows from the fact that the two eigenvalues of a 2-by-2 matrix A are given by a formula \(\frac{1}{2} \hbox {Tr }A \pm \frac{1}{2}\sqrt{(\hbox {Tr }A)^2 - 4\det A}\). By Corollary 1, we can conclude that \(W(E_{12}|E_{11} - \lambda _{11}) = z(\varphi _E, \lambda _{11})\) and \(W(F_{12}|F_{11} - \lambda _{11}) = z(\varphi _F, \lambda _{11})\), where function \(z(\varphi , \lambda )\) is defined in Eq. (9). The intersection \({\varGamma }(\lambda _{11})\) is determined by:

$$\begin{aligned} {\varGamma }(\lambda ) = \{z \in \mathbb {C}|z = z(\varphi _E, \lambda _{11}) = z(\varphi _F, \lambda _{11}),\ \varphi _E, \varphi _F\in [0, 2 \pi ), \ \lambda _{11} \in {\varOmega }\}. \end{aligned}$$
(21)

In order to determine the vectors and , which are useful to construct projectors on the respective correction code bases [recall Eq. (15)], one can in principle find points in \({\varGamma }(\lambda _{11})\) determined by \(\lambda _{11} \in {\varOmega }\) and their phase angles \(\varphi _E\) and \(\varphi _F\). Following the proof of Theorem 1, one diagonalizes \(E_{11}\) and \(F_{11}\) using orthogonal matrices \(U_E\) and \(U_F\):

$$\begin{aligned} U_E = \begin{bmatrix} \cos \alpha _E&\sin \alpha _E \\ \sin \alpha _E&- \cos \alpha _E \end{bmatrix}, \qquad U_F = \begin{bmatrix} \cos \alpha _F&\sin \alpha _F \\ \sin \alpha _F&- \cos \alpha _F \end{bmatrix}. \end{aligned}$$
(22)

Now one introduces the transformed states: and , where and are given by:

(23)

where \(i = \sqrt{-1}\) denotes the complex imaginary unit. In order to determine these projection vectors, one has to find angles \(\theta _{\sigma }, \varphi _{\sigma }, \alpha _{\sigma }\) for \(\sigma \in \{E, F\}\) in terms of a point \(\widetilde{z} = \widetilde{x} + i \widetilde{y} \in {\varGamma }(\lambda )\) in the complex plane and a real parameter \(\lambda \). Let us suppose that we have found a point \(\widetilde{z} \in {\varGamma }(\lambda )\) for some value of \(\lambda \). We shall consider only the case for matrices \(E_{11}\) and \(E_{12}\) since for the matrices \(F_{11}\) and \(F_{12}\), the same pattern can be applied. Our current task is to find the projection vectors . Using Eq. (42), one can determine the azimuthal angle for projection \(\theta _E\) in terms of \(\lambda \):

$$\begin{aligned} \cos \theta _{E} = \frac{1}{\epsilon _E}\left( \lambda - \frac{1}{2} \hbox {Tr }E_{11}\right) , \end{aligned}$$
(24)

where \(\pm \epsilon _E\) denotes the eigenvalues of matrix \(E_{11}\) (defined in “Appendix”, in the text above Eq. 39). Using the definition of \(\widetilde{z} = \widetilde{z}(\varphi , \lambda ) \in {\varGamma }(\lambda )\) given in Eq. (45) with the following substitution: \((A, Z) \rightarrow (E_{12}, E_{11})\), one finds:

$$\begin{aligned} x(\varphi _E) = \frac{\widetilde{x} - x_0(\lambda )}{p(\lambda )}, \qquad y(\varphi _E) = \frac{\widetilde{y} - y_0(\lambda )}{p(\lambda )}. \end{aligned}$$

Using the above expressions and Eqs. (46) and (47), one can determine the polar angle \(\varphi _E\),

$$\begin{aligned} \cos \varphi _E = \frac{x r_2 - y r_1}{q_1 r_2 - r_1 q_2} = \frac{1}{p(\lambda )}\frac{r_2\left[ \widetilde{x} - x_0(\lambda )\right] - r_1\left[ \widetilde{y} - y_0(\lambda ))\right] }{q_1 r_2 - r_1 q_2}. \end{aligned}$$
(25)

The angle \(\alpha _E\) can be computed by recalling that \(\tan 2\alpha _E\) is related by Eq. (40) to elements of matrix \(E_{11}\). The family of states is then given by:

(26)

In a similar manner, we can determine the projection vectors for matrices \(F_{11}\) and \(F_{12}\). Thus, our initial problem of finding error correction code and solving the compression Eq. (3), equivalent to the set of Eq. (17) for projectors and , is reduced into a geometric problem of finding intersection points of two elliptic curves in the complex plane. The quantum code subspace—the projection \(P_2\) from Eq. (3)—is then given by a matrix of order four, given by the direct sum of two projectors of size two:

(27)

5 Exemplary quantum error correction codes

In this section, we present two examples of non-unitary quantum channels and determine their quantum error correction code subspaces \(\mathcal {H}^n_C\) using the method described in the previous section.

5.1 Simplified two-qubit amplitude-damping channel

The amplitude-damping channel (AD channel) is an important channel describing effects due to loss of energy of a quantum system [6]. Here, we consider two-level systems (qubits), but channels describing arbitrary n-level systems are also known [12]. Moreover, an interesting study of generalized amplitude-damping channels based on approximate quantum error correction schemes appeared recently [17]. Exemplary physical processes which can be described by this channel include the relaxation of atom from its excited state to the ground state [13], sending a quantum state from one location to another using a spin chain [14, 15] and attenuation of a photon in a cavity [16]. The Kraus representation of one-qubit amplitude-damping channel acting on state \(\rho \) with probability p, where \(0 \le p \le 1\) is given by:

$$\begin{aligned} {\varPhi }^{1AD}(\rho ) = B_1(p) \rho B_1^{\dagger }(p) + B_2(p) \rho B_2^{\dagger }(p), \end{aligned}$$

where the Kraus operators \(B_1(p)\) and \(B_2(p)\) are defined:

$$\begin{aligned} B_1(p) = \begin{bmatrix} 1&0 \\ 0&\sqrt{1-p} \end{bmatrix},\qquad B_2(p) = \begin{bmatrix} 0&\sqrt{p} \\ 0&0 \end{bmatrix}. \end{aligned}$$

To extend this channel into the two-qubit system, one can consider a product of two one-qubit channels:

$$\begin{aligned} {\varPhi }^{2AD} = {\varPhi }_1^{1AD} \otimes {\varPhi }_2^{1AD}. \end{aligned}$$

This channel is given by four Kraus operators. The bi-partite channel will be described by a sum of four new Kraus operators (all possible tensor products of \(B_1\) and \(B_2\)) acting on a two-qubit state \(\rho _1 \otimes \rho _2\). Assuming that the damping occurs with probability \(p_1\) on the first qubit and \(p_2\) on the second, we may write:

$$\begin{aligned} C_{ij} = B_i(p_1)\otimes B_j(p_2), \qquad 1 \le i,\quad j \le 2. \end{aligned}$$

To simplify the model, we consider a two-qubit channel defined by two Kraus operators: \(A_1 = \sqrt{\mathbb {I} - A_2^{\dagger } A_2}\) and \(A_2 = C_{12} = B_1 (p_1) \otimes B_2 (p_2)\), where \(A_1\) is defined up to a unitary transformation. We make the following choice of Kraus operators:

The trace-preserving channel analyzed here is then given by:

$$\begin{aligned} {\varPhi }^{2AD}(\rho _1 \otimes \rho _2) = A_{1}^{\dagger }(\rho _1 \otimes \rho _2)A_{1} + A_{2}^{\dagger }(\rho _1 \otimes \rho _2)A_{2}. \end{aligned}$$
(28)

One can determine the quantum error correction code for this channel by solving the compression problem given in Eq. (3). Let us solve it using the geometrical method presented in Sect. 4. In the first step, we compute matrices \(T_{11}\) and \(T_{12}\) given in Eq. (4), from which we obtain matrices \(E_{ij}\) and \(F_{ij}\) for \(1 \le i, j \le 2\):

Our aim is to find projection operator \(P_2 = P_E \oplus P_F\) from Eq. (14), where \(P_E\) and \(P_F\) are given by (15), which is equivalent to set of Eq. (17). To find these operators, we use notion of nuclear numerical range. We first compute the intersection \({\varGamma }(\lambda _{11})\), where \(\lambda _{11} \in {\varOmega }= W(E_{11})\cap W(F_{11})\), and from this intersection, we determine projection operators \(P_E\) and \(P_F\).

Fig. 1
figure 1

Standard numerical ranges \(W(E_{11})\) and \(W(F_{11})\) for the noise model given by (28) form line segments on the real axis. For clarity, the range \(W(E_{11})\) is lifted up. Their intersection forms the set \({\varOmega }\), which determines all possible values of the parameter \(\lambda _{11}\)

Let us introduce kets and , where diagonalizing unitary matrices \(U_E\) and \(U_F\) are given by Eq. (22). Since \(E_{11}\) and \(F_{11}\) are already diagonal, we may assume that \(U_E = U_F = \mathbb {I}\) and thus and , where and are yet arbitrary and parametrized according to Eq. (23). In order to find the azimuthal angles \(\theta _{E}\) and \(\theta _{F}\), one has to assure that the set \({\varOmega }\) is not empty, which means that the overlapping condition (19) holds. In this case, this condition reduces to the following two expressions:

$$\begin{aligned} p_1 \le 1 \quad \text {and} \quad p_2 \ge 0. \end{aligned}$$

Both of above conditions are satisfied by \(p_1\) and \(p_2\). We can now write respective standard numerical ranges as:

$$\begin{aligned} W(E_{11}) = 1 - \frac{p_2}{2} \left( 1 - \cos \theta _E\right) , \qquad W(F_{11}) =1 - \frac{p_2}{2}(1-p_1)(1-\cos \theta _F). \end{aligned}$$

This allows us to conclude that the set \({\varOmega }\) is given by the intersection of two sets, \({\varOmega }= W(E_{11})\cap W(F_{11}) = \bigl [ 1-(1-p_1)p_2, 1\bigr ]\), as shown in Fig. 1. If we now treat \(\theta _E\) and \(\theta _F\) as functions of the parameter \(\lambda _{11}\in {\varOmega }\), we get:

$$\begin{aligned} \cos \theta _E = 1 - \frac{2}{p_2}(1-\lambda _{11}) \quad \text {and} \quad \cos \theta _F = 1 - \frac{2}{p_2 (1-p_1)} \left( 1 -\lambda _{11}\right) . \end{aligned}$$
(29)

In order to find \(W(E_{12}|E_{11}-\lambda _{11}\mathbb {I})\) and \(W(F_{12}|F_{11}-\lambda _{11}\mathbb {I})\), we first compute the standard numerical range of matrices \(E_{12}\) and \(F_{12}\):

Using Eq. (29) in the above expressions, we obtain respective nuclear numerical ranges:

$$\begin{aligned} W(E_{12}|E_{11}-\lambda _{11} \mathbb {I})&= (1-\lambda _{11})\sqrt{\frac{1}{p_2} - 1}, \\ W(F_{12}|F_{11}-\lambda _{11} \mathbb {I})&= \frac{1}{\sqrt{p_2(1-p_1)}} e^{i \varphi _F} \left[ (1-\lambda _{11})(\lambda _{11} + p_2 - p_1 p_2 -1)\right] ^{1/2}. \end{aligned}$$

Here, \(\lambda \in {\varOmega }= [1 - (1-p_1)p_2, 1]\) can be treated as a free parameter. First of the above two sets is a line segment placed on the real axis while the second one forms a circle centered at zero. The non-trivial intersection of these two sets is possible only if \(\varphi _F = 0\) and:

$$\begin{aligned} \lambda _{11} = 1 - \frac{p_2(1-p_1)}{2-p_1-p_2+p_1 p_2}. \end{aligned}$$

Thus, the error correction code for this specific model is determined by the projector \(P_2\) of the form (27) with vectors and defined in Eq. (23) with the following horizontal angles:

$$\begin{aligned} \theta _E&= \cos ^{-1} \left( \frac{p_1 - p_2 + p_1 p_2}{2-p_1-p_2+p_1 p_2} \right) ,\\ \theta _F&= \cos ^{-1} \left( \frac{-p_1 - p_2 + p_1 p_2}{2-p_1-p_2+p_1 p_2} \right) . \end{aligned}$$

The polar angle \(\varphi _E\) can be chosen arbitrarily: \(\varphi _E \in [0, 2 \pi )\), as shown in Fig. 2.

Fig. 2
figure 2

Standard numerical ranges \(W(E_{12})\) (thick line segment on real axis) and \(W(F_{12})\) (dark disk) for \(p_1 = 0.5\) and \(p_2 = 0.7\). Intersection of their respective nuclear numerical ranges \({\varGamma }(\lambda _{11}) = W(E_{12}|E_{11}-\lambda _{11} \mathbb {I}) \cap W(F_{12}|F_{11}-\lambda _{11} \mathbb {I}) \) is given by a single point (dot) at approximately \({\varGamma }(0.696\ldots )\approx (0.2, 0)\)

5.2 General block-diagonal channel of length two

Let us now consider a more general case of a quantum channel of length two (\(k = 2\)) acting on system consisting of two qubits:

$$\begin{aligned} {\varPhi }^{G}(\rho ) = A_1 \rho A_1^{\dagger } + A_2 \rho A_2^{\dagger }. \end{aligned}$$
(30)

Our motivation is to find the most general noise model with a maximal number of free parameters, whose Kraus representation consists of two block-diagonal matrices \(\{ A_1, A_2 \}\). Operators \(A_i\) in general contain 16 free variables \(\{a_1,\, \ldots a_{16}\}\), where \(a_i \in \left( 0, 1\right) \) and \(1 \le i \le 16\). The condition that \({\varPhi }^{G}(\rho )\) preserves the trace:

$$\begin{aligned} A_1 A_1^{\dagger } + A_2 A_2^{\dagger } = \mathbb {I} \end{aligned}$$
(31)

imposes additional six constraints (there are eight equations for nonzero block-diagonal elements, from which two are not independent), so that we have 10 free parameters in total. Without loss of generality, we can choose any 6 parameters \(a_i\) from the set \(\{a_1,\, a_2, \ldots a_{16}\}\) and express them in terms of the remaining ones. We label them by \(b_j\) for \(1 \le j \le 6\) so that \(\{b_j\}\) are all dependent. Let us label the vector of free parameters by \(\mathbf {a} = (a_1,\, a_2, \ldots a_{10})\). Having this in mind, we can write Kraus operators \(A_i\) in the following form:

$$\begin{aligned} A_1 = \begin{bmatrix} a_1&a_2&0&0 \\ a_3&a_4&0&0 \\ 0&0&a_5&a_6 \\ 0&0&a_7&a_8 \\ \end{bmatrix}, \qquad A_2 = \begin{bmatrix} b_1&b_2&0&0 \\ a_9&b_3&0&0 \\ 0&0&b_4&b_5 \\ 0&0&a_{10}&b_6 \\ \end{bmatrix},&\end{aligned}$$
(32)

where the trace-preserving condition (31) implies:

$$\begin{aligned} b_1&= \sqrt{1-a_1^2-a_3^2-a_9^2}, \quad b_2 = \frac{\left( a_1 a_2+a_3 a_4\right) b_1 - c_1a_9}{a_1^2+a_3^2-1},&\nonumber \\ b_3&= \frac{\left( a_1 a_2+a_3 a_4\right) a_9+ c_1 b_1}{a_1^2+a_3^2-1}, \quad b_4 = \sqrt{1-a_5^2-a_7^2-a_{10}^2},&\nonumber \\ b_5&= \frac{\left( a_5 a_6+a_7 a_8\right) b_4-c_2a_{10}}{a_5^2+a_7^2-1}, \quad b_6 =\frac{\left( a_5 a_6+a_7 a_8\right) a_{10}+ c_2 b_4}{a_5^2+a_7^2-1}.&\end{aligned}$$
(33)

To simplify notation, we introduce variables \(c_i\) which are functions of the independent parameters:

$$\begin{aligned} c_1&= \sqrt{\left( a_4^2-1\right) a_1^2-2 a_1 a_2 a_3 a_4-a_4^2+\left( a_2^2-1\right) \left( a_3^2-1\right) }&\nonumber \\ c_2&= \sqrt{\left( a_8^2-1\right) a_5^2-2 a_5 a_6 a_7 a_8-a_8^2+\left( a_6^2-1\right) \left( a_7^2-1\right) }&\end{aligned}$$
(34)

In order to find QECC for the map \({\varPhi }^{G}\), we proceed with the method described in Sect. 4. Let us begin by computing matrices \(T_{11}\) and \(T_{12}\). Defining:

figure a

one can write matrices \(T_{11}\) and \(T_{12}\) as:

$$\begin{aligned} T_{11}&= \begin{bmatrix} e_1&e_2&0&0 \\ e_2&e_3&0&0 \\ 0&0&f_1&f_2 \\ 0&0&f_2&f_3 \\ \end{bmatrix}, \qquad T_{12} = \begin{bmatrix} e_4&e_5&0&0 \\ e_6&e_7&0&0 \\ 0&0&f_4&f_5 \\ 0&0&f_5&f_6 \\ \end{bmatrix}. \end{aligned}$$
(35)

Matrices \(E_{11}\), \(F_{11}\), \(E_{12}\) and \(F_{12}\) are then:

$$\begin{aligned} E_{11}&= \begin{bmatrix} e_1&e_2 \\ e_2&e_3 \end{bmatrix}, \quad F_{11} = \begin{bmatrix} f_1&f_2 \\ f_2&f_3 \end{bmatrix}, \end{aligned}$$
(36)
$$\begin{aligned} E_{12}&= \begin{bmatrix} e_4&e_5 \\ e_6&e_7 \end{bmatrix}, \quad F_{12} = \begin{bmatrix} f_4&f_5 \\ f_6&f_7 \end{bmatrix}. \end{aligned}$$
(37)

Once again our aim is to find projection operator \(P_2\) satisfying Eq. (3). We will do it by first computing the intersection of nuclear numerical ranges as explained in Sect. 4. To do so, we first compute the intersection \({\varOmega }= W(E_{11})\cap W(F_{11})\) which, according to Eq. (19), is completely determined by the eigenvalues of matrices \(E_{11}\) and \(F_{11}\), denoted by \(\nu _{1}\), \(\nu _{2}\) and \(\mu _{1}\), \(\mu _{2}\), respectively. Without loss of generality, we may assume \(\nu _{2} > \nu _{1}\) and \(\mu _{2} > \mu _{1}\). The set \({\varOmega }\) is then given by:

$$\begin{aligned} {\varOmega }= [\max (\nu _1,\mu _1), \min (\nu _2, \mu _2)]. \end{aligned}$$
(38)

This set is not empty if condition (19) holds for matrices \(E_{11}\) and \(F_{11}\). Let us denote by \(\lambda _{11}\) a free parameter contained in \({\varOmega }\). In order to find an appropriate QECC, one can determine the set \({\varGamma }(\lambda _{11})\) defined in Eq. (21). We plot the set \({\varGamma }(\lambda _{11})\) for some convenient choice of parameters \(\{ a_i\}\) in Fig. 3.

Fig. 3
figure 3

Standard numerical ranges \(W(E_{12})\) and \(W(F_{12})\) (denoted by dark elliptic disks) for the choice of parameter vector \(\mathbf {a} = (a_1,\, a_2, \ldots a_{10}) = (0.9,\ 0.7,\ 0.2,\ 0.9,\ 0.6,\ 0.7,\ 0.9,\ 0.1,\ 0.6,\ 0.5)\). Nuclear numerical ranges in this example are given by two elliptic curves in bold parametrized by \(\lambda _{11}\). Their intersection \(W(E_{12}|E_{11}-\lambda _{11}\mathbb {I})\cap W(F_{12}|F_{11}-\lambda _{11}\mathbb {I})\) for a given value of \(\lambda _{11}\) forms the set \({\varGamma }(\lambda _{11})\), which in this case consists of two points (dots)

Following the reasoning from Sect. 4 and methodology present in the proof of Theorem 1, we conclude that for a given value of \(\lambda _{11}\) one can construct elliptic curves in the complex plane and determine their intersection points \(\widetilde{z} = \widetilde{x} + i \widetilde{y} \in {\varGamma }(\lambda _{11})\) using Theorem 1. Having obtained the set \({\varGamma }(\lambda )\) and using Eqs. (24) and (25) one can then determine the azimuthal and polar angles \((\theta _E,\ \varphi _E)\) for , and analogous angles \((\theta _F,\ \varphi _F)\) for vector , respectively. The projection vectors and , which form the projection operator \(P_2\), can be computed using Eq. (26). The correction code subspace for this particular noise model is then given by Eq. (27).

6 Concluding remarks

In this work, we have introduced the notion of nuclear numerical range W(A|Z) of an operator A with respect to an auxiliary operator Z and demonstrated that it allows one to find quantum error correction codes protecting against noise. In particular, this technique works for models of quantum errors with non-unitary noise operators. Using a simple geometric construction involving an intersection of two ellipses in the complex plane, we found such a quantum error correction code for a simplified model of two-qubit amplitude-damping channel and a general noise model with two Kraus operators of size \(n = 4\) with block-diagonal structure. Note that the method used here for the two-qubit system is straightforward to generalize for larger dimensions. We expect that further development of this technique will allow for effective construction of quantum error correction codes protecting information against more general non-unitary noise models.