Here is an attempt. Start from your integral expression.
Writing $\frac{1}{(x_1+\ldots+x_n)^{\alpha_1+\ldots+\alpha_n}}=\int_0^{\infty}\frac{ t^{\alpha_1+\ldots+\alpha_n-1}}{\Gamma(\alpha_1+\ldots+\alpha_n)}\,e^{-t(x_1+\ldots+x_n)}\,dt$ and changing the order of integration (and denoting by $p$ the distribution of $\tau(N)(P)$) leads to $$p(m_1,\ldots,m_n)=\sum_{j=1}^n \int_0^\infty \frac{m_j^{\alpha_j}t^{\alpha_j-1}}{\Gamma(\alpha_j)}\prod_{i\neq j}\big (R_{\alpha_i}(m_it)-R_{\alpha_i}((m_i+1)t)\big) e^{-m_jt}\,dt$$ where $R_{\alpha} (t):=\int_t^\infty \frac{y^{\alpha-1}}{\Gamma(\alpha)} e^{-y}\,dy$.
For $\alpha=k$ (a positive integer) $$R_k(t)=q_k(t) e^{-t} \mbox{ with } q_k(t)=\sum_{i=0}^{k-1}\frac{t^i}{i!}\;\;.$$ If all $\alpha_i$ are all integers the integrals above therefore lead to explicit rational expressions for the probabilities $p(m_1,\ldots,m_n)$. For the case $\alpha_1=\ldots=\alpha_{n-1}=1, \alpha_n=k$ (positive integer) one gets \begin{align*} p(m_1,\ldots,m_n)= &\frac{m_n^k}{(k-1)!}\int_0^\infty t^{k-1} e^{-Nt} (1-e^{-t})^{n-1}\,dt\\& + (N-m_n)\int_0^\infty \big(q_k(m_nt)-e^{-t}q_k((m_n+1)t)\big) e^{-Nt}(1-e^{-t})^{n-2}\,dt\end{align*} (and this can easily be worked out explicitly using $\Gamma$-integrals). Since $p(m_1,\ldots,m_n)$ depends on $m_1,\ldots,m_{n-1}$ only through $m_1+\ldots+m_{n-1}=N-m_n$ the distribution of the last coordinate $M_n$ is proportional to $p$. I assume (but haven't checked) that computer algebra can now find explicit expressions for the expectation of $M_n$.