I think it might be conceptually easier to start backwards. Let us assume we have $N$ decoupled harmonic oscillators:
$$
H = \sum_{a=1}^N \left[ \frac{(\pi_a)^2}{2m} + \frac{1}{2} m \omega_i^2 \phi_a^2 \right]
$$
Now imagine applying the unitary transformation $x_j = {U_j}^a \phi_a$ and $p_j = {U_j}^a \pi_a$ where $U^T U = 1$ in order that both $[x,p]$ and $[\phi,\pi]$ have the standard commutation relation. Such a transformation will act simply on the $\pi_i^2$ term, but the $\omega_i$'s in the second contribution to $H$ will lead to interesting quadratic couplings between the $x_j$. (If you further want unequal masses (I'm not sure you really do in the context of entanglement entropy of a free scalar field), you can perform the scaling $x_i \to \sqrt{m_i} x_i$ and $p_i \to p_i / \sqrt{m_i}$ without changing the canonical commutation relation between the $x$'s and $p$'s.)
Now the question is about a very particular nearest neighbor coupling between the $x_i$. The usual way to think about this transformation involves Fourier series. If I'm not mistaken, something like
$$
x_j = \sum_{a=1}^N e^{2 \pi i j a /N} \phi_a
$$
should act to diagonalize your Hamiltonian if I further make the assumption that $x_{N+1} = x_1$, i.e. I'm living on a circle.
In the context of the Srednicki paper and entanglement entropy, there are some more modern approaches involving two-point functions. You might want to look at section 2.2 of the review http://arxiv.org/abs/0905.2562 or section III of http://arxiv.org/abs/0906.1663.
@Andrew's answer provides the big picture, but I'd like to give a few more specific pointers that hopefully may help.
Questions 1-2: Does everything look correct so far? How to express these operators in terms of one-particle operators?
So you want to set up single-particle analogues of the ladder operators using eigenstates of the 1st quantized Hamiltonian, then use those to construct the second quantization ladder operators in the symmetrized multi-particle space. The problem is that such a procedure may not be possible. Here is why:
The entire second quantization framework rests on an isomorphism between the (anti)symmetric subspace of the N-particle Hilbert space and an abstract direct product of "mode" Hilbert spaces, each constructed around its own ladder operator algebra. The ladder operators ${\hat a}_n$, ${\hat a}^\dagger_n$ in the abstract/"mode" Hilbert space obviously do have equivalents on the original N-particle (anti)symmetric subspace. But cannot be expressed as symmetrized sums of similar single-particle operators. To see why, let us suppose the ${\hat a}_n$, ${\hat a}^\dagger_n$ can indeed be expressed as such symmetrized sums, reading
$$
{\hat a}_n \sim \sum_{k=1}^N{ {\hat \alpha}^{(k)}_n},\;\;\; {\hat a}^\dagger_n \sim \sum_{k=1}^N{ \left( {\hat \alpha}^{(k)}_n\right)^\dagger}
$$
up to some suitable normalization factor, where ${\hat \alpha}^{(k)}_n$, $\left( {\hat \alpha}^{(k)}_n\right)^\dagger$ are the desired "single-particle ladder operators" for particle $k$ and "mode"/eigenstate $n$, and each term is to be understood in the sense of ${\hat \alpha}^{(k)}_n\otimes \left[\bigotimes_{j\neq k}{\hat I}^{(j)}\right]$. The latter can be naturally assumed to (anti)commute for different particles and eigenstates (pre-symmetrization), so $\left[{\hat \alpha}^{(k)}_m, {\hat \alpha}^{(j)}_n\right]_\pm = \left[{\hat \alpha}^{(k)}_m, \left( {\hat \alpha}^{(j)}_n\right)^\dagger\right]_\pm = 0$ for any $j\neq k$ and for $n\neq m$ when $j = k$. Then the (anti)commutation relations for the ${\hat a}_n$-s,
$$
\left[{\hat a}_m, {\hat a}_n\right]_\pm = 0,\;\;\; \left[ {\hat a}_m, {\hat a}^\dagger_n\right]_\pm = \delta_{mn}{\hat I}
$$
require
$$
0 = \left[{\hat a}_m, {\hat a}_n\right]_\pm \sim \sum_{k=0}^N{\left[{\hat \alpha}^{(k)}_m, {\hat \alpha}^{(k)}_n\right]_\pm }, \;\;\;\;
\delta_{mn}{\hat I} = \left[{\hat a}_m, {\hat a}^\dagger_n\right]_\pm \sim \sum_{k=0}^N{\left[{\hat \alpha}^{(k)}_m,\left( {\hat \alpha}^{(k)}_n\right)^\dagger\right]_\pm }
$$
The important point here is that in either case the lhs does not depend on $N$. Then the first eq. above implies that each term on the rhs must vanish identically, which is great. But things are no longer clear cut for the 2nd eq. Whatever prescription we might propose for $\left[{\hat \alpha}^{(k)}_m,\left( {\hat \alpha}^{(k)}_n\right)^\dagger\right]_\pm $, there is no way to normalize the sum such that the result is non-zero yet independent of $N$ for any $N$.
Bottom line: however intuitive it may seem at first sight, this is not the way to go.
Question 3: In QFT why do we treat those excitations as multi-particle states?
Short answer: Due to the isomorphism with the many-particle framework. Sometimes, as in solid state physics, the excitations are referred to as quasi-particles for this exact reason. Think phonons and excitons. Same goes for photons and any other field quanta, but for historical reasons they are referred to as "particles".
Question 4: Can we 'second-quantise' the EM field by treating the Maxwell equations as a Schrödinger equation (e.g. see books by Fushchich and Nikitin) and then considering the multi-particle states of those?
Unfortunately I'm not familiar with the book you mention, but you may want to google "Maxwell equations in Dirac form", for instance this paper and this Wikipedia page (especially refs. within). Never saw this used as a starting point for second quantization though, but why not? Perhaps an interesting idea?
Question 5: Why are the two procedures equivalent and lead to the same result?
Another short answer, same as for Question 3: Because both procedures rely on an isomorphism with the same type of abstract Hilbert space and the associated operator algebra. As mentioned before, the "Classical Mechanics"/"particle" procedure observes an isomorphism between the finite (anti)symmetric subspace of the N-particle Hilbert space and a "fixed particle number" subspace of the abstract second quantized space, while the "Field Theory" procedure yields a true isomorphism (the "postulate" you mention).
Best Answer
Indirectly.
To start with, set the pesky and useless constants $\hbar/m\omega \mapsto 1 $, to work in natural units. You then recognize the matrix, as Eric emphasizes by identifying to his (5), as written in the number basis of oscillator discrete energy eigenstates $|n\rangle$. Now for a celebrated basis change. Infinity should not bother you.
Consider the state $$ |q\rangle\equiv \sum_n \psi^*_n(x) |n\rangle , $$ where $\psi_n(x)= \langle x|n\rangle $, are (real) Hermite functions, eigenfunctions of the oscillator Schroedinger operator in the x representation-- dot with $\langle m|$: get it? Dotting with $\langle y|$ and exploiting the completeness of said Hermite functions, $$\langle y|q\rangle =\sum_n \psi^*_n(x) \psi_n(y)=\delta(x-y), $$ so you identify this state with $|q\rangle=|x\rangle$.
The pros normally use simpler operators instead of complete sums; that is, they sum the sum to $$ |x\rangle =\frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2}|0\rangle , $$ so that , as you should check, $$ \hat{x} |x\rangle= \frac{(a+a^\dagger)}{\sqrt 2} ~ \frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle= x~\left ( \frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle \right ) =x|x\rangle ~. $$
So the answer to your question is the eigenvalues are all continuous values x, for eigenvectors $e^{x^2/2} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle /\pi^{1/4} $.
To check your algebra, try the evident unnormalized eigenvector corresponding to x=0, $(1,0,-1/\sqrt{2},0,\sqrt{3/8},0,-\sqrt{5}/4,0,...)$, associated to Hermite numbers.
For a quick crib-sheet on the number basis ops cf Messiah, Quantum Mechanics vI Ch XII §5.