How does the normalization of an $n$-particle state $|n_{\mathbf{k}}\rangle$ work?

Question

You can expand the free, real scalar field in the following manner $$ \phi(x) = \int \frac{d^{3}\mathbf{k}}{(2\pi)^3} \frac{1}{\sqrt{2\omega_{\mathbf{k}}}} \bigg[ e^{- i \omega_{\mathbf{k}}x^0+i \mathbf{k} \cdot \mathbf{x}} a_{\mathbf{k}} + e^{+ i \omega_{\mathbf{k}}x^0-i \mathbf{k} \cdot \mathbf{x}} a_{\mathbf{k}}^{\dagger} \bigg] $$ where $\omega_{\mathbf{k}} = \sqrt{ |\mathbf{k}|^2 + m^2 }$ and the operators $a_{\mathbf{k}}$ and $a_{\mathbf{k}}^{\dagger}$ obey the bosonic canonical commutation relations (CCRs): $$ [a_{\mathbf{k}},a_{\mathbf{p}}]=[a_{\mathbf{k}}^{\dagger},a_{\mathbf{p}}^{\dagger}]=0 \ \ \ \ \ \ \ , \ \ [a_{\mathbf{k}},a_{\mathbf{p}}^{\dagger}] = \delta(\mathbf{k}-\mathbf{p}) $$

From here you build a symmetric Fock space as follows: You define a vacuum state $|0\rangle$ such that: $$ a_{\mathbf{k}} |0\rangle = 0 $$ which is assumed to be normalized $\langle0|0\rangle=1$. Then a single-particle state can be defined as $$ | 1_{\mathbf{k}} \rangle \equiv a_{\mathbf{k}}^{\dagger}|0\rangle $$ which are orthonormal $\langle 1_{\mathbf{k}} | 1_{\mathbf{p}} \rangle = \delta(\mathbf{k} - \mathbf{p})$ on account of the CCRs. From here you can construct a state containing $\sum_{s=1}^{N} n_s$ particles: $$ | (n_{1})_{\mathbf{k}_1} (n_{2})_{\mathbf{k}_2} \ldots (n_{N})_{\mathbf{k}_N} \rangle \equiv \frac{ \big( a_{\mathbf{k}_1}^{\dagger} \big)^{n_{1}} }{\sqrt{n_1!}}\frac{ \big( a_{\mathbf{k}_2}^{\dagger} \big)^{n_{2}} }{\sqrt{n_2!}} \cdots \frac{ \big( a_{\mathbf{k}_N}^{\dagger} \big)^{n_{N}} }{\sqrt{n_N!}} | 0 \rangle $$ where there are $n_s$ particles with the momentum $\mathbf{k}_{s}$.

I am confused about the way the following states are normalized: $$ | n_{\mathbf{k}} \rangle = \frac{ \big( a_{\mathbf{k}}^{\dagger} \big)^{n} }{\sqrt{n!}} | 0 \rangle $$ What is the normalization property of $\langle m_{\mathbf{p}} | n_{\mathbf{k}}\rangle$? Or also $\langle n_{\mathbf{p}} | n_{\mathbf{k}}\rangle$? For example, I find that \begin{align} \langle 2_{\mathbf{p}} | 2_{\mathbf{k}} \rangle &= \tfrac{1}{2} \langle 0 | a_{\mathbf{p}} a_{\mathbf{p}} a^{\dagger}_{\mathbf{k}} a^{\dagger}_{\mathbf{k}} | 0 \rangle\\ &= \tfrac{1}{2} \langle 0 | a_{\mathbf{p}} \big[ a^{\dagger}_{\mathbf{k}} a_{\mathbf{p}} + \delta(\mathbf{p} - \mathbf{k}) \big] a^{\dagger}_{\mathbf{k}} | 0 \rangle\\ &= \tfrac{1}{2} \langle 0 | a_{\mathbf{p}} a^{\dagger}_{\mathbf{k}} a_{\mathbf{p}} a^{\dagger}_{\mathbf{k}} | 0 \rangle + \tfrac{1}{2} \delta(\mathbf{p} - \mathbf{k}) \langle 0 | a_{\mathbf{p}} a^{\dagger}_{\mathbf{k}} | 0 \rangle\\ &= \tfrac{1}{2} \langle 0 | a_{\mathbf{p}} a^{\dagger}_{\mathbf{k}} \big[ a^{\dagger}_{\mathbf{k}} a_{\mathbf{p}} + \delta(\mathbf{p} - \mathbf{k}) \big] | 0 \rangle + \tfrac{1}{2} \delta(\mathbf{p} - \mathbf{k}) \langle 1_{\mathbf{p}}| 1_{\mathbf{k}} \rangle\\ &= 0 + \tfrac{1}{2} \delta(\mathbf{p} - \mathbf{k}) \langle 1_{\mathbf{p}} | 1_{\mathbf{k}} \rangle + \tfrac{1}{2} \delta(\mathbf{p} - \mathbf{k}) \langle 1_{\mathbf{p}}| 1_{\mathbf{k}} \rangle\\ &= \delta(\mathbf{p} - \mathbf{k})^{2} \end{align}

