I was reading through this paper on thermal transport in crystals, and saw that a primary function of their interest was $n_{\mu}(\mathbf{x},t)$, which they state is the phonon excitation/occupation number at position $\mathbf{x}$ and time $t$, where $\mu=(\mathbf{k},s)$ is the index of the phonon ($\mathbf{k}$=wavenumber, $s$=polarization).
How can you have a phonon occupation number at a particular location and time? I thought phonons were lattice excitations that necessarily extended across the entire crystal, i.e. they were excitations of a sort of "phonon field". I know phonons are bosons, and that for a crystal in thermal equilibrium with its environment the phonon occupation number will be given by standard Bose-Einstein occupation number factor:
$$\bar{n} (\epsilon)=\frac{1}{e^{(\epsilon-\mu)/T}-1}$$
where $\epsilon$ is the energy of the phonon, $\mu$ is the chemical potential, and $T$ is the temperature.
I understand that this is for out of equilibrium dynamics, and that typically one describes the relaxation of the spacial distribution (density) of constituent particles. Quantum-mechanically, even "classical" particles extend across the entire containing volume just as phonons do (through their wavefunction), but in the macroscopic limit we can faithfully assign a position and momentum to each of them without much loss of information. However, particles are collections/superpositions of harmonic modes of the system (think particle in a box), conspiring in such a way to give a classically definite position and momentum. Phonons, on the other hand, are the harmonic excitations themselves.
So,... how?