Consider a (free and massless) real scalar field $\phi(x)$ with Hamiltonian $$ H := \int d^3\mathbf{x}\; \bigg[ \frac{1}{2} \pi^2(\mathbf{x}) + |\nabla\phi(\mathbf{x})|^2 \bigg] $$ where $\pi(\mathbf{x})$ is the conjugate momentum field (where we have $[\phi(\mathbf{x}), \pi(\mathbf{y}) ] = i \delta^{(3)}(\mathbf{x} - \mathbf{y})$). We work in the Schrodinger picture here.
The density matrix for the field $\rho(t)$ satisfies the von Neumann equation $$ \frac{d \rho(t)}{d t} = - i [ H, \rho(t) ] \ . $$ Let $| \phi_a \rangle$ be the eigenkets of the field operator $\phi(\mathbf{x})$ with eigenvalue $\phi_a(\mathbf{x})$. This means $$ \hat{\phi}(\mathbf{x}) | \phi_a \rangle = \phi_a(\mathbf{x}) | \phi_a \rangle \ . $$ From here we will derive an equation of motion for the components of the density matrix in the field basis, where $$ \rho_{\phi_1 \phi_2}(t) := \langle \phi_1 | \rho(t) | \phi_2 \rangle \ . $$ By using the fact that for any state $| \chi \rangle$ we have $$ \langle \phi_1 | \pi(\mathbf{x}) | \chi \rangle = - i \frac{\delta}{\delta \phi_1(\mathbf{x})} \langle \phi_1 | \chi \rangle \ , $$ it is straightforward to use the von Neumann equation to derive the equation for the off-diagonal components $\phi_1 \neq \phi_2$ that $$ \frac{d \rho_{\phi_1\phi_2}}{dt} = - \frac{i}{2} \int d^3\mathbf{x} \left( - \frac{\delta^2}{\delta \phi_1(\mathbf{x})^2} + |\nabla \phi_1(\mathbf{x})|^2 + \frac{\delta^2}{\delta \phi_2(\mathbf{x})^2} - |\nabla \phi_2(\mathbf{x})|^2 \right) \rho_{\phi_1\phi_2}(t) \ . $$ However strangely for diagonal components $\phi_1 = \phi_2$ I am finding that you get $$ \frac{d \rho_{\phi_1\phi_1}}{dt} = 0 \ . $$ What is the meaning of this? This seems to suggest that the diagonal components are always constant, no matter what bizarre state the field starts out in. This somehow seems wrong.