Show that $I = frac{rho + sigma_xrhosigma_x +sigma_yrhosigma_y + sigma_zrhosigma_z}{2}$ for all states $rho$

For arbitrary one qubit density matrix we have:

$$rho = frac{I}{2} + frac{r_x sigma_x + r_y sigma_y + r_z sigma_z}{2}$$

where $|r| le 1$. Here we should take into account that $sigma_i sigma_j sigma_i = -sigma_j$ where $i ne j$ and $i, j in {x, y, z}$, and, also $sigma_isigma_i=I$. With this we will obtain the equality presented in the question. Let's see, for example what will be equal the $sigma_x rho sigma_x$ term:

$$sigma_x rho sigma_x = frac{I}{2} + frac{r_x sigma_x - r_y sigma_y - r_z sigma_z}{2}$$

Similarly, we will obtain:

$$frac{rho + sigma_xrhosigma_x +sigma_yrhosigma_y + sigma_zrhosigma_z}{2} = =I + frac{2r_x sigma_x + 2r_y sigma_y + 2r_z sigma_z -2r_x sigma_x - 2r_y sigma_y - 2r_z sigma_z}{2} = I $$

Given that you're only working with a single qubit state it is possible to also show this by direct calculation on a parametrised state. That is, we can write any single qubit $rho$ as $$ rho = begin{pmatrix} a & beta overline{beta} & 1-a end{pmatrix} $$ with $ain[0,1]$ and $beta in mathbb{C}$ such that $(1-2a)^2 + 4 |beta|^2 leq 1$. Then we can directly compute the action of the Pauli conjugation $$ sigma_x rho sigma_x = begin{pmatrix} 1-a & overline{beta} beta & a end{pmatrix} $$ $$ sigma_y rho sigma_y = begin{pmatrix} 1-a & -overline{beta} -beta & a end{pmatrix} $$ $$ sigma_z rho sigma_z = begin{pmatrix} a & -beta -overline{beta} & 1-a end{pmatrix}. $$ Summing these up with $rho$ and dividing through by $2$ we get the desired result.

Assuming $rho$ is a pure state, these are the explicit calculations. You can easily generalize to mixed states.