This makes no sense, as the square of the Dirac delta function is ill-defined. What am I doing wrong? How does the normalization actually work here?

EDIT: I have been using definitions for particle states that correspond to discrete labels - it is wrong to use them for continuous state labels. You cannot have more than one particle with exactly the same momentum. As usual Weinberg (Chapter 4) says it the clearest, he calls an $N$-particle state; $$ | \mathbf{k}_1\mathbf{k}_2 \cdots \mathbf{k}_N \rangle = a_{\mathbf{k}_1}^{\dagger} a_{\mathbf{k}_2}^{\dagger} \cdots a_{\mathbf{k}_N}^{\dagger} | 0 \rangle $$ Which has the normalization condition $$ \langle \mathbf{p}_1\mathbf{p}_2 \cdots \mathbf{p}_N | \mathbf{k}_1\mathbf{k}_2 \cdots \mathbf{k}_M \rangle = \delta_{NM} \sum_{\mathscr{P}} \prod_{j=1}^{N} \delta(\mathbf{k}_{j} - \mathbf{p}_{\mathscr{P}(j)}) $$ where $\mathscr{P}$ is the permutation over the integers $1,\ldots, N$. This is symmetric under interchange of any of the labels and 0 if $N \neq M$.

Answer 1

There is absolutely nothing wrong, you've simply abused the notation by setting different momenta exactly equal to each other. As a simple example, consider just one-particle states, where $$\langle \mathbf{p} | \mathbf{q} \rangle \sim \delta(\mathbf{p} - \mathbf{q}).$$ If you naively set the two momenta equal to each other, you would get $$\langle \mathbf{p} | \mathbf{p} \rangle \sim \delta(0).$$ Then you might complain this is not mathematically correct, because $\delta(0)$ is not a valid function or even a functional; it's not even defined at all. But the general overlap is perfectly well-defined as a functional of two variables, i.e. $\delta(\mathbf{p} - \mathbf{q})$. Your mistake was simply to set those variables equal, when you should have kept them explicit to integrate over later.

Similarly, for two-particle states the normalization is $$\langle \mathbf{p}_1 \mathbf{p}_2 | \mathbf{q}_1 \mathbf{q}_2 \rangle \sim \delta(\mathbf{p}_1 - \mathbf{q}_1) \delta(\mathbf{p}_2 - \mathbf{q}_2) \pm \delta(\mathbf{p}_1 - \mathbf{q}_2) \delta(\mathbf{p}_2 - \mathbf{q}_1)$$ with the sign depending on the particle statistics. The right-hand side is a perfectly well-defined functional of four variables, which you must integrate over later. However, if you naively start setting these variables equal to each other you'll get, for bosons, $$\langle \mathbf{p} \mathbf{p} | \mathbf{q} \mathbf{q} \rangle \sim \delta(\mathbf{p} - \mathbf{q})^2$$ which is not well-defined.