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 the mathematics operators involved. 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.