$newcommand{ket}[1]{|{#1}rangle}$ $newcommand{bra}[1]{langle{#1}|}$ Let $ket{psi} = alphaket{0}+ betaket{1}$, where $alpha,beta in mathbb{C}^2$ and $|alpha|^2 + |beta|^2 = 1$.

Thus $rho = ket{psi}bra{psi} = begin{pmatrix} alphaalpha^* & alphabeta^* betaalpha^* & betabeta^*end{pmatrix}$.

Now $sigma_x = begin{pmatrix} 0 & 1 1 & 0end{pmatrix}$, so $ sigma_x rho sigma_x = begin{pmatrix} betabeta^* & betaalpha^* alphabeta^* & alphaalpha^*end{pmatrix}$.

Similarly, you can compute $sigma_yrhosigma_y = begin{pmatrix} betabeta^* & -betaalpha^* -alphabeta^* & alphaalpha^*end{pmatrix}$ and $sigma_zrhosigma_z = begin{pmatrix} alphaalpha* & -alphabeta* -betaalpha* & betabeta*end{pmatrix}$.

Finally, summing up:

$frac{1}{2} (rho + sigma_x rho sigma_x + sigma_y rho sigma_y + sigma_z rho sigma_z) = frac{1}{2}begin{pmatrix} 2(alphaalpha^* + betabeta^*) & 0 0 & 2(alphaalpha^* + betabeta^*) end{pmatrix} = I$.

This method is largely similar to Davit's (this covers a slightly more general case where $rho$ is any arbitrary matrix with trace 1, and you easily see how to adjust it without the trace 1 condition). Any $2times 2$ matrix can be decomposed as $aI+vec{n}cdotvec{sigma}$ if we allow $a$ and $vec{n}$ to take on arbitrary complex values. Moreover, two $2times 2$ matrices are equal if an only if their values of $a$ and $vec{n}$ are equal. So, let $$ tau=frac{rho + sigma_xrhosigma_x +sigma_yrhosigma_y + sigma_zrhosigma_z}{2}. $$ We want to show that $a=1$ and $vec{n}=0$. Now, $$ a=text{Tr}(tau)/2,qquad n_i=text{Tr}(sigma_itau)/2. $$ Remember that trace is invariant under cyclic permutations, so $$ a=frac{1}{4}text{Tr}(rho + sigma_xrhosigma_x +sigma_yrhosigma_y + sigma_zrhosigma_z)=frac{1}{4}text{Tr}(rho + rhosigma_x^2 +rhosigma_y^2 + rhosigma_z^2)=text{Tr}(rho)=1. $$ Similarly, $$ n_x=frac12text{Tr}(sigma_xrho + rhosigma_x +sigma_xsigma_yrhosigma_y + sigma_xsigma_zrhosigma_z)=frac12text{Tr}(2sigma_xrho +rhosigma_ysigma_xsigma_y + rhosigma_zsigma_xsigma_z). $$ Now use the anti-commutation properties of the Pauli matrices to get $$ n_x=frac12text{Tr}(2rhosigma_x -rhosigma_x - rhosigma_x)=0. $$ The other two components are just the same.

Chapter VII. E. in Daniel Lidar's notes. Use $rho = frac{1}{2}(I + vec{v}cdotvec{sigma})$ and products of Pauli matrices:

Check for each pair that: $sigma_i sigma_j = delta_{ij} I + i epsilon_{ijk}sigma_k$

Use it to show: $ sigma_i sigma_j sigma_k = delta_{ij} sigma_k - delta_{ik} sigma_j + delta_{jk} sigma_i + i epsilon_{ijk} I $

One more step $ sigma_i sigma_j sigma_i = 2delta_{ij} sigma_i - sigma_j = begin{cases} +sigma_j &, i = j -sigma_j &, i neq j end{cases} $

with this go to eq. 189 from Daniel Lidar: $$ sigma_x(I + vec{v}cdot vec{sigma}) sigma_x = I + v_x sigma_x - v_y sigma_y - v_z sigma_z $$ $$ sigma_y(I + vec{v}cdot vec{sigma}) sigma_y = I - v_x sigma_x + v_y sigma_y - v_z sigma_z $$ $$ sigma_z(I + vec{v}cdot vec{sigma}) sigma_z = I - v_x sigma_x - v_y sigma_y + v_z sigma_z $$

add it together with $$ I(I + vec{v}cdot vec{sigma}) I = I + v_x sigma_x + v_y sigma_y + v_z sigma_z $$

to get $$ 2(rho + sigma_x rho sigma_x + sigma_y rho sigma_y + sigma_z rho sigma_z) = 4I $$

This is a special instance of a general linear algebra result.

Note that the identity matrix $newcommand{vec}{operatorname{vec}}I$ can be decomposed as $I=sum_k v_kotimes v_k^*$ for any orthonormal basis ${v_k}_k$, and vice versa any such decomposition identifies the identity matrix.

Now notice that the Pauli matrices are an orthonormal basis in an enlarged Hilbert space, meaning that $$operatorname{Tr}[(sigma_i/sqrt2)(sigma_j/sqrt2)]=delta_{ij}.$$ More explicitly, this is saying that we can think of the matrices $sigma_i$ as orthonormal vectors in some space, $i.e.$ we have $langle vec(sigma_i/sqrt2),vec(sigma_j/sqrt2)rangle=delta_{ij}$, where $vec(B)$ is the vectorisation of the operator $B$.

---

If $A_ainmathrm{Lin}(mathcal X,mathcal Y)$ are a set of such orthonormal operators, we have $$mathrm{tr}(A_a^dagger A_b)=delta_{ab} Longleftrightarrow langlemathrm{vec}(A_a),mathrm{vec}(A_b)rangle=delta_{ab},$$ where $vec(A_a)inmathcal Yotimesmathcal X$ is the vectorisation of $A_a$. If the set is a basis, then we also have $$sum_a (A_a)_{12} (A_a^*)_{34} = delta_{13}delta_{24} Longleftrightarrow sum_avec(A_a)vec(A_a)^dagger = I_{mathcal Yotimesmathcal X} $$ Now, the statement we are interested in is of the form $sum_a A_a rho A_a^dagger = I$. This amounts to $$ sum_{a34} (A_a)_{13} (A_a^*)_{24} rho_{34} = delta_{12} Longleftrightarrow sum_a (A_aotimes A_a^*)vec(rho) = lvert mrangle, $$ where $lvert mrangleequivsum_k lvert k,krangle$. The question is thus, what type of operator is $sum_a A_aotimes A_a^*$? Componentwise, the relation with $sum_a vec(A_a)vec(A_a)^dagger$ is clear: $$(A_aotimes A_a^*)_{ij,nm} = (A_a)_{in} (A_a^*)_{jm} = (vec(A_a)vec(A_a)^dagger)_{in,jm} = delta_{ij}delta_{nm},$$ that is, $A_aotimes A_a^*$ is the Choi of $vec(A_a)vec(A_a)^dagger$. Summing over $a$ this is the identity, which means that $sum_a A_aotimes A_a^*$ is the Choi of the identity, which is the projector over the maximally entangled state: $$sum_a A_aotimes A_a^*=lvert mrangle!langle mrvert.$$ We conclude that $$sum_A A_a rho A_a^dagger = operatorname{unvec}left(sum_a (A_aotimes A_a^*) vec(rho)right) = operatorname{unvec}(lvert mrangle ) = I_{mathcal X}.$$

I know this is an old question, but I feel like giving, imo, the simplest answer following the question specified in N&C.

Firstly as specified in the question define:

$ mathcal{E}(A) = frac{A + XAX + YAY + ZAZ}{4}$ .

It it is easy to see that

$ mathcal{E}(I) = frac{I + XIX + YIY + ZIZ}{4} = frac{I + XX + YY + ZZ}{4} = I$

For the other three quantities $mathcal{E}(X),mathcal{E}(Y),mathcal{E}(Z)$ we can make use of the elementary identities:

$sigma_isigma_isigma_i = sigma_i$

$sigma_isigma_jsigma_i = -sigma_j$

Plugging in these identities we can see $mathcal{E}(X),mathcal{E}(Y),mathcal{E}(Z) = 0$. More generally, (after grouping the terms), we can see that:

$ mathcal{E}(sigma_i) = frac{2sigma_i - 2sigma_i}{4} = 0$

Finally we know from e.q 2.175 that

$rho = frac{I + vec{r} cdotvec{sigma}}{2} = frac{I + r_xsigma_x + r_ysigma_y + r_zsigma_z}{2}$,

and plugging this into

$mathcal{E}(rho) = frac{mathcal{E}(I) + mathcal{E}(r_xsigma_x) + mathcal{E}(r_ysigma_y) + mathcal{E}(r_zsigma_z)}{2}$,

using the results from above we see all the $mathcal{E}(r_isigma_i)=0$, leaving us just with

$mathcal{E}(rho) = frac{mathcal{E}left({I}right)}{2} = frac{I}{2}$,

finally to complete the proof

$2mathcal{E}(rho) = I$