There are (at least) two ways to perturbatively solve a matrix initial value problem: the Dyson expansion and the Magnus expansion
To be explicit, suppose you're solving for a (interaction-picture) density matrix $\rho(t)$ at time $t$ characterized by a Hamiltonian $H = H_0 + H_{\mathrm{interaction}}$, where your initial state is $\rho(0)\equiv \rho_0$. The Dyson expansion gives a perturbative solution: $$ \rho(t) = \rho_0 - i \int_0^t ds\ \left[ H_{\mathrm{I}}(s), \rho_0 \right] + (-i)^2 \int_0^t ds \int_0^s dr\ \left[ H_{\mathrm{I}}(s), \left[ H_{\mathrm{I}}(r), \rho_0 \right] \right] + \ldots $$
Where $H_{\mathrm{I}}(t)=e^{+iH_0 t} H_{\mathrm{interaction}} e^{-iH_0 t}$ in the interaction picture.
I have heard that the Magnus expansion is useful in that it preserves the symplectic form of the time evolution at every order. I heard that this means that at every order, the volume of phase space is preserved.
What does this mean precisely?
(In contrast, the utility of the Dyson series is that probabilities are preserved at every order, meaning the trace of each contribution order-by-order remains $1$, this I can understand)
