# 9. Classical and quantum dynamics of density matrices


Statistical mechanics makes the connection between macroscopic dynamics and equilibriums states based on microscopic dynamics. For example, while thermodynamics can manipulate equations of state and fundamental relations, it cannot be used to derive them. Statistical mechanics can derive such equations and relations from first principles.

Before we study statistical mechanics, we need to introduce the concept of the density, referred to in classical mechanics as density function, and in quantum mechanics as density operator or density matrix. The key idea in statistical mechanics is that the system can have “microstates,” and these microstates have a probability. For example, there may be a certain probability that all gas atoms are in a corner of the room, and this is probably much lower than the probability that they are evenly distributed throughout the room. Statistical mechanics deals with these probabilities, rather than with individual particles. Three general contexts of probability are in common use:

a. Discrete systems: and example would be rolling dice. A die has 6 faces, each of which is a “microstate.” Each outcome is equally likely if the die is not loaded, so

$p_{i}=\dfrac{1}{W}=\dfrac{1}{6}$

is the probability of being in microstate “i” where W=6 is the total number of microstates.

b. Classical systems: here the microstate has to be specified by the positions xi and the momenta pi of all particles in the system, so $$\rho\left ( x_{i},p_{i},t \right )$$ is the time dependent probability of finding the particles at xi and pi. Do not confuse the momentum here with the probability in a.! It should be clear from the context. r is the classical density function.

c. Quantum systems: here the microstate is specified by the density operator or density matrix $$\hat{\rho}\left ( t \right )$$. The probability that the system is in quantum state “i” with state |i> is given by

$p_{i}(t)=\int d x \Psi_{i}^{*}(x) \hat{\rho}(x, t) \Psi_{i}(x)=\langle i|\hat{\rho}(t)| i\rangle$.

Of course the probability does not have to depend on time if we are in an equilibrium state.

In all three cases, statistical mechanics attempts to evaluate the probability from first principles, using the Hamiltonian of the closed system. In an equilibrium state, the probability does not depend on time, but still depends on x and p (classical) or just x (quantum).

We need to briefly review basic concepts in classical and quantum dynamics to see how the probability evolves in time, and when it does not evolve (reach equilibrium).

Mechanics: classical

Definition of phase space: Phase space is the 6n dimensional space of the 3n coordinates and 3n momenta of a set of n particles, which, taken together constitute the system.

Definition of a trajectory: The dynamics of a system of N degrees of freedom (N = 6n for n particles in 3-D) are specified by a trajectory {xj(t), pj(t)}, j = 13n in phase space.

Note: a system of n particles is phase space is defined by a single point {x1(t), x2(t), ... x3n(t), p1(t), p2(t), ... p3n(t)} that evolves in time. For that specific system, the density function is a delta function in time centered at the single point, and moving along in time as the phase spacve trajectory trajectory moves along. The phase space trajectory is not to be confused with the 3-D trajectories of individual particles.

In statistical mechanics, the system must satisfy certain constraints: (e.g. all xi(t) must lie within a box of volume V; all speeds must be less than the speed of light; etc). The density function $$\rho\left ( x_{i},p_{i},t_{i};constraints \right )$$ that we usually care about in statistical mechanics is the average probability density of ALL systems satisfying the constraints. We call the group of systems satisfying the same constraints an “ensemble”, and the average density matrix the “ensemble denisty matrix”:

$\rho=\dfrac{1}{W} \sum_{i=1}^{W} \rho_{i}$,

where the sum is over all W possible systems in the microstates “m” satisfying the constraint. Unlike an individual rm, which is a moving spike in phase space, r sums over all microstates and looks continuous (at least after a small amount of smoothing). For example, consider a single atom with position and momentum {x, p} in a small box. Each different rm for each microstate “m” is a spike moving about in the phase space {x, p}. Averaging over all such microstates yields a r that is uniformly spread over the positions in the box (independent of position x), with a Maxwell-Boltzmann distribution of velocities: $$\rho \sim e^{-p^{2}/2mk_{B}T}$$ (assuming the walls of the box can equilibrate the particle to a temperature T). How does r evolve in time?

Each coordinate qi = xi evolves according to Newton’s law, which can be recast as

