 Solution

Write ${\rho }_0=|\varphi ⟩⟨\varphi |=(I+r_{x}X+r_{y}Y+r_{z}Z)/2$. Conjugating by a Pauli matrix keeps its own Bloch component and reverses the other two; for example,In the average of the four conjugations, the coefficient of every nonidentity Pauli matrix is zero. ConsequentlyThe Pauli operator basis gives the Bloch representation of a qubit density operator.^[ref]

This is the qubit Pauli twirl, the dimension-two case of the Heisenberg-Weyl twirling channel. The argument uses only trace one, so it also works for mixed qubit inputs. The preparation label must be forgotten: conditioning on the label would retain a known pure branch rather than the maximally mixed state obtained by the ensemble average.