This is the first time in several years that I've done a post about physics that didn't have to do with my research. This came about from thinking about applying techniques in statistical physics to game theory; although I still have a lot more to learn about that and need to do more to flesh out those ideas, it occurred to me in the process that I never had such a good intuition for the phase space density in classical mechanics, and notes that I've found online focus almost exclusively on the phase space density of a large number of particles in an explicitly statistical treatment. I intend to use this post to shed light on why this may be the case, help build intuition for how things like the Liouville equation work for simple systems of one or a few particles, and reinforce the notion that there is no classical analogue to the phenomenon of a multi-particle entangled quantum state yielding a mixed single-particle state under a partial trace. Follow the jump to see more.
Consider the case of a single particle whose Hamiltonian is generically given by \(H(x, p)\) and whose phase space density is posited to be \(\rho(t, x, p) = \delta(x - X(t))\delta(p - P(t))\) where \((X(t), P(t))\) are functions of time \(t\) but are independent of the phase space coordinates \((x, p)\); this should ultimately be self-consistent as the phase space density corresponding to a deterministic phase space trajectory \((X(t), P(t))\). The Liouville equation is \(\frac{\partial \rho}{\partial t} = -\frac{\partial \rho}{\partial x} \frac{\partial H}{\partial p} + \frac{\partial \rho}{\partial p} \frac{\partial H}{\partial x}\). Evaluating the partial time derivative gives \(\frac{\partial \rho}{\partial t} = -\dot{X}(t) \delta'(x - X(t))\delta(p - P(t)) - \dot{P}(t) \delta(x - X(t))\delta'(p - P(t))\), where \(\delta'\) is the derivative of the Dirac delta function with respect to its argument; the only parts of the phase space density that explicitly depend on \(t\) are the functions \((X(t), P(t))\), so their time derivatives (which can be written as total time derivatives) multiply the Dirac delta functions and their derivatives. Evaluating the partial phase space coordinate derivatives gives \(\frac{\partial \rho}{\partial x} = \delta'(x - X(t))\delta(p - P(t))\) and \(\frac{\partial \rho}{\partial p} = \delta(x - X(t))\delta'(p - P(t))\). This overall yields a single equation of motion \(-\dot{X}(t) \delta'(x - X(t))\delta(p - P(t)) - \dot{P}(t) \delta(x - X(t))\delta'(p - P(t)) = \) \( -\delta'(x - X(t))\delta(p - P(t))\frac{\partial H}{\partial p} + \delta(x - X(t))\delta'(p - P(t))\frac{\partial H}{\partial x}\). This yields the usual Hamiltonian equations of motion as follows. If an integration is performed over \(x\), then the result is multiplied by \(p\), and then that is integrated over \(p\), the result is \(\dot{P}(t) = -\frac{\partial H}{\partial x}\). If instead an integration is performed over \(p\), then the result is multiplied by \(x\), and then that is integrated over \(x\), the result is \(\dot{X}(t) = \frac{\partial H}{\partial p}\). Importantly in each case, the partial derivatives of the Hamiltonian \(\frac{\partial H}{\partial x}\) and \(\frac{\partial H}{\partial p}\) are evaluated making the replacements \((x, p) \to (X(t), P(t))\). These functions of time, once solved for, are plugged back into the definition of the phase space density. Thus, it is rare for a treatment of the classical dynamics of a single deterministic particle to consider the phase space density, because it yields nothing new compared to the standard Hamiltonian equations of motion, and it just adds more complicated mathematical baggage to yield the same results. As an example, it is easy to check that the standard results are recovered for the harmonic oscillator Hamiltonian \(H(x, p) = \frac{p^{2}}{2m} + \frac{kx^{2}}{2}\).
Consider the case of two particles whose Hamiltonian is generically given by \(H(x_{1}, p_{1}, x_{2}, p_{2})\) and whose phase space density is posited to be \(\rho(t, x_{1}, p_{1}, x_{2}, p_{2}) = \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\) where \((X_{1}(t), P_{1}(t), X_{2}(t), P_{2}(t))\) are functions of time \(t\) but are independent of the phase space coordinates \((x_{1}, p_{1}, x_{2}, p_{2})\); this should ultimately be self-consistent as the phase space density corresponding to a deterministic phase space trajectory \((X_{1}(t), P_{1}(t), X_{2}(t), P_{2}(t))\). The Liouville equation is \(\frac{\partial \rho}{\partial t} = -\frac{\partial \rho}{\partial x_{1}} \frac{\partial H}{\partial p_{1}} + \frac{\partial \rho}{\partial p_{1}} \frac{\partial H}{\partial x_{1}} - \frac{\partial \rho}{\partial x_{2}} \frac{\partial H}{\partial p_{2}} + \frac{\partial \rho}{\partial p_{2}} \frac{\partial H}{\partial x_{2}}\). Evaluating the partial time derivative gives \(\frac{\partial \rho}{\partial t} = -\dot{X}_{1}(t) \delta'(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \( - \dot{P}_{1}(t) \delta(x_{1} - X_{1}(t))\delta'(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{X}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta'(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{P}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta'(p_{2} - P_{2}(t))\). The partial phase space coordinate derivatives are straightforward, analogous to the case for one particle with appropriate replacements. Therefore, the Liouville equation becomes \(-\dot{X}_{1}(t) \delta'(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \( - \dot{P}_{1}(t) \delta(x_{1} - X_{1}(t))\delta'(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{X}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta'(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{P}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta'(p_{2} - P_{2}(t))\) \( = -\delta'(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\frac{\partial H}{\partial p_{1}} \) \(+ \delta(x_{1} - X_{1}(t))\delta'(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\frac{\partial H}{\partial x_{1}} \) \(- \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta'(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\frac{\partial H}{\partial p_{2}} \) \(+ \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta'(p_{2} - P_{2}(t))\frac{\partial H}{\partial x_{2}} \). Integrating this equation over \(p_{1}\), \(x_{2}\), and \(p_{2}\), then multiplying by and integrating over \(x_{1}\) yields \(\dot{X}_{1}(t) = \frac{\partial H}{\partial p_{1}}\). Integrating this equation over \(x_{1}\), \(x_{2}\), and \(p_{2}\), then multiplying by and integrating over \(p_{1}\) yields \(\dot{P}_{1}(t) = -\frac{\partial H}{\partial x_{1}}\). Integrating this equation over \(x_{1}\), \(p_{1}\), and \(p_{2}\), then multiplying by and integrating over \(x_{2}\) yields \(\dot{X}_{2}(t) = \frac{\partial H}{\partial p_{2}}\). Integrating this equation over \(x_{1}\), \(p_{1}\), and \(x_{2}\), then multiplying by and integrating over \(p_{2}\) yields \(\dot{P}_{2}(t) = -\frac{\partial H}{\partial x_{2}}\). These functions of time, once solved for, are plugged back into the definition of the phase space density.
Once again, it seems like this has yielded no useful results beyond what can be derived through standard Hamiltonian mechanics, and only produces more mathematical baggage. However, there are questions of interpretation for multi-particle states that are not available for single-particle states. In particular, it is clear that the single-particle phase space density can be obtained as \( \int \rho(t, x_{1}, p_{1}, x_{2}, p_{2}) \operatorname{d}x_{2} \operatorname{d}p_{2} = \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t)) \), which is exactly as expected for a deterministic classical trajectory for a particle among a deterministic collection therein. This is very different from the quantum case, where taking the partial trace over an entangled pure state in the basis of separable states yields a mixed state for the remaining particles. This reinforces the notion that there is no way to reproduce quantum entanglement using a classical multi-particle phase space density. (I'm vaguely aware of work from the last several years about correlations between droplets of water carried by surface waves being mathematically analogous to quantum entanglement as predicted by the deBroglie-Bohm pilot wave interpretation of quantum mechanics, but don't know enough about it to comment further.) The fact that the single-particle classical state is still separable holds despite the fact that the equations of motion for \((X_{1}(t), P_{1}(t))\) may have some complicated dependence on \((X_{2}(t), P_{2}(t))\); those are all simply functions of time, and are not the phase space coordinates defining the classical state of the system, and this can be seen using the example of two coupled harmonic oscillators.
A friend of mine sent an easier validation of this classical separability for the case of two identical coupled harmonic oscillators. Skipping many steps, it is possible to write the phase space density generically as \( \rho(t, x_{+}, p_{+}, x_{-}, p_{-}) = \delta(x_{+} - X_{+}(t))\delta(p_{+} - P_{+}(t))\delta(x_{-} - X_{-}(t))\delta(p_{-} - P_{-}(t))\), having defined the phase space coordinates \(x_{\pm} = (x_{1} \pm x_{2})/\sqrt{2}\) and \(p_{\pm} = (p_{1} \pm p_{2})/\sqrt{2}\). Using the rules of integration over Dirac delta functions, it can be shown that the single-particle phase space density is \( \int \rho(t, (x_{1} + x_{2})/\sqrt{2}, (p_{1} + p_{2})/\sqrt{2}, (x_{1} - x_{2})/\sqrt{2}, (p_{1} - p_{2})/\sqrt{2}) \operatorname{d}x_{2} \operatorname{d}p_{2} \) \( = \delta(x_{1} - (X_{+}(t) + X_{-}(t))/\sqrt{2})\delta(p_{1} - (P_{+}(t) - P_{-}(t))/\sqrt{2}) \). Not only is this a valid single-particle phase space density (helped in part by Liouville's theorem, guaranteeing the constancy of infinitesimal phase space volumes under coordinate transformations as well as time evolution), but it exactly yields the deterministic dynamics of particle 1; the linear combinations only affect the functions of time and do not yield a nontrivial probabilistic spread of states for particle 1 from the pure classical state of the two-particle system.
Single Particle
Consider the case of a single particle whose Hamiltonian is generically given by \(H(x, p)\) and whose phase space density is posited to be \(\rho(t, x, p) = \delta(x - X(t))\delta(p - P(t))\) where \((X(t), P(t))\) are functions of time \(t\) but are independent of the phase space coordinates \((x, p)\); this should ultimately be self-consistent as the phase space density corresponding to a deterministic phase space trajectory \((X(t), P(t))\). The Liouville equation is \(\frac{\partial \rho}{\partial t} = -\frac{\partial \rho}{\partial x} \frac{\partial H}{\partial p} + \frac{\partial \rho}{\partial p} \frac{\partial H}{\partial x}\). Evaluating the partial time derivative gives \(\frac{\partial \rho}{\partial t} = -\dot{X}(t) \delta'(x - X(t))\delta(p - P(t)) - \dot{P}(t) \delta(x - X(t))\delta'(p - P(t))\), where \(\delta'\) is the derivative of the Dirac delta function with respect to its argument; the only parts of the phase space density that explicitly depend on \(t\) are the functions \((X(t), P(t))\), so their time derivatives (which can be written as total time derivatives) multiply the Dirac delta functions and their derivatives. Evaluating the partial phase space coordinate derivatives gives \(\frac{\partial \rho}{\partial x} = \delta'(x - X(t))\delta(p - P(t))\) and \(\frac{\partial \rho}{\partial p} = \delta(x - X(t))\delta'(p - P(t))\). This overall yields a single equation of motion \(-\dot{X}(t) \delta'(x - X(t))\delta(p - P(t)) - \dot{P}(t) \delta(x - X(t))\delta'(p - P(t)) = \) \( -\delta'(x - X(t))\delta(p - P(t))\frac{\partial H}{\partial p} + \delta(x - X(t))\delta'(p - P(t))\frac{\partial H}{\partial x}\). This yields the usual Hamiltonian equations of motion as follows. If an integration is performed over \(x\), then the result is multiplied by \(p\), and then that is integrated over \(p\), the result is \(\dot{P}(t) = -\frac{\partial H}{\partial x}\). If instead an integration is performed over \(p\), then the result is multiplied by \(x\), and then that is integrated over \(x\), the result is \(\dot{X}(t) = \frac{\partial H}{\partial p}\). Importantly in each case, the partial derivatives of the Hamiltonian \(\frac{\partial H}{\partial x}\) and \(\frac{\partial H}{\partial p}\) are evaluated making the replacements \((x, p) \to (X(t), P(t))\). These functions of time, once solved for, are plugged back into the definition of the phase space density. Thus, it is rare for a treatment of the classical dynamics of a single deterministic particle to consider the phase space density, because it yields nothing new compared to the standard Hamiltonian equations of motion, and it just adds more complicated mathematical baggage to yield the same results. As an example, it is easy to check that the standard results are recovered for the harmonic oscillator Hamiltonian \(H(x, p) = \frac{p^{2}}{2m} + \frac{kx^{2}}{2}\).
Two particles
Consider the case of two particles whose Hamiltonian is generically given by \(H(x_{1}, p_{1}, x_{2}, p_{2})\) and whose phase space density is posited to be \(\rho(t, x_{1}, p_{1}, x_{2}, p_{2}) = \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\) where \((X_{1}(t), P_{1}(t), X_{2}(t), P_{2}(t))\) are functions of time \(t\) but are independent of the phase space coordinates \((x_{1}, p_{1}, x_{2}, p_{2})\); this should ultimately be self-consistent as the phase space density corresponding to a deterministic phase space trajectory \((X_{1}(t), P_{1}(t), X_{2}(t), P_{2}(t))\). The Liouville equation is \(\frac{\partial \rho}{\partial t} = -\frac{\partial \rho}{\partial x_{1}} \frac{\partial H}{\partial p_{1}} + \frac{\partial \rho}{\partial p_{1}} \frac{\partial H}{\partial x_{1}} - \frac{\partial \rho}{\partial x_{2}} \frac{\partial H}{\partial p_{2}} + \frac{\partial \rho}{\partial p_{2}} \frac{\partial H}{\partial x_{2}}\). Evaluating the partial time derivative gives \(\frac{\partial \rho}{\partial t} = -\dot{X}_{1}(t) \delta'(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \( - \dot{P}_{1}(t) \delta(x_{1} - X_{1}(t))\delta'(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{X}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta'(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{P}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta'(p_{2} - P_{2}(t))\). The partial phase space coordinate derivatives are straightforward, analogous to the case for one particle with appropriate replacements. Therefore, the Liouville equation becomes \(-\dot{X}_{1}(t) \delta'(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \( - \dot{P}_{1}(t) \delta(x_{1} - X_{1}(t))\delta'(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{X}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta'(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t)) \) \(- \dot{P}_{2}(t) \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta'(p_{2} - P_{2}(t))\) \( = -\delta'(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\frac{\partial H}{\partial p_{1}} \) \(+ \delta(x_{1} - X_{1}(t))\delta'(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\frac{\partial H}{\partial x_{1}} \) \(- \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta'(x_{2} - X_{2}(t))\delta(p_{2} - P_{2}(t))\frac{\partial H}{\partial p_{2}} \) \(+ \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t))\delta(x_{2} - X_{2}(t))\delta'(p_{2} - P_{2}(t))\frac{\partial H}{\partial x_{2}} \). Integrating this equation over \(p_{1}\), \(x_{2}\), and \(p_{2}\), then multiplying by and integrating over \(x_{1}\) yields \(\dot{X}_{1}(t) = \frac{\partial H}{\partial p_{1}}\). Integrating this equation over \(x_{1}\), \(x_{2}\), and \(p_{2}\), then multiplying by and integrating over \(p_{1}\) yields \(\dot{P}_{1}(t) = -\frac{\partial H}{\partial x_{1}}\). Integrating this equation over \(x_{1}\), \(p_{1}\), and \(p_{2}\), then multiplying by and integrating over \(x_{2}\) yields \(\dot{X}_{2}(t) = \frac{\partial H}{\partial p_{2}}\). Integrating this equation over \(x_{1}\), \(p_{1}\), and \(x_{2}\), then multiplying by and integrating over \(p_{2}\) yields \(\dot{P}_{2}(t) = -\frac{\partial H}{\partial x_{2}}\). These functions of time, once solved for, are plugged back into the definition of the phase space density.
Once again, it seems like this has yielded no useful results beyond what can be derived through standard Hamiltonian mechanics, and only produces more mathematical baggage. However, there are questions of interpretation for multi-particle states that are not available for single-particle states. In particular, it is clear that the single-particle phase space density can be obtained as \( \int \rho(t, x_{1}, p_{1}, x_{2}, p_{2}) \operatorname{d}x_{2} \operatorname{d}p_{2} = \delta(x_{1} - X_{1}(t))\delta(p_{1} - P_{1}(t)) \), which is exactly as expected for a deterministic classical trajectory for a particle among a deterministic collection therein. This is very different from the quantum case, where taking the partial trace over an entangled pure state in the basis of separable states yields a mixed state for the remaining particles. This reinforces the notion that there is no way to reproduce quantum entanglement using a classical multi-particle phase space density. (I'm vaguely aware of work from the last several years about correlations between droplets of water carried by surface waves being mathematically analogous to quantum entanglement as predicted by the deBroglie-Bohm pilot wave interpretation of quantum mechanics, but don't know enough about it to comment further.) The fact that the single-particle classical state is still separable holds despite the fact that the equations of motion for \((X_{1}(t), P_{1}(t))\) may have some complicated dependence on \((X_{2}(t), P_{2}(t))\); those are all simply functions of time, and are not the phase space coordinates defining the classical state of the system, and this can be seen using the example of two coupled harmonic oscillators.
A friend of mine sent an easier validation of this classical separability for the case of two identical coupled harmonic oscillators. Skipping many steps, it is possible to write the phase space density generically as \( \rho(t, x_{+}, p_{+}, x_{-}, p_{-}) = \delta(x_{+} - X_{+}(t))\delta(p_{+} - P_{+}(t))\delta(x_{-} - X_{-}(t))\delta(p_{-} - P_{-}(t))\), having defined the phase space coordinates \(x_{\pm} = (x_{1} \pm x_{2})/\sqrt{2}\) and \(p_{\pm} = (p_{1} \pm p_{2})/\sqrt{2}\). Using the rules of integration over Dirac delta functions, it can be shown that the single-particle phase space density is \( \int \rho(t, (x_{1} + x_{2})/\sqrt{2}, (p_{1} + p_{2})/\sqrt{2}, (x_{1} - x_{2})/\sqrt{2}, (p_{1} - p_{2})/\sqrt{2}) \operatorname{d}x_{2} \operatorname{d}p_{2} \) \( = \delta(x_{1} - (X_{+}(t) + X_{-}(t))/\sqrt{2})\delta(p_{1} - (P_{+}(t) - P_{-}(t))/\sqrt{2}) \). Not only is this a valid single-particle phase space density (helped in part by Liouville's theorem, guaranteeing the constancy of infinitesimal phase space volumes under coordinate transformations as well as time evolution), but it exactly yields the deterministic dynamics of particle 1; the linear combinations only affect the functions of time and do not yield a nontrivial probabilistic spread of states for particle 1 from the pure classical state of the two-particle system.