$\begin{gathered} F_{i}=m \ddot{x}_{i} \\ -\dfrac{\partial V}{\partial x_{i}}=\dfrac{d}{d t}\left(m \dot{x}_{i}\right)=\dfrac{d}{d t}\left(\dfrac{\partial K}{\partial \dot{x}_{i}}\right) \end{gathered}$

if the force is derived from a potential V and where

$K=\dfrac{1}{2} \sum_{i} m_{i} \dot{x}_{i}^{2}=\dfrac{1}{2} \sum_{i} \dfrac{p_{i}^{2}}{2 m_{i}}$

is the kinetics energy. We can define the Lagrangian L = KV and rewrite the above equation as

$\dfrac{d}{d t}\left(\dfrac{\partial L\left(x_{i}, \dot{x}_{i}, t\right)}{\partial \dot{x}_{i}}\right)-\dfrac{\partial L}{\partial x_{i}}=0$.

One can prove using variational calculus (see appendix A) that this differential equation (Lagrange’s equation) is valid in any coordinate system, and is equivalent to the statement

$\min \left\{S=\int_{t_{0}}^{t_{f}} L\left(x_{i}, \dot{x}_{i}\right) d t\right\}$

where S is the action (not to be confused with the entropy!). Let’s say the particle moves from from positions xi(t0) at t=t0 to xi(tf) at tf. Guess a trajectory xi(t). The trajectory xi(t) can be used to compute velocity $$= )\partial x_{i}/\partial t =\dot{x} \;and\;L\left ( x_{i},\dot{x}_{i} \right )$$. The actual trajectory followed by the classical particles is the one that minimizes the above integral, called “the action.”

$\dfrac{\partial L}{\partial \dot{x}_{i}}=p_{i} \quad\left(p_{i}=m \dot{x}_{i}\right)$

Thinking of $$L$$ as $$L\left(\dot{x}_{i}\right)$$ and of $$p_{i}$$ as the derivatives, we can Legendre transform to a new representation $$H$$

$-H \equiv L-\sum_{i=1}^{N} \dot{x}_{i} p_{i} .$

It will become obvious shortly why we define $$H$$ with a minus sign. According to the rules for Legendre transforms,

$\dfrac{\partial H}{\partial p_{i}}=\dot{x}_{i} \text { and } \dfrac{\partial H}{\partial \dot{x}_{i}}=p_{i}$

These are Hamilton's equation of motion for a trajectory in phase space. They are equivalent to solving Newton's equation. Evaluating $$H$$,

\begin{aligned} &H\left(\dot{x}_{i}, p_{i}\right)=-K+V+\sum_{i} \dfrac{p_{i}}{m_{i}} p_{i} \\ &=K+V \end{aligned}

The Hamiltonian is sum of kinetic and potential energy, i.e. the total energy, and is conserved if $$H$$ is not explicitly time-dependent:

$\dfrac{d H(x, p)}{d t}=\dfrac{\partial H}{\partial x} \dfrac{d x}{d t}+\dfrac{\partial H}{\partial p} \dfrac{d p}{d t}=\dfrac{\partial H}{\partial x} \dfrac{d H}{d p}-\dfrac{\partial H}{\partial p} \dfrac{d H}{d x}=0 .$

Thus Newton's equations conserve energy. This is because all particles are accounted for in the Hamiltonian $$H$$ (closed system). Note that Lagrange's and Hamilton's equations hold in any coordinate system, so from now on we will write $$H\left(q_{i}, p_{i}\right)$$ instead of using $$x_{i}$$ (cartesian coordinates).

Let $$\hat{A}\left(q_{i}, p_{i}, t\right)$$ be any dynamical variable (many $$\hat{A} \mathrm{~s}$$ of interest do not depend explicitly on t, but we include it here for generality). Then

$\dfrac{d \hat{A}}{d t}=\sum_{i} \dfrac{\partial \hat{A}}{\partial q_{i}} \dfrac{d q_{i}}{d t}+\dfrac{\partial \hat{A}}{\partial p_{i}} \dfrac{d p_{i}}{d t}+\dfrac{\partial \hat{A}}{\partial t}=\sum_{i}\left(\dfrac{\partial \hat{A}}{\partial q_{i}} \dfrac{d H}{d p_{i}}-\dfrac{\partial \hat{A}}{\partial p_{i}} \dfrac{d H}{d q_{i}}\right)+\dfrac{\partial \hat{A}}{\partial t}=[\hat{A}, H]_{p}+\dfrac{\partial \hat{A}}{\partial t}$

