Solution (source code)

= Solution

Write $\rho_0=|\phi\rangle\langle\phi|=(I+r_xX+r_yY+r_zZ)/2$. Conjugating by a <Pauli matrix> keeps its own Bloch component and reverses the other two; for example,
$$
X\rho_0X=\tfrac12(I+r_xX-r_yY-r_zZ).
$$
In the average of the four conjugations, the coefficient of every nonidentity <Pauli matrix> is zero. Consequently
$$
\boxed{\frac{\rho_0+X\rho_0X+Y\rho_0Y+Z\rho_0Z}{4}=\frac I2.}
$$
The <Pauli operator basis> gives the Bloch representation of a <qubit> <density operator>.<Foundations I: States and Ensembles by John Preskill>[§2.3.2, printed pp. 18–19, equations (2.62)–(2.65)]{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.