I'm a computer scientist! So let's consider the finite dimensional case first. Let's try n=4, even. Suppose you have a particle that can be in one 4 states, say, |a>, |b>, |c>, and |d>. Pick any ol' 4x4 hermitian matrix you like for the hamiltonian. This says how much the particle likes being in different states, in the sense of how much enegery it has in them. (At thermal equilibrium, more energetic things are less likely; things like dropping down to lower energy states) For example I could say |a> is 5 Joules, |b> is 17, |c> is -3, |d> is 0. But this is quantum physics, so I could also pick some other basis: I could say (|a>+|b>)/\2 is 50 Joules, (|a>-|b>)/\2 is 1/2, (|c>+3i|d>)/\10 is 3.9, (3i|c>+|d>) is 100.

But no matter what I say, it's hermitian, so it has 4 nice orthonormal real-eigenvalued eigenvectors, and the eigenvectors are the states that have a definite energy, and the eigenvalues are the corresponding energies. Call the eigenvectors U_1, ... U_4. Their coordinates in terms of the basis vectors |a>, |b>, |c>, and |d> I'm going to write like U_{1a}, U_{2a}, etc., so that U_n = U_{na}|a> + U_{nb}|b> + U_{nc}|c> + U_{nd}|d>.

I'm going to now consider a larger state space, but still finite, of dimension 2^4, namely the corresponding fermionic fock space. You can have anywhere between 0 and 4 particles, and any of |a>, |b>, |c>, and |d> may be present or absent, independently. There's 16 basis states, some examples of which are |>, |a>, |ac>, |abc>, |bcd>, |abcd>. The only trick is that the permutation of the labels actually matters, (this is the "fermion" part) and swapping any two letters yields the same state only with the label flipped. So |ab> = -|ba>, |ac> = -|ca>, etc.

I'm going to define the dynamics --- i.e. define a Hamiltonian --- on the big Fock space in terms of the one I already picked for the single-particle case. What are the eigenvectors of the big Hamiltonian? Tensor products --- either 0-ary or 1-ary or binary or 3-ary or whatever you want --- of eigenvectors of the one-particle Hamiltonian, subject to order-flipping proviso. And the energies of these 'composite' eigenvectors are just the sums of the energies that went into them.

For an example, imagine the 1-particle hamiltonian with

eigenvector (|a>+|b>)/\2, eigenvalue 50

eigenvector (|a>-|b>)/\2, eigenvalue 1

eigenvector |c>, eigenvalue 9

eigenvector |d>, eigenvalue 0

Then a few examples of eigenvectors in Fock space are:

If we take (|a>+|b>)/\2 tensor |c>, we get an eigenvector (|ac>+|bc>)/\2 with eigenvalue 50+9=59

If we take (|a>+|b>)/\2 tensor (|a>-|b>)/\2, we get an eigenvector (|a>+|b>)(|a>-|b>)/2 = (|aa>+|ba>-|ab>-|bb>)/2 = (|ba>+|ba>)/2 = |ba> with eigenvalue 50+1=51

If we take (|a>-|b>)/\2 tensor |c> tensor |d> we get (|acd>-|bcd>)/\2 with eigenvalue 9+1+0=10

Now we have 16 orthogonal energy eigenstates; we let the system come to thermal equilibrium and sample it in the |a>,|b>,|c>,|d> basis. What does this mean? We pick an energy eigenstate at random with probability inversely proportional to the exponential of its energy. Then we observe some subset of |a>, |b>, |c>, |d> according to whatever the Born rule says about that eigenstate.

I think the punchline is: this whole game ends up being a determinantal point process with an L-ensemble kernel which has the same eigenvectors as the Hamiltonian, but the nth eigenvalue is not the energy E_n, but rather exp(-E_n) (or, like, exp(-E_n/β) if you want to actually think about the Boltzmann distribution temperature).