gives the time dependence of $$\hat{A}$$. [] $$]_{\mathrm{P}}$$ is the Poisson bracket.

Now consider the ensemble density $$\rho\left ( x_{i},p_{i},t \right )$$ as a specific example of a dynamical variable. Because trajectories cannot be destroyed, we can normalize

$\iint d q_{i} d p_{i} \rho\left(q_{i}, p_{i}, t\right)=1$

Integrating the probability over all state space, we are guaranteed to find the system somewhere subject to the constraints. Since the above integral is a constant, we have

\begin{aligned} \dfrac{d}{d t} \iint d q_{i} d p_{i} \rho=0=& \iint d q_{i} d p_{i} \dfrac{d \rho}{d t}=\iint d q_{i} d p_{i}\left\{[\rho, H]+\dfrac{\partial \rho}{\partial t}\right\}=0 \\ & \Rightarrow \dfrac{\partial \rho}{\partial t}=-[\rho, H]_{p} \end{aligned}

This is the Liouville equation, it describes how the density propagates in time. To calculate the average value of an observable $$\hat{A}\left(q_{i}, p_{i}\right)$$ in the ensemble of systems described by $$\rho$$, we calculate

$A(t)=\iint d q_{i} d p_{i} \rho\left(q_{i}, p_{i}, t\right) \hat{A}\left(q_{i}, p_{i}\right)=\langle A\rangle_{\text {ens }}$

Thus if we know $$\rho$$, we can calculate any average observable. For certain systems which are left unperturbed by outside influences (closed systems)

$\lim _{t \rightarrow \infty} \dfrac{\partial \rho}{\partial t}=0 \text {. }$

$$\rho$$ reaches an equilibrium distribution $$\rho_{e q}\left(q_{i}, p_{i}\right)$$ and

$\left[\rho_{e q}, H\right]_{p}=0 \text { (definition of equilibrium) }$

In such a case, $$A(\mathrm{t}) \rightarrow \mathrm{A}$$, the equilibrium value of the observable. The goal of equilibrium statistical mechanics is to find the values A for a $$\rho$$ subject to certain imposed constraints; e.g. $$\rho_{e q}\left(q_{i}, p_{i}, U, V\right)$$ is the set of all possible system trajectories such that $$U$$ and $$V$$ are constant. The more general goal of non-equilibrium statistical mechanics is to find $$A(\mathrm{t})$$ given an initial condition $$\rho_{0}\left(q_{i}, p_{i}\right.$$; constraints). So much for the classical picture.

## Mechanics: quantum

Now let us rehearse the whole situation again for quantum mechanics. The quantum formulation is the one best suited to systems where the energy available to a degree of freedom becomes comparable or smaller than the characteristic energy gap of the degree of freedom. Classical and quantum formulations are highly analogous.
A fundamental quantity in quantum mechanics is the density operator

$\hat{\rho}_{i}(t)=\left|\psi_{i}(t)\right\rangle\left\langle\psi_{i}(t)|=| t\right\rangle\langle t| \text {. }$

This density operator projects onto the microstate "i" of the system, $$\left|\psi_{i}(t)\right\rangle$$. If we have an ensemble of W systems, we can define the ensemble density operator $$\hat{\rho}$$ subject to some constraints as

$\hat{\rho}(t ; \text { constraints })=\dfrac{1}{W} \sum_{I=1}^{\mathrm{W}} \hat{\rho}_{I}(t ; \text { constraints })$

For example, let the constraint be $$U=$$ const. Then we would sum over all microstates that are degenerate at the same energy $$U$$. This average is analogous to averaging the classical probability density over microstates subject to constraints. To obtain the equation of motion for $$\hat{\rho}$$, we first look at the wavefunction.

Its equation of motion is

$H \psi=i \hbar \dfrac{\partial}{\partial T} \psi$

which by splitting it into $$\psi_{r}$$ and $$i \psi_{i}$$, can be written

$\dfrac{1}{\hbar} \hat{H} \psi_{r}=\psi_{i} \text { and } \dfrac{1}{\hbar} \hat{H} \psi_{i}=-\psi_{r} .$

