Reference: Fetter and Walecka, Quantum Theory of Many Particle Systems, Ch. 1
The Hamiltonian for a SHO is:
$$ H = \sum_{i = 0}^{\infty}\hbar \omega ( a_i^{+} a_i + \frac{1}{2} ) $$
where $\{a^+_i, a_i\}$ are the creation and annihilation operators for the $i^\textrm{th}$ eigenstate (momentum mode). The Fock space $\mathbf{F}$ consists of states of the form:
$$ \vert n_{a_0},n_{a_1}, ...,n_{a_N} \rangle $$
which are obtained by repeatedly acting on the vacuum $\vert 0 \rangle $ by the ladder operators:
$$ \Psi = \vert n_{i_0},n_{i_1}, ...,n_{i_N} \rangle = (a_0^+)^{i_0} (a_1^+)^{i_1} \ldots (a_N^+)^{i_N} \vert 0 \rangle $$
The interpretation of $\Psi$ is as the state which contains $i_k$ quanta of the $k^\textrm{th}$ eigenstate created by application of $(a^+_k)^{i_k}$ on the vacuum.
The above state is not normalized until multiplied by factor of the form $\prod_{k=0}^N \frac{1}{\sqrt{k+1}}$. If your excitations are bosonic you are done, because the commutator of the ladder operators $[a^+_i,a_j] = \delta_{ij}$ vanishes for $i\ne j$. However if the statistics of your particles are non-bosonic (fermionic or anyonic) then the order, in which you act on the vacuum with the ladder operators, matters.
Of course, to construct a Fock space $\mathbf{F}$ you do not need to specify a Hamiltonian. Only the ladder operators with their commutation/anti-commutation relations are needed. In usual flat-space problems the ladder operators correspond to our usual fourier modes $ a^+_k \Rightarrow \exp ^{i k x} $. For curved spacetimes this can procedure can be generalized by defining our ladder operators to correspond to suitable positive (negative) frequency solutions of a laplacian on that space. For details, see Wald, QFT in Curved Spacetimes.
Now, given any Hamiltonian of the form:
$$ H = \sum_{k=1}^{N} T(x_k) + \frac{1}{2} \sum_{k \ne l = 1}^N V(x_k,x_l) $$
with a kinetic term $T$ for a particle at $x_k$ and a pairwise potential term $V(x_k,x_l)$, one can write down the quantum Hamiltonian in terms of matrix elements of these operators:
$$ H = \sum_{ij} a^+_i \langle i \vert T \vert j \rangle a_i + \frac{1}{2}a^+_i a^+_j \langle ij \vert V \vert kl \rangle a_l a_k $$
where $|i\rangle$ is the state with a single excited quantum corresponding the action of $a^+_i$ on the vacuum. (For details, steps, see Fetter & Walecka, Ch. 1).
I hope this helps resolves some of your doubts. Being as you are from math, there are bound to be semantic differences between my language and yours so if you have any questions at all please don't hesitate to ask.
(Not really an answer, but as one should not state such things in comments, I'm putting it here)
You commented: "This seems to boil down to the relationship between the phase space and the Hilbert space."
That's a deep question. I recommend reading Urs Schreiber's excellent post on how one gets from the phase space to the operators on a Hilbert space in a natural fashion. I'm not certain how the Wigner/Moyal picture of QM relates to quantum statistical mechanics, since we define the quantum canonical partition function to be $Z(\beta) := \mathrm{Tr}(\mathrm{e}^{\beta H})$ on the Hilbert space of states, as we basically draw the analogy that the classical phase space is the "space of states" for our classical theory, and the integral the trace over it, and generalize that to the quantum theory.
Also note that, in a quantum world, $\int\mathrm{d}x\mathrm{d}p\mathrm{e}^{-\beta H}$ is a bit of a non-sensical expression, since $H$ is an operator - the result of this would not be a number, which the partition function certainly should be.
Best Answer
The hypothetical balls are part of a single quantum system, i.e., there can be (and indeed are) quantum mechanical correlations between them.
If the system is in a state representing a single particle, then it is known that only one ball is excited, but it is uncertain which ball it is.
For each ball, there is a probability amplitude that it is the one that is excited. If you write a function for the probability amplitude that the ball at a particular position is excited, that gives you the quantum wavefunction of the particle.