In my reference, Page 375, Chapter 8, Quantum Computation and Quantum Information by Nielsen and Chuang, it is given that
The Pauli matrices, along with the identity matrix $I$, form an orthonormal basis for the Hilbert space of all $2\times 2$ complex matrices, $\mathbb{M}_{2\times 2}(\mathbb{C})$.
$\therefore$ The operators $E_i$ can be written in the form, $$ E_i=\alpha_i I+\sum_{k=1}^3a_{ik}\sigma_k $$ $$ \rho=\frac{I+\vec{r}.\vec{\sigma}}{2}=\frac{1}{2}\begin{bmatrix}1+r_z&r_x-ir_y\\r_x+ir_y&1-r_z\end{bmatrix} $$ where $\vec{r}=(r_x,r_y,r_z)$ and $\vec{\sigma}=(\sigma_1,\sigma_2,\sigma_3)=\Big(\begin{bmatrix}0&1\\1&0\end{bmatrix},\begin{bmatrix}0&-i\\i&0\end{bmatrix},\begin{bmatrix}1&0\\0&-1\end{bmatrix}\Big)$
The Pauli matrices satisfy the following properties, $$ \sigma_j\sigma_k=\delta_{jk}I+i\epsilon_{jkl}\sigma_l\\ \sigma_1^2=\sigma_2^2=\sigma_3^2=-i\sigma_1\sigma_2\sigma_3=I $$
\begin{align} \mathcal{E}&(\rho)=\sum_l E_l\rho E_l^\dagger=\frac{I+\sum_l E_l(\vec{r}.\vec{\sigma})E_l^\dagger}{2}\\ &=\frac{1}{2}\bigg[I+\sum_l\Big[ \Big(\alpha_lI+\sum_{j=1}^3a_{lj}\sigma_j\Big)\Big(\sum_{k=1}^3r_k\sigma_k\Big)\Big(\alpha_l^*I+\sum_{p=1}^3a^*_{lp}\sigma_p\Big) \Big]\bigg]\\ &=\frac{1}{2}\bigg[I+\sum_l\Big[ |\alpha_l|^2\sum_{k=1}^3r_k\sigma_k+\alpha_l\sum_{k=1}^3r_k\sigma_k\sum_{p=1}^3a_{lp}^*\sigma_p+\alpha_l^*\sum_{j=1}^3a_{lj}\sigma_j \sum_{k=1}^3r_k\sigma_k+\sum_{j=1}^3a_{lj}\sigma_j\sum_{k=1}^3r_k\sigma_k\sum_{p=1}^3a_{lp}^*\sigma_p \Big]\bigg]\\ &=\frac{1}{2}\bigg[I+\sum_l\Big[ |\alpha_l|^2\sum_{k=1}^3r_k\sigma_k+\alpha_l\sum_{k=1}^3\sum_{p=1}^3r_ka_{lp}^*\sigma_k\sigma_p+\alpha_l^*\sum_{j=1}^3\sum_{k=1}^3r_ka_{lj}\sigma_j \sigma_k+\sum_{j=1}^3\sum_{k=1}^3\sum_{p=1}^3r_ka_{lj}a_{lp}^*\sigma_j\sigma_k\sigma_p \Big]\bigg]\\ \end{align}
I think I need to rearrange this in the form $\mathcal{E}(\rho)=\frac{I+(M\vec{r}+\vec{c}).\vec{\sigma}}{2}$ to obtain $M_{jk}=\sum_{l}\bigg[ a_{lj}a_{lk}^*+a_{lj}^*a_{lk}+(|a_l|^2-\sum_p a_{lp}a_{lp}^*)+i\sum_p\epsilon_{jkp}(\alpha_la_{lp}^*-\alpha_l^*a_{lp}) \bigg]$ and $c_{_k}=2i\sum_l\sum_{jp}\epsilon_{jkp}a_{lj}a_{lp}^*$. How do I proceed further and obtain the required expression for $M_{jk}$ and $c_k$ ?
This is asked in Affine Map of the Bloch sphere but it does not contain any solution.
This is also given as in Exercise 6.1, Page 338, Principles of Quantum Computation II by Giuliano Benenti and Giulio Casat