Using any complete basis $$H\left|\varphi_{j}\right\rangle=E_{j}\left|\varphi_{j}\right\rangle$$, the trace of $$\hat{\rho}$$ is conserved

$\operatorname{Tr}\left\{\rho_{i}(t)\right\}=\sum_{j}<\varphi_{j}\left|\psi_{i}(t)\right\rangle\left\langle\psi_{i}(t) \mid \varphi_{j}\right\rangle=\sum_{j}\left\langle\psi_{i} \mid \varphi_{j}\right\rangle\left\langle\varphi_{j} \mid \psi_{i}\right\rangle=\left\langle\psi_{i}(t) \mid \psi_{i}(t)\right\rangle=1$

if $$\psi_{j}(t)$$ is normalized, or

$\operatorname{Tr}\left\{\hat{\rho}_{i}\right\}=1 \text {, and } \operatorname{Tr}\{\hat{\rho}\}=1 \text {. }$

This basically means that probability density cannot be destroyed, in analogy to trajectory conservation. Note that if

$\begin{gathered} \hat{\rho}_{i}=|\psi\rangle\langle\psi| \\ \Rightarrow \hat{\rho}_{i}^{2}=\hat{\rho}_{i} \text { so } \operatorname{Tr}\left(\hat{\rho}_{i}^{2}\right)=1 \end{gathered}$

A state described by a wavefunction $$\mid \Psi>$$ that satisfies the latter equation is a pure state. Most states of interest in statistical mechanics are NOT pure states. If

$\hat{\rho}=\dfrac{1}{W} \sum_{i=1}^{W} \hat{\rho}_{i},$

the complex off-diagonal elements tend to cancel because of random phases and $$\operatorname{Tr}\left(\hat{\rho}^{2}\right)<1$$, this an impure state.

### Example:

Let $$\psi=c_{0}^{i}|0\rangle+c_{1}^{i}|1\rangle$$ be an arbitrary wavefunction for a two-level system.

$\Rightarrow \rho=|\Psi\rangle\left\langle\Psi\left|=c_{0}^{i} c_{0}^{i^{*}}\right| 0\right\rangle\left\langle 0\left|+c_{0}^{i} c_{1}^{i^{*}}\right| 0\right\rangle\left\langle 1\left|+c_{1}^{i} c_{0}^{i^{*}}\right| 1\right\rangle\left\langle 0\left|+c_{1}^{i} c_{1}^{i^{*}}\right| 1\right\rangle\langle 1|$

or, in matrix form,

$\hat{\rho}_{i}=\left(\begin{array}{cc}\left|c_{0}^{i}\right|^{2} & c_{0}^{i} c_{1}^{i^{*}} \\c_{0}^{i} c_{1}^{i^{*}} & \left|c_{1}^{i}\right|^{2} \end{array}\right) \quad T_{r}\left(\hat{\rho}_{i}\right)=\left|c_{0}^{i}\right|^{2}+\left|c_{1}^{i}\right|^{2}=1$

Generally, macroscopic constraints (volume, spin population, etc.) do not constrain the phases of $$c_{0}=\left|c_{0}\right| e^{i \varphi_{0}}$$ and $$c_{1}=\left|c_{1}\right| e^{i \varphi_{1}}$$. Thus, ensemble averaging

\begin{aligned} &\hat{\rho}=\lim _{W \rightarrow \infty} \dfrac{1}{W} \sum_{i=1}^{W} \hat{\rho}_{I}=\lim _{W \rightarrow \infty}\left(\begin{array}{cc} \left|c_{0}\right|^{2} & \dfrac{e^{i \bar{\varphi}}}{\sqrt{W}} \\ \dfrac{e^{-i \bar{\varphi}}}{\sqrt{W}} & \left|c_{1}\right|^{2} \end{array}\right)=\left(\begin{array}{cc} \left|c_{0}\right|^{2} & 0 \\ 0 & \left|c_{1}\right|^{2} \end{array}\right) \\ &\Rightarrow T_{r}\left(\hat{\rho}^{2}\right)=\left|c_{0}\right|^{4}+\left|c_{1}\right|^{4}<1 \end{aligned}

Such a state is known as an impure state.

To obtain the equation of motion of $$\hat{\rho}$$, use the time-dependent Schrödinger equation in propagator form,

