An analytic approximation for large $k$ and $a \ge 2$ is
$$ x=\frac{-k-a\cos(\pi/a)\sin(\pi/a) + a^2\,k\,\sin^2(\pi/a)}
{\sin^2(\pi/a) + a\,k\,\cos(\pi/a)\,\sin(\pi/a)} .$$
What's nice about it is that it reduces to $x=3k$ as given in Claude's initial exposition for $a=2.$ This is a large $k$ expansion. For a numerical check I let $a$ run from 2 to 3 by 0.05, and $k=5+1.5^m$ for $m$ = 0 to 20. The largest relative error of about 1.2% was obtained on the 'small k' side of $k=6.$
The technique was to take the tangent of eq. (1) above and write it, using the tangent addition formula, in the equivalent form
$$x=\tan(a\,\tan^{-1}(x/a)) - k\big(1+x\tan(a\,\tan^{-1}(x/a))\big). $$
At first I tried expanding $\tan(a\,\tan^{-1}(x/a))$ in a Laurent series around its singular point to get an equation in $x$ that could be solved analytically. Unfortunately it takes too many terms to give an accurate representation of the RHS of the equation. They will have to be included to go beyond the given approximation. In making the plots I noticed that a ballpark approximation to the solution was $x=a \tan(\pi/a).$
Furthermore the curve is decently approximated by a line over a fairly large region. Finally, a $k \to \infty$ analysis suggests that the point $x=a \tan(\pi/a)$ is the best. Thus the RHS of the previous equation is expanded around this point to first order in $(x-a \tan(\pi/a))$ and, solve for $x$ to get the equation
at the top of this answer.
In an article by Ejsmont and Lehner a method using operator trace is applied to evaluate finite sums involving the cotangent. The paper gives numerous references concerning finite trigonometric sums which seem to have a rich recent history. We use here a result obtained by Berndt and Yeap by a contour integration technique. It reads
\begin{equation}
\sum_{n=1}^{q-1}\cot^{2k}\frac{n\pi}{q}=(-1)^kq-(-1)^k2^{2k}\sum_{\substack{j_0,j_1,\ldots,j_{2k}\ge0 \\j_0+j_1+\ldots+j_{2k}=k}}q^{2j_0}\prod_{p=0}^{2k}\frac{B_{2j_p}}{(2j_p)!}\tag{BY}\label{BY}
\end{equation}
where $B_s$ are the Bernoulli numbers. It seems that no similar results with odd powers of the function were found.
To evaluate
\begin{equation}
S_{2k}=\sum _{n=1}^m \left[\cot \left(\frac{n \,\pi }{2 m+1}\right)\right]^{2k}
\end{equation}
we remark that, by symmetry,
\begin{equation}
S_{2k}=\frac{1}{2}\sum _{n=1}^{2m} \left[\cot \left(\frac{n \,\pi }{2 m+1}\right)\right]^{2k}
\end{equation}
then, by (\ref{BY}),
\begin{equation}
(-1)^kS_{2k}=m+\frac{1}{2}-2^{2k-1}\sum_{\substack{j_0,j_1,\ldots,j_{2k}\ge0 \\j_0+j_1+\ldots+j_{2k}=k}}(2m+1)^{2j_0}\prod_{p=0}^{2k}\frac{B_{2j_p}}{(2j_p)!}
\end{equation}
Under this form, it appears that the sum is polynomial in $m$.
The summation is however difficult to perform. It can be transformed by the use of the series expansion of $\coth Z$
\begin{equation}
Z\coth Z=1+\sum_{s=1}^\infty\frac{B_{2s}}{(2s)!}(2Z)^{2s}
\end{equation}
One has then
\begin{align}
S_{2k}&=(-1)^k\left( m+\frac{1}{2}\right)\left[ 1-\left[Z^{-1}\right]\left( \coth\left( (2m+1)Z \right)\coth^{2k}(Z) \right) \right]\\
&=(-1)^k\left( m+\frac{1}{2}\right)\left[ 1-\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \coth\left( (2m+1)z \right)\coth^{2k}(z)\,dz \right]\label{eq:Cauchy}
\end{align}
where $\left[z^r\right]\left( F(z) \right)$ represents the coefficient of $z^r$ in the Laurent series of $F(z)$ at $z=0$.
The latter expression for $S_{2k}$ vanishes when $m=0$, as $\operatorname{Res}\left( \coth^{2k+1}z,z=0 \right)=1$. Moreover, from the identity $\coth 2z=(\coth z+(\coth z)^{-1})/2$, it also vanishes for $m=1/2$ as proposed. Then, polynomials $Q_k(m)$ exist such that
\begin{equation}
S_{2k}=m(2m-1)Q_k(m)
\end{equation}
where
\begin{align}
Q_k(m)&=\frac{(-1)^k(2m+1)}{2m(2m-1)}\left[ 1-\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \coth\left( (2m+1)z \right)\coth^{2k}(z)\,dz \right]\\
&=\frac{(-1)^k(2m+1)}{2m(2m-1)}\left[ 1-\left[Z^{-1}\right]\left( \coth\left( (2m+1)Z \right)\coth^{2k}(Z) \right) \right]
\end{align}
Evaluated with a CAS, the latter formula gives the results which are given in the OP up to $k=5$. Higher values of $k$ were numerically tested up to $S_{20}$ for several values of $m$.
Some values of $m$ are remarkable. For instance, it is not difficult to show that $Q_k(-1)=(-1)^{k+1}/3$ and $Q_k\left( -\frac{3}{2} \right)=\frac{(-1)^{k+1}}{3}$. Moreover $Q_k\left( \frac{3}{2} \right)=\frac{1}{3}$. Indeed, for $m=3/2$,
\begin{align}
\coth\left( (2m+1)z \right)&=\coth(4z)\\
&=\frac{1}{2}\left( \coth 2z+\frac{1}{\coth 2z} \right)\\
&=\frac{1}{4}\left( \coth z+\frac{1}{\coth z} \right)+\frac{\coth z}{1+\coth^2z}
\end{align}
the integral reads
\begin{align}
\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \coth\left( 4z \right)\coth^{2k}(z)\,dz &=\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \frac14\left( \coth^{2k+1}(z)+\coth^{2k-1}(z) \right)\,dz+I_k\\
&=\frac{1}{2}+I_k
\end{align}
where
\begin{equation}
I_k=\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \frac{\coth^{2k+1}(z)}{1+\coth^2z}\,dz
\end{equation}
We have the recurrence relation
\begin{equation}
I_k=1-I_{k-1}
\end{equation}
and
\begin{align}
I_0&=\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon} \frac{\coth(z)}{1+\coth^2z}\,dz\\
&=\frac{1}{2}\frac{1}{2i\pi}\int_{\left|z\right|=\varepsilon}\tanh(2z)\,dz\\
&=0
\end{align}
Thus $I_{2s}=0$ and $I_{2s+1}=1$. Finally,
\begin{equation}
Q_k\left( \frac32 \right)=
\begin{cases}
\frac23\left( 1-\frac{1}{2}-0\right)=\frac{1}{3}&\text{ when } k=2s\\
-\frac23\left( 1-\frac{1}{2}-1 \right)=\frac{1}{3}&\text{ when } k=2s+1
\end{cases}
\end{equation}
Thus $Q_k\left( \frac{3}{2} \right)=\frac{1}{3}$.
This also shows, by adapting the signs, that, for $m=-5/2$, the values of the polynomials oscillate between $1/3$ and $-1/5$. Adapting to the decomposition of $\coth 3z$, should also make it possible to show (as suggested by @Blue in a comment) that $Q_k(1)=3^{-k}$, which is numerically observed.
Best Answer
If instead of using $a_m = \mathrm{lcm} \{ 1, 3, \dotsc, 2m-1 \}$ you use $b_m = (2m-1)!!$, you can write $$x=q-\sum _{m=1}^{\infty }\frac{Q_{m}(k) }{b_m } \frac 1 {(kq)^{2 m-1}}$$ and then $$Q_{m+1}(k) = \sum_{i=0}^m (-1)^{m-i} C(m, i) k^i$$ where the coefficients $C(m, i)$ satisfy $C(0, i) = (2i-1)!!$ and $$C(m, i) = \frac {(2m-1)2m(2m+1)} {(2m-i)(2m-i+1)} [C(m-1, i-1) + C(m-1, i)]$$
For a reference, see L. Comtet, Advanced Combinatorics, Reidel, 1974, p. 170.