With another approach we get
$(n, k) = (k+1, k)\;\implies\;\mathrm{E}[m_0] - \frac{kN}{n} = $ $$\frac{N}{k+1} - \frac{1}{(1+N)^k}\sum_{i=1}^N i^k$$ $(n, k) = (k+2, k)\;\implies\;\mathrm{E}[m_0] - \frac{kN}{n} = $ $$\frac{2N}{k+2} - \left(\frac{N+2}{(N+1)^k}-\frac{N+1}{(N+2)^k}\right)\sum_{i=1}^N i^k + \left(\frac{1}{(N+1)^k}-\frac{1}{(N+2)^k}\right) \sum_{i=1}^N i^{k+1} $$ where $m_0$ denotes your $m_{n-k+1}+\dots+m_n\,$. Applying Faulhaber's formula yields your five first results.
This is because $-$ ties put aside $-$ the numbers of seats $(m_0,\ldots,m_d)\in\mathbb{N}^{d+1}$ alloted to $d+1$ parties by the Jefferson-D'Hondt method is known to occur when the proportions of votes $(x_0, \dots, x_d)\in\mathbb{R}^{d+1}$ satisfy
$$\frac{m_i}{m_j+1}x_j \le x_i\qquad\forall i\ne j=0,\dots,d\;.$$
Choosing $x_d$ to stand for $1-x_0-\dots-x_{d-1}\,,\;$ these $d(d+1)$ inequalities upon $x_0,$ $\dots,$ $x_{d-1}$ intersect half-spaces of $\mathbb{R}^d$ and draw a rational bounded convex polytope $\mathcal{P} = \mathcal{P}(m_0,\dots,m_d)$ of $\mathbb{R}^d\,$.
The probability for $(m_0,\dots,m_d)$ to occur when the $x_i$'s follow a Dirichlet distribution with parameters $(\alpha_0,\dots,\alpha_d)$ is then $\mathrm{Pr}[m_0,\dots,m_d]=$
$$\frac{\Gamma(\alpha_0+\dots+\alpha_d)}{\Gamma(\alpha_0)\dots\Gamma(\alpha_d)}\int_{\mathcal{P}} x_0^{\alpha_0-1}\cdots x_{d-1}^{\alpha_{d-1}-1}(1-x_0-\dots-x_{d-1})^{\alpha_d-1} \,\mathrm{d}x_0\cdots \mathrm{d}x_{d-1}\;. $$
When $d=1$, writing $(x,y)$ or $(x,1-x)$ for $(x_0,x_1)\,$, and $(a,b)$ or $(a,N-a)$ for $(m_0,m_1)\,$, this is $$\frac{a}{b+1}(1-x)\le x\quad\text{and}\quad\frac{b}{a+1}x\le 1-x\;.$$ In other words $$x\in\mathcal{P} = \mathcal{P}(a,b) = \left[\frac{a}{N+1}, \frac{a+1}{N+1}\right]$$ so that $$ \mathrm{Pr}[a,b]= \frac{\Gamma(\alpha_0+\alpha_1)}{\Gamma(\alpha_0)\Gamma(\alpha_1)}\int_\frac{a}{N+1}^\frac{a+1}{N+1} x^{\alpha_0-1}(1-x)^{\alpha_1-1}\,\mathrm{d}x\;. $$
In the simpler version $n=1+k$ and $(\alpha_0,\alpha_1)=(k,1)$, we get $$ \mathrm{Pr}[a,b]= k\int_\frac{a}{N+1}^\frac{a+1}{N+1} x^{k-1}\,\mathrm{d}x = \left(\frac{a+1}{N+1}\right)^k - \left(\frac{a}{N+1}\right)^k $$ so that $$\begin{align*} \mathrm{E}[m_0]& = \sum_{a+b=N}a\,\mathrm{Pr}[a,b] = \frac{1}{(N+1)^k}\sum_{a=0}^N a\,\left((a+1)^k-a^k\right) = N - \frac{1}{(N+1)^k}\sum_{i=1}^N i^k \end{align*}$$
When $d=2$, indifferently writing $(x,y,z)$ or $(x,y,1-x-y)$ for $(x_0,x_1,x_2)\,$, and $(a,b,c)$ or $(a,b,N-a-b)$ for $(m_0,m_1,m_2)\,$, this is $$ \begin{array}{l} &\dfrac{a}{1+b}y\le x&\dfrac{a}{1+c}(1-x-y)\le x\\ \dfrac{b}{1+a}x\le y&&\dfrac{b}{1+c}(1-x-y)\le y\\ \dfrac{c}{1+a}x\le 1-x-y&\dfrac{c}{1+b}y\le 1-x-y& \end{array} $$
Here is a picture of the resulting $2$-dimensional $\mathcal{P}(a,b,c)$ polytopes when $N=a+b+c=6$:
$N=a+b+c=6$ seats" />
The vertices of $\mathcal{P}(a,b,c)$ are related to the following six binary strings and points of $\mathbb{R}^3$
$$\mathsf{100}\to\left(\tfrac{a+1}{N+1},\tfrac{b}{N+1},\tfrac{c}{N+1}\right),\\
\mathsf{101}\to\left(\tfrac{a+1}{N+2},\tfrac{b}{N+2},\tfrac{c+1}{N+2}\right),\qquad\mathsf{110}\to\left(\tfrac{a+1}{N+2},\tfrac{b+1}{N+2},\tfrac{c}{N+2}\right),\\
\mathsf{001}\to\left(\tfrac{a}{N+1},\tfrac{b}{N+1},\tfrac{c+1}{N+1}\right),\qquad\mathsf{010}\to\left(\tfrac{a}{N+1},\tfrac{b+1}{N+1},\tfrac{c}{N+1}\right),\\
\mathsf{011}\to\left(\tfrac{a}{N+2},\tfrac{b+1}{N+2},\tfrac{c+1}{N+2}\right).$$
They may not all be extreme points but they are always on the boundary. Edges are between any two points that happen to be at Hamming distance $1$.
In the simpler version $n=2+k$ and $(\alpha_0,\alpha_1,\alpha_2)=(k,1,1)$, the density is constant at $x$ constant and the length of the slice through the polytope $\mathcal{P}$ for $x$ constant is a piecewise linear function of $x$. Working out those linear terms yields the elementary integration $$\begin{align} \mathrm{Pr}[a,b,c] = k(k+1)\int_{\frac{a}{N+2}}^{\frac{a}{N+1}}x^{k-1}\left(\tfrac{N+2}{a}x-1\right)\mathrm{d}x + k(k+1)\int_{\frac{a}{N+1}}^{\frac{a+1}{N+2}}x^{k-1}\tfrac{1}{1+N-a}(1-x)\mathrm{d}x\qquad\\ + k(k+1)\int_{\frac{a+1}{N+2}}^{\frac{a+1}{N+1}}x^{k-1}\left(1-\tfrac{N+1}{a+1}x\right)\mathrm{d}x\\ = \left(\frac{a+1}{N+1}\right)^k - \frac{N-a}{N-a+1}\left(\frac{a+1}{N+2}\right)^k - \frac{N-a+2}{N-a+1}\left(\frac{a}{N+1}\right)^k + \left(\frac{a}{N+2}\right)^k \end{align}$$ which doesn't depend on $b$ or $c$ individually, just like the density of $(x,y,z)$ doesn't depend on $y$ or $z$ individually. For $N$ fixed we then get $$E[m_0] = \sum_{a+b+c=N} a\,\mathrm{Pr}[a,b,c] = \sum_{a=0}^Na\,(N-a+1)\,\mathrm{Pr}[a,0,N-a]$$ which yields the formula given in introduction.
I didn't figure out $d=3$ yet while experimenting with polymake. It is hard to come up with a description without many subcases, such as whether each $m_i$ is smaller or greater than $\frac{N+1}{2}$.