# 12.5: Perturbative solution of the Liouville equation


Substituting the perturbative form for $$f (x, t )$$ into the Liouville equation, one obtains

${\partial \over \partial t}(f_0({\rm x}) + \Delta f({\rm x}, t ) ) + (iL_0 + i\Delta L(t))(f_0({\rm x}) + \Delta f({\rm x},t)) = 0 \nonumber$

Recall $$\partial f_0/\partial t=0$$. Thus, working to linear order in small quantities, one obtains the following equation for $$\Delta f({\rm x},t)$$:

$\left({\partial \over \partial t} + iL_0\right)\Delta f({\rm x},t) =-i\Delta L f_0({\rm x}) \nonumber$

which is just a first-order inhomogeneous differential equation. This can easily be solved using an integrating factor, and one obtains the result

$\Delta f({\rm x},t) = -\int_0^t ds e^{-iL_0(t-s)}i\Delta L(s)f_0({\rm x}) \nonumber$

Note that

$i\Delta L f_0({\rm x}) = iL f_0({\rm x}) - iL_0 f_0({\rm x}) = iL f_0({\rm x})=\dot{\rm x}\cdot\nabla_{\rm x}f_0({\rm x}) \nonumber$

But, using the chain rule, we have

\begin{align*} \dot{\rm x}\cdot\nabla_{\rm x}f_0({\rm x}) &= \dot{\rm x}\cdot{\partial f_0 \over \partial H}{\partial H \over \partial {\rm x}} \\[4pt] &= {\partial f_0 \over \partial H}\sum_{i=1}^{3N}\left[\dot{p}_i{\partial H \over \partial p_i} +\dot{q}_i {\partial H \over \partial q_i}\right] \\[4pt] &={\partial f_0 \over \partial H}\sum_{i=1}^{3N}\left[{\partial H \over \partial p_i} \left (-{\partial H \over \partial q_i} + D_iF_e (t) \right ) +{\partial H \over \partial q_i}\left({\partial H \over \partial p_i} + C_i F_e(t)\right)\right] \\[4pt] &={\partial f_0 \over \partial H}\sum_{i=1}^{3N}\left[D_i({\rm x}) {\partial H \over \partial p_i} + C_i({\rm x}){\partial H \over \partial q_i}\right]F_e(t) \end{align*}

Define

$j({\rm x}) = -\sum_{i=1}^{3N}\left[D_i({\rm x}){\partial H \over \partial p_i} + C_i({\rm x}){\partial H \over \partial q_i}\right] \nonumber$

which is known as the dissipative flux. Thus, for a Cartesian Hamiltonian

$H = \sum_{i=1}^N {{\textbf p}_i^2 \over 2m_i} + U({\textbf r}_1,...,{\textbf r}_N) \nonumber$

where $${\textbf F}_i({\textbf r}_1,...,{\textbf r}_N) = -\nabla_iU$$ is the force on the $${i}$$ th particle, the dissipative flux becomes:

$j({\rm x}) = \sum_{i=1}^N\left[{\bf C}_i({\rm x})\cdot{\textbf F}_i-{\textbf D}_i({\rm x})\cdot {{\textbf p}_i \over m_i}\right] \nonumber$

In general,

$\dot{\rm x}\cdot\nabla_{\rm x}f_0({\rm x}) = -{\partial f_0 \over \partial H}j({\rm x})F_e(t) \nonumber$

Now, suppose $$f_0 (x)$$ is a canonical distribution function

$f_0(H({\rm x})) = {1 \over Q(N,V,T)}e^{-\beta H({\rm x})} \nonumber$

then

${\partial f_0 \over \partial H} = -\beta f_0(H) \nonumber$

so that

$\dot{\rm x}\cdot\nabla_{\rm x}f_0({\rm x}) = \beta f_0({\rm x})j({\rm x})F_e(t) \nonumber$

Thus, the solution for $$\Delta f (x, t)$$ is

$\Delta f({\rm x},t) = -\beta\int_0^t ds e^{-iL_0(t-s)}f_0({\rm x})j({\rm x})F_e(s) \nonumber$

The ensemble average of the observable $$A (x)$$ now becomes

\begin{align*} \langle A (t) \rangle &=\langle A \rangle_0 - \beta\int d{\rm x}A({\rm x})\int_0^t ds e^{-iL_0(t-s)}f_0({\rm x})j({\rm x})F_e(s) \\[4pt] &=\langle A \rangle_0 - \beta\int_0^t ds \int d{\rm x}A({\rm x})e^{-iL_0(t-s)}f_0({\rm x})j({\rm x})F_e(s) \\[4pt] &=\langle A \rangle_0 - \beta\int_0^t ds \int d{\rm x}f_0({\rm x}) A({\rm x})e^{-iL_0(t-s)}j({\rm x})F_e(s) \end{align*}

Recall that the classical propagator is $$exp(iLt)$$. Thus the operator appearing in the above expression is a classical propagator of the unperturbed system for propagating backwards in time to $$- (t - s)$$. An observable $$A (x)$$ evolves in time according to

\begin{align*} \dfrac{dA}{dt} &= iLA \\[4pt] A (t) &= e^{iLt}A(0) \\[4pt] A (-t ) &= e^{-iLt}A(0) \end{align*}

Now, if we take the complex conjugate of both sides, we find

$A^*(t) = A^*(0)e^{-iLt} \nonumber$

where now the operator acts to the left on $$A^* (0)$$. However, since observables are real, we have

$A(t) = A(0) e^{-iLt} \nonumber$

which implies that forward evolution in time can be achieved by acting to the left on an observable with the time reversed classical propagator. Thus, the ensemble average of $$A$$ becomes

\begin{align*} \langle A(t)\rangle &=\langle A \rangle_0 - \beta\int_0^t ds F_e(s) \int d{\rm x}_0 f_0({\rm x}_0)A({\rm x}_{t-s}({\rm x}_0))j({\rm x}_0) \\[4pt] &=\langle A \rangle_0 - \beta\int_0^t ds F_e(s)\langle j(0)A(t-s)\rangle_0 \end{align*}

where the quantity on the last line is an object we have not encountered yet before. It is known as an equilibrium time correlation function. An equilibrium time correlation function is an ensemble average over the unperturbed (canonical) ensemble of the product of the dissipative flux at $$t = 0$$ with an observable $$A$$ evolved to a time $${t - s}$$. Several things are worth noting:

1. The nonequilibrium average $$\langle A(t)\rangle$$, in the linear response regime, can be expressed solely in terms of equilibrium averages.
2. The propagator used to evolve $$A (x)$$ to $$A (x, t - s )$$ is the operator $$exp (iL_0 (t - s ) )$$, which is the propagator for the unperturbed, Hamiltonian dynamics with $$C_i = D_i = 0$$. That is, it is just the dynamics determined by $$H$$.
3. Since $$A({\rm x},t-s) = A({\rm x}(t-s))$$ is a function of the phase space variables evolved to a time $${t - s }$$, we must now specify over which set of phase space variables the integration $$\int dx$$ is taken. The choice is actually arbitrary, and for convenience, we choose the initial conditions. Since $$x (t)$$ is a function of the initial conditions $$x (0)$$, we can write the time correlation function as $\langle j(0)A(t-s)\rangle_0 = {1 \over Q}\int d{\rm x}_0 e^{-\beta H({\rm x}_0)}j({\rm x}_0)A({\rm x}_{t-s}({\rm x}_0)) \nonumber$

This page titled 12.5: Perturbative solution of the Liouville equation is shared under a CC BY-NC-SA 4.0 license and was authored, remixed, and/or curated by Mark Tuckerman.