$\left|\psi_{i}(t)\right\rangle=e^{-\dfrac{i}{\hbar} \hat{H} t}\left|\psi_{i}(0)\right\rangle$

and the definition of $$\rho$$,

$\hat{\rho}_{i}(t)=\left|\psi_{i}(t)\right\rangle\left\langle\psi_{i}(t)\right|$

to obtain the equation of motion

\begin{aligned} \dfrac{\partial}{\partial t} \hat{\rho}_{i}(t) &=\dfrac{\partial}{\partial t}\left\{e^{-\dfrac{i}{\hbar} \hat{H} t}\left|\psi_{i}(0)\right\rangle\left\langle\psi_{i}(0)\right| e^{+\dfrac{i}{\hbar} \hat{H} t}\right\} \\ &=-\dfrac{i}{\hbar} \hat{H} e^{-\dfrac{i}{\hbar} \hat{H} t}\left|\psi_{i}(0)\right\rangle\left\langle\psi_{i}(0)\left|e^{+\dfrac{i}{\hbar} \hat{H} t}+e^{-\dfrac{i}{\hbar} \hat{H} t}\right| \psi_{i}(0)\right\rangle\left\langle\psi_{i}(0)\right| e^{+\dfrac{i}{\hbar} \hat{H} t}\left(+\dfrac{i}{\hbar} \hat{H}\right) \\ &=-\dfrac{i}{\hbar} H \rho_{i}+\dfrac{i}{\hbar} \rho_{i} H \\ &=\dfrac{1}{i \hbar}\left[\rho_{i}, H\right] \end{aligned}

This is known as the Liouville-von Neumann equation. The commutator defined in the last line is equivalent to the Poisson bracket in classical dynamics. Summing over all microstates to obtain the average density operator

$\dfrac{1}{W} \sum_{i=1}^{W} \hat{\rho}_{i} \Rightarrow \dfrac{\partial \hat{\rho}}{\partial t}=-\dfrac{1}{i \hbar}[\hat{\rho}, \hat{H}]$

This von-Neumann equation is the quantum equation of motion for $$\hat{\rho}$$. If $$\rho$$ represent an impure state, this propagation cannot be represented by the time-dependent Schrödinger equation.

We are interested in average values of observables in an ensemble of systems. Starting with a pure state,

$A(t)=\left\langle\psi_{i}(t)|\hat{A}| \psi_{i}(t)\right\rangle=\sum_{j}\left\langle\psi_{i}(t) \mid \varphi_{j}\right\rangle\left\langle\varphi_{j}|\hat{A}| \psi_{i}(t)\right\rangle=\sum_{j}\left\langle\varphi_{j}|\hat{A}| \psi_{i}(t)\right\rangle\left\langle\psi_{i}(t) \mid \varphi_{j}\right\rangle=T_{r}\left\{\hat{A} \hat{\rho}_{i}(t)\right\}$

Summing over ensembles,

$\dfrac{1}{W} \sum_{i=1}^{W} T_{r}\left\{\hat{A} \rho_{i}\right\}=T_{r}\left\{\hat{A} \dfrac{1}{W} \sum_{i=1}^{W} \hat{\rho}_{i}\right\} \Rightarrow A(t)=T_{r}\{\hat{A} \hat{\rho}(t)\}$

In particular, $$P_{j}=T_{r}\left\{\hat{P}_{j} \hat{\rho}(t)\right\}=\operatorname{Tr}\left\{\left|\varphi_{j}\right\rangle\left\langle\varphi_{j}\right| \hat{\rho}(t)\right\}=\left\langle\varphi_{j}|\hat{\rho}(t)| \varphi_{j}\right\rangle=\rho_{j j}(t) \quad$$ is $$\quad$$ the probability of being in state $$\mathrm{j}$$ at time $$\mathrm{t}$$.

Finally, $$\hat{\rho}$$ may evolve to long-time solutions $$\hat{\rho}_{e q}$$ such that $$\dfrac{\partial \rho}{\partial t}=0 \Rightarrow$$

$\left[\hat{\rho}_{e q}, H\right]=0 \text{(condition for equilibrium)}$.

In that case, the density matrix has relaxed to the equilibrium density matrix, which no longer evolves in time.

### Example:

Consider a two-level system again. Let

$H|j\rangle=E|j\rangle \text { for } j=0,1 \text { or } H=\left(\begin{array}{cc} E_{0} & 0 \\ 0 & E_{1} \end{array}\right) ; \text { if } \hat{p}=\left(\begin{array}{cc} \rho_{00} & 0 \\ 0 & \rho_{11} \end{array}\right) \Rightarrow[\hat{\rho}, H]=0$

Thus a diagonal density matrix of a closed system does not evolve. The equilibrium density matrix must be diagonal so $$\left[\rho_{e q}, H\right]=0$$. This corresponds to a completely impure state. Note that unitary evolution cannot change the purity of any closed system. Thus, the density matrix of a single closed system cannot evolve to diagonality unless
$\hat{\rho}=\dfrac{1}{W} \sum_{i} \rho_{i}$
for the ensemble is already diagonal. In reality, single systems still decohere because they are open to the environment: let i denote the degrees of freedom of the system, and $$\mathrm{j}$$ of a bath (e.g. a heat reservoir). $$\hat{\rho}$$ depends on both $$i$$ and $$j$$ and can be written as a matrix $$\hat{\rho}_{i j, i^{1} j^{1}}(t)$$. For example, for a two-level system coupled to a large bath, indices i,j only go from 1 to 2 , but indices i' and j' could go to $$10^{20}$$. We can average over the bath by letting

$\hat{\rho}^{(r e d)}=\operatorname{Tr}_{j}\{\hat{\rho}\}$

$$\hat{\rho}^{(r d)}$$ only has matrix elements for the system degrees of freedom, e.g. for a two level system in contact with a bath of $$10^{20}$$ states, $$\hat{\rho}^{(c e d)}$$ is still a $$2 \times 2$$ matrix.

We will show later that for a bath is at constant $$T$$.

$\hat{\rho}^{(r e d)} \rightarrow \rho_{e q}=\left(\begin{array}{cc} \dfrac{e^{-E_{0} / k_{B} T}}{Q} & 0 \\ 0 & \dfrac{e^{-E_{1} / k_{B} T}}{Q} \end{array}\right)$

Note that "reducing" was not necessary in the classical discussion because quantum coherence and phases do not exist there. To summarize analogous entities:

Classical Quantum
$${\rho\left(q_{i}, p_{i}\text { constr. }\right)}$$ $${\hat{\rho}(t ; \text { constr. })}$$ Density function or operator
$$\hat{A}\left(q_{i}, p_{i}\right)$$ operator $$\hat{A}$$ Observable, hermition
$$\dfrac{\partial}{\partial p} H=\dot{q} \dfrac{\partial}{\partial q} H=-\dot{p}$$ $$\dfrac{1}{\hbar} H \psi_{r}=\dot{\psi}_{i} \dfrac{1}{\hbar} H \psi_{i}=-\dot{\psi}_{r}$$

Eq. of motion for traj. ψ

$$\text { Averaging } \\ \iint d p_{i} d q_{i}$$ $$T_{r}\{\}$$ Averaging
$$\iint d p_{i} d q_{i} \rho=1$$ $$T_{r}\{\hat{\rho}\}=1$$ Conservation of probability
$$A(t)=\iint d p_{i} d q_{i} \hat{A}\left(q_{i}, p_{i}\right) \rho\left(q_{i}, p_{i}, t\right)$$ $$A(t)=T_{r}(\hat{A}, \hat{\rho})$$ Expectation value in ensemble
$$\dfrac{\partial \rho}{\partial t}=-[\rho, H]_{P}$$ $$\dfrac{\partial \hat{\rho}}{\partial t}=\dfrac{1}{i \hbar}[\rho, H]$$ Equation of motion for p
$${\left[\rho_{e q}, H\right]_{P}=0}$$ $${\left[\hat{\rho}_{e q}, H\right]=0}$$ Necessary equlibrium condition (closed system)

To illustrate basic ideas, we will often go back to simple discrete models with probability for each microstate $$p_{i}$$, but to get accurate answers, one may have to work with the full classical or quantum probability $$\rho$$.

9. Classical and quantum dynamics of density matrices is shared under a CC BY-NC-SA 4.0 license and was authored, remixed, and/or curated by LibreTexts.