Real scalar field and its quantum state: why are the diagonal components static here? 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.
 A: To answer your comment on the simplest NR QM  pure state illustration you chose,
$$
\rho (t)= \frac{1}{2} (e^{-it\omega/2}|0\rangle +e^{-it3\omega/2}|1\rangle)(e^{it\omega/2}\langle 0| +e^{it3\omega/2}\langle 1 |),
$$
express it in the coordinate basis, yielding Hermite functions,
$$
\rho_{xx'}(t)= \frac{1}{2} \left ( \psi_0(x)  +e^{-it\omega} \psi_1(x)\right )~~\left ( \psi^*_0(x')+e^{it\omega} \psi_1^*  (x')   \right  ),
$$
so $x=x'$  collapses to a real  function, a probability density.
Its trace, if you will, in this basis, is the $\int dx \rho=1$,   naturally yielding unity. The Landau-Lifshitz NR QM text covers this nicely.
Your point in the comment that this collapsed real function is time dependent
$$
\frac{1}{2} ( |\psi_0(x)|^2 + |\psi_1(x)|^2+ \cos (\omega t)   \psi_0(x)   \psi_1(x) ),
$$
is well put.  It oscillates with time, as it should. You took the $x'\to x$ limit incorrectly.
The vanishing of the r.h.side of this equation as $x'\to x$ is not sound, given differential operators... there are singular terms omitted. Here is a formal wisecrack possibly evoking this:
$$
(\partial_x - \partial_{x'}) e^{ic x} e^{-ic x'}= 2ic  e^{ic (x-x')}.
$$
But, according to your $x'\to x$  limit the left hand side vanishes, whereas the right hand side not, $=2ic$.  Your vanishing quantity is a mirage.
But... you really have seen this before in the diagonal sector, with the probability current! $\rho(x,x)=\rho(x)= |\psi(x)|^2$, the probability density, which is not time invariant at all: probability flows!
$$
\partial_t \rho(x) = \partial_x\left (  \frac{i\hbar}{2m} (\psi^*(x)\partial_x \psi(x) - \psi(x) \partial_x \psi^*(x) )  \right )\\ 
= \frac{i\hbar}{2m}  \left ( \psi^*(x)\partial_x^2 \psi(x) - \psi(x) \partial_x^2 \psi^*(x)   \right ) \neq 0 !
$$
For your flip-flop wave function, you may confirm on both sides that it is the oscillating quantity
$$
-\omega\sin(\omega t) \psi_0 (x) \psi_1(x) .
$$

*

*Takeaway: the diagonal components are not static. The limit $(\partial_x^2-\partial^2_{x'})\rho_{xx'}\neq 0$, when taken carefully. First the derivatives operate on "half" of $\rho_{xx'}$, with the time variation encoded in the imaginary part, and only then is $x'\to x$ taken.


*It is the probability current inflow. So the premise of the question, "why", is off: should be missing.
A: It is not true that $\dot\rho(\phi_1,\phi_1)=0$: the object $\rho(\phi_1,\phi_1)$ does not exist.
It doesn't make sense to talk about $\rho$ for $\phi_1=\phi_2$. The object $\rho$ is a distribution (as a function of $\phi$, not $x$) so it is meaningless to "evaluate" it on a given function.
Indeed,
$$
\langle\phi_1|\phi_2\rangle=\delta(\phi_1-\phi_2)
$$
where $\delta$ is the Dirac delta in the space of functionals (i.e., with respect to functional integration, not regular integration). See this PSE post for some (formal) applications of this formula.
The density matrix $\rho(\phi_1,\phi_2)$ only makes sense under the integral sign:
$$
\int_{\mathcal F^2}\rho(\phi_1,\phi_2)F[\phi_1,\phi_2]\,\mathrm d\phi_1\mathrm d\phi_1
$$
where $\mathcal F$ denotes some reasonable function space and $F$ is a healthy functional on it.
