Suppose we have a depolarizing channel operation $$E(\rho)=\frac{p}{2}\textbf{1}+(1-p)\rho$$ acting on a Spin$\frac{1}{2}$ density matrix of the form $\rho=\frac{1}{2}(\textbf{1}+\textbf{s}\cdot\textbf{$\sigma$})$. I have found the Kraus operators to be: $$E_1=\sqrt{\left(1-\frac{3}{4}p\right)}\textbf{1}, E_2=\frac{\sqrt{p}}{2}\sigma_x,E_3=\frac{\sqrt{p}}{2}\sigma_y \text{ and } E_4=\frac{\sqrt{p}}{2}\sigma_z$$ I am now supposed to find the unitary matrix U such that the Operation can be expressed in a bigger system i.e. after adding a System S. As far as I understand it, the new operation can be written as: $$E(\rho)=\sum_kE_k\rho E_k^\dagger=\text{Tr}_S(U\rho\otimes\rho_EU^\dagger)$$ Supposing the new system S is prepared on a state $|e_0\rangle$, How do I find the correct unitary matrix?

I appreciate your cooperation.

Crossposted to qc: https://quantumcomputing.stackexchange.com/questions/13548/how-to-find-the-unitary-operation-of-a-depolarizing-channel