Deriving Nakajima–Zwanzig

Question

I would like help on one step of the derivation of the Nakajima–Zwanzig equations. Specifically given the below equations: $${\partial_t}\left( \begin{matrix} \mathcal{P} \\ \mathcal{Q} \\ \end{matrix} \right)\rho =\left( \begin{matrix} \mathcal{P} \\ \mathcal{Q} \\ \end{matrix} \right)L\left( \begin{matrix} \mathcal{P} \\ \mathcal{Q} \\ \end{matrix} \right)\rho +\left( \begin{matrix} \mathcal{P} \\ \mathcal{Q} \\ \end{matrix} \right)L\left( \begin{matrix} \mathcal{Q} \\ \mathcal{P} \\ \end{matrix} \right)\rho,$$

I would like to understand how to formally solve them to get the below equation: $$\mathcal{Q}\rho ={{e}^{\mathcal{Q}Lt}}Q\rho (t=0)+\int_{0}^{t}dt'{e}^{\mathcal{Q}Lt'}\mathcal{Q}L\mathcal{P}\rho (t-{t}') $$

Many places just say the equations can be formally solved, such as:

1. This SE question here, which gives me the impression that I could check this by substituting in the solution. However I would like to derive the formal solution (i.e., if I was not given the solution, how would one go about deriving it).
2. This reference simplifies the problem down to the core constituents, which is a coupled first order set of differential equations. However it then introduces Green functions, which I am not familiar with in the context of solving coupled first order equations.
3. The Wikipedia reference here includes a bit more information in their note 3 that includes an exponential, but I don't know where that exponential comes from.