No need to iterate 
 \begin{equation}
    H_I(t)=H_I(0)+\frac{i}{\hbar}\int\limits_0^t dt' \, [H_0,H_I(t')] .
\end{equation}
Just observe that it is readily solved by the Ansatz
\begin{equation}
    H_I(t)=e^{\frac{i}{\hbar}tH_0 }H_I(0)e^{-\frac{i}{\hbar} t H_0  }
\end{equation}
which yields the r.h.side
$$
H_I(0)+\frac{i}{\hbar}\int\limits_0^t dt'  ~~ [H_0,e^{\frac{i}{\hbar}t' H_0  }H_I(0)e^{-\frac{i}{\hbar}t' H_0  }  ]\\ 
= H_I(0)+\frac{i}{\hbar}\int\limits_0^t dt'  ~~e^{\frac{i}{\hbar}t'[H_0  }[H_0,H_I(0)] =  H_I(0)+ \left ( e^{\frac{i}{\hbar}t[H_0  } H_I (0) - H_I(0)  \right )\\  =  e^{\frac{i}{\hbar}t[H_0  } H_I (0)   \equiv  e^{\frac{i}{\hbar}t H_0 } H_I (0)  e^{-\frac{i}{\hbar}t H_I }.
$$
This last line is a famous [Lemma (4)](https://en.wikipedia.org/wiki/Derivative_of_the_exponential_map), namely that 
$$
e^A B e^{-A}= B+ [A,B]+ [A,[A,B]]/2!+ [A,[A,[A,B]]]/3!+...
$$
where I have used the notation 
$$
e^{[A} \equiv e^{\operatorname{ad}_A } = \operatorname{Ad}_{e^A}
$$
since most physicists are unfamiliar with it. The stunt then is integration of a simple exponential of the ad operator.

The very same integration of ad maneuver will net you the density matrix.