The presumed equivalence between the canonical quantization and the Fock space representation is only a particular case.
The canonical formalism provides only with canonical Poisson brackets. The first step according to Dirac's axioms is to replace the Poisson brackets by commutators and since these commutators satisfy the Jacobi identity, they can be represented by linear operators on a Hilbert space.
Canonical quantization does not specify the Hilbert space.
Finding a Hilbert space where the operators acts linearly and satisfy the commutation relations is a problem in representation theory. This task is referred to as "quantization" in the modern literature.
The problem is that in the case of free fields, this problem does not have a unique solution (up to a unitary transformation in the Hilbert space). This situation is referred to as the existence of inequivalent quantizations or inequivalent representations. The Fock representation is only a special case. Some of the quantizations are called "non-Fock", because the Hilbert space does not have an underlying Fock space structure (i.e., cannot be interpreted as free particles), but there can even be inquivalent Fock representations.
Before, proceeding, let me tell you that inequivalent quantizations may be the areas where "new physics" can emerge because they can correspond to different quantum systems.
Also, let me emphasize, that the situation is completely different in the finite dimensional case. This is because that due to the Stone-von Neumann theorem, any representation of the canonical commutation relations in quantum mechanics is unitarily equivalent to the harmonic oscillator representation. Thus the issue of inquivalent representations of the canonical commutation relations occurs only due to the infinite dimensionality.
For a few examples of inquivalent quantizations of the canonical commutation relations of a scalar field on a Minkowski space-time, please see the following article by: Moschella and Schaeffer. In this article, they construct inequivalent representations by means of Bogoliubov transformation which changes the vacuum and they also present a thermofield representation. In all these representations the canonical operators are represented on a Hilbert space and the canonical commutation relations are satisfied. The Bogoliubov shifted vacuum cases correspond to broken Poincare' symmetries. One can argue that these solutions are unphysical, but the symmetry argument will not be enough in the case of quantization on a general curved nonhomogenous manifold. In this case we will not have a "physical" argument to dismiss some of the inequivalent representations.
The phenomena of inequivalent quantizations can be present even in the case of finite number of degrees of freedom on non-flat phase spaces.
Having said all that, I want nevertheless to provide you a more direct answer to your question (although it will not be unique due to the reasons listed above). As I understand the question, it can be stated that there is an algorithm for passing from the single particle Hilbert space to the Fock space. This algorithm can be summarized by the Fock factorization:
$$ \mathcal{F} = e^{\otimes \mathcal{h}}$$
Where $\mathcal{h}$ is the single particle Hilbert space and $\mathcal{F}$ is the Fock space. As stated before canonical quantization provides us only with the canonical commutaion relations:
$$[a_{\mathbf{k}}, a^{\dagger}_{\mathbf{l}}] = \delta^3(\mathbf{k} - \mathbf{l}) \mathbf{1}$$
At this stage we have only an ($C^{*}$)algebra of operators. The reverse question about the existence of an algorithm starting from the canonical commutation relations and ending with the Fock space (or equivalently, the answer to the question where is the Hilbert space?) is provided by the Gelfand -Naimark-Segal construction (GNS), which provides representations of $C^{*}$ algebras in terms of bounded operators on a Hilbert space.
The GNS construction starts from a state $\omega$ which is a positive linear functional on the algebra $ \mathcal{A}$ (in our case the algebra is the completion of all possible products of any number creation and annihilation operators).
The second step is choosing the whole algebra as an initial linear space $ \mathcal{A}$. In general, there will be null elements satisfying:
$$\omega (A^{\dagger}{A}) = 0$$
The Hilbert space is obtained by identifying elements differing by a null vector:
$$ \mathcal{H} = \mathcal{A} / \mathcal{N} $$
($\mathcal{N} $ is the space of null vectors).
The inner product on this Hilbert space is given by:
$$(A, B) = \omega (A^{\dagger}{B}) $$
It can be proved that the GNS construction is a cyclic representation where the Hilbert space is given by the action of operators on a cyclic "vacuum vector". The GNS construction gives all inequivalent representations of a given $C^{*}$ algebra (by bounded operators). In the case of a free scalar field the choice of a Gaussian state defined by its characteristic function:
$$ \omega_{\mathcal{F} }(e^{\int\frac{d^3k}{E_k} z_{\mathbf{k}}a^{\dagger}_{\mathbf{k}} + \bar{z}_{\mathbf{k}}a^{\mathbf{k}} }) = e^{\int\frac{d^3k}{E_k} \bar{z}_{\mathbf{k}} z_{\mathbf{k}}}$$
Where $z_{\mathbf{k}}$ are indeterminates which can be differentiated by to obtain the result for any product of operators.
The null vectors of this construction will be just combinations vanishing due to the canonical commutation relations (like $a_1 a_2 - a_2 a_1$). Thus this choice has Bose statistics. Also subspaces spanned by a product of a given number of creation operators will be the number subspaces.
The state of this specific construction is denoted by:
$\omega_{\mathcal{F}}$, since it produces the usual Fock space. Different state choices may result inequivalent quantizations.
@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
You are right. The $i$ indices of the occupation numbers $n_i$ can be thought of as multi-indices meaning that each $i$ encodes both the wavevector $k$ (which is quantized once boundary conditions are imposed) and the spin sign $\sigma$. In general, $i$ will enumerate all the allowed combinations of quantum numbers that characterize the single particle states.
It is also true that only the entries $n_i$ of those states which are populated will be non-zero. If the total particle number $N$ is fixed, the constraint $\sum_in_i=N$ has to hold.
I added this in response to your comment:
The space the occupation number states naturally live in is referred to as Fock space. It can be written as
$\mathcal F^\pm = \bigoplus_{N=0}^\infty\mathcal H_N^\pm = \mathcal H_0 \oplus \mathcal H_1 \oplus \mathcal H_2^\pm \oplus \dots$
where $\mathcal H_N^\pm$ refers to the usually $N$ particle Hilbert spaces and the $\pm$ indicates whether we restict ourselves to symmetric or antisymmetric states, i.e. bosons and fermions respectively. It is convenient to work in Fock space because annihilation and creation operators will take you from one Hilbert space with fixed particle number to another.
Now, assuming that you mean 'dimension' when you talk about the 'size' of a Hilbert space, keep in mind that you have to obey the (anti-)symmetry constaint. This reduces the dimension of each of the $N$ particle subspaces.
As an example, consider a two state fermionic system, e.g. electrons with spin either $\downarrow$ ($i=1$) or $\uparrow$ ($i=2$). Due to the Pauli principle (i.e. as a consequence of the antisymmetry requirement) the only occupation numbers allowed are 0 and 1. $\mathcal H_0$ is one-dimensional as always, it contains only the vacuum state $|00\rangle$ (and multiples thereof, but we are interested in normalized linear independent states). The one particle Hilbert space is two-dimensional as the electron can be in either one of the states (and symmetry still doesn't matter since there are no two particles to swap). The two particle Hilbert space $\mathcal H_2^-$ has again only one linear independent state, $|11\rangle$, because Mr. Pauli doesn't allow us to populate any of the two states twice. Thus, the two-state fermionic Fock space is four-dimenional.
More generally, the $M$-state fermionic Fock space will be $\sum_{N=0}^M\binom{M}{N}=2^M$ dimensional (which follows from simple combinatorics) and bosonic Fock spaces are infinite dimensional because bosons can occupy the same state as many times as they want.