Perturbation in Linear Response Theory (classical formalism) A N-particles system is described by the following Hamiltonian
$$H=H_0+H'(t)$$
where $H_0$ is the unperturbated Hamiltonian and $H'(t)$ is the perturbation
$$H'(f)=-A\cdot\mathcal{F}(t)$$ 
written as the coupling of a "characteristic lenght" related to the system ($A$) and a "generalized external force" ($\mathcal{F}(t)$): infact, the product of the two has to be energy-dimensional.
The phase-space distribution function of this N-particles system is written as follows:
$$\rho(t)=\rho_0(t)+\Delta\rho(t)$$
where the first contribution is the one related to the unperturbated system and the second one is the one due to the perturbation (first-order only, being interested in LINEAR response theory!).
Now, defining the Liouville Operator as
$$\mathcal{L}\equiv i\{H,\cdot\}$$
it is possible to to write the equation of motion of $\rho(t)$ as
$$\frac{\partial\rho(t)}{\partial t}=-i\mathcal{L_0}\rho(t)-\{A,\rho(t)\}\mathcal{F}(t)$$
($\mathcal{L_0}\equiv i\{H_0,\cdot\}$) and - substituing the expression $\rho(t)=\rho_0(t)+\Delta\rho(t)$ in it, simplifying and neglecting higher-order terms - one can derive the following equation:
$$\frac{\partial\Delta\rho(t)}{\partial t}=-i\mathcal{L_0}\Delta\rho(t)-\{A,\rho_0(t)\}\mathcal{F}(t)$$
My book says that integrating the last differential equation (having set the origin at $t=-\infty$) one can obtain the solution
$$\Delta\rho(t)=-\int_{-\infty}^{t} e^{-i(t-s)\mathcal{L_0}} \{A,\rho_0(s)\}\mathcal{F}(s) ds$$
I think I understand the physical meaning of the last equation but here comes my question: how can I derive it myself? I don't seem to be able to do the math, because I don't know how to work with operators (the Liouville one) in this type of problems.
Hope someone will help me, thank you in advance for your time!
 A: We'll give some details based on Andrew's comments.
We are trying to solve
$$
\frac{ \partial \Delta \rho(t) }{ \partial t }
=
-i\mathcal L_0 \Delta \rho(t)
-\{A, \rho_0(t) \} \, \mathcal F(t).
\qquad (1)
$$
Homogeneous solution
To develop some familiarity with the Liouville operator, let us first consider the special case of $A = 0$.
Then we have a homogeneous equation:
$$
\frac{ \partial \Delta \rho(t) }{ \partial t }
=
-i\mathcal L_0 \Delta \rho(t)
\qquad (2)
$$
Note that the Liouville operator $\mathcal L_0$ does not explicitly depend on $t$
because the unperturbed Hamiltonian $H_0(q, p)$ is not time-dependent.
That is, for any function $f(q, p)$, we have
$$
\begin{align}
-i\mathcal L_0 f(q, p)
&= \{H_0(q, p), f(q, p)\} \\
&= \sum_i \frac{ \partial H_0} { \partial q_i } \frac{ \partial f} { \partial p_i }
- \frac{ \partial H_0} { \partial p_i } \frac{ \partial f} { \partial q_i },
\end{align}
$$
which is a time-independent expression.
In short, $\mathcal L_0$, albeit its operator nature, is a constant of time.
This means that we can integrate $(2)$ as if it were a regular ODE:
$$
\begin{aligned}
\Delta \rho(t)
=
e^{-i t\mathcal L_0} g(q, p),
\end{aligned}
\qquad (3)
$$
where $g(q, p)$ is some time-independent function.  To be safe, we can plug this into the $(2)$ to directly verify that this is indeed a solution:
$$
\begin{aligned}
\mbox{l.h.s.}=
\frac{\partial \Delta \rho(t) }{\partial t}
=
i \mathcal L_0 \,e^{-i t\mathcal L_0} g(q, p)
=
i \mathcal L_0 \, \Delta \rho(t) = \mbox{r.h.s.}.
\end{aligned}
$$
Special solution
We can similar borrow the technique for regular ODEs to solve $(1)$.  Denote
$$
B(t) \equiv -\{A, \rho_0(t)\} F(t),
$$
then $(1)$ is equivalent to
$$
\frac{ \partial \Delta \rho(t) }{ \partial t }
=
-i\mathcal L_0 \Delta \rho(t)
+B(t).
\qquad (4)
$$
To solve $(4)$, we shall now treat $g(q, p)$ in the homogeneous solution to be time dependent, $g(q, p, t)$, i.e., we postulate a trial solution
$$
\begin{aligned}
\Delta \rho(t)
=
e^{-i t\mathcal L_0} g(q, p, t).
\end{aligned}
\qquad (5)
$$
Then
$$
-i\mathcal L_0 e^{-i t\mathcal L_0} g(q, p, t)
+
e^{-i t\mathcal L_0} \frac{\partial}{\partial t} g(q, p, t)
=
-i\mathcal L_0 e^{-i t\mathcal L_0} g(q, p, t)
+B(t),
$$
or
$$
\frac{\partial}{\partial t} g(q, p, t)
=
e^{i t\mathcal L_0} B(t).
$$
The right hand side is just a function, so
$$
g(q, p, t)
=
\int_{-\infty}^t e^{i s\mathcal L_0} B(s) \, ds,
$$
assuming $B(-\infty) = 0$.
Using this in $(5)$, we get
$$
\begin{aligned}
\Delta \rho(t)
&=
e^{-i t\mathcal L_0} \int_{-\infty}^t e^{i s\mathcal L_0} B(s) \, ds \\
&=
\int_{-\infty}^t e^{-i (t - s) \mathcal L_0} B(s) \, ds \\
&=
-\int_{-\infty}^t e^{-i (t - s) \mathcal L_0} \{A, \rho_0(s)\} \mathcal F(s) \, ds.
\end{aligned}
$$
This is the desired result.  To ensure nothing went wrong, we can always explicitly plug this solution in $(1)$.
In summary, our solution is a simple copy of the general solution of a regular ODE.  The key is to remember that $\mathcal L_0(q, p)$, unlike $\mathcal L(q, p, t)$, does not explicitly depend on time.
