I finally did it!
What follows is the process that leads to the final result in terms of hyperbolic functions. Plus, two additional very interesting side-results.
I will start from the formula in $(\text{*})$: for $z=\frac{1+i}{2}$ the following holds
$$\sum_{n=1}^\infty \left( \frac{1}{(k+n)^2+n^2}-\frac{1}{(k+1-n)^2+n^2} \right)=$$
$$=\frac{i}{2k}\left(\psi\left(1+\bar{z}k\right)-\psi\left(1+zk\right)\right)+\frac{i}{2(k+1)}\left(\psi(z-\bar{z}k)-\psi(\bar{z}-zk)\right)$$
and then work my way through the cancelation of the digamma functions.
Tools:
The main tools we are going to use are the following:
$$\psi(1+u)=\psi(u)+\frac{1}{u} \tag{1}$$
$$\psi(1-u)-\psi(u)=\pi\cot(\pi u) \tag{2} $$
$$\psi(u+k)=\psi(u)+\sum_{n=0}^{k-1}\frac{1}{u+n}\tag{3}$$
$$\psi(-u)=\psi(u)+\frac1u+\pi\cot(\pi u)\tag{4}$$
Now, for reference purpose, let's name
$$\color{red}{S_1}=\frac{i}{2k}\left(\psi\left(1+\bar{z}k\right)-\psi\left(1+zk\right)\right) \hspace{1cm} \color{blue}{S_2}=\frac{i}{2(k+1)}\left(\psi(z-\bar{z}k)-\psi(\bar{z}-zk)\right)$$
and let's start with $S_1$.
PART $1$: $S_1$
It all starts by noticing that since $z=\frac{1+i}{2}$, we have that $\bar{z}=1-z$, therefore:
\begin{align}
\psi(1+\bar{z}k)-\psi(1+zk)&=\psi(1+k-zk)-\psi(1+zk) \\
&\stackrel{\text{(1)}}{=}\psi(1+k-zk)-\psi(zk)-\frac{1}{zk} \\
&\stackrel{\text{(3)}}{=}\psi(1-zk)+\sum_{n=0}^{k-1}\frac{1}{1-zk+n}-\psi(zk)-\frac{1}{zk} \\
&\stackrel{\text{(2)}}{=}\pi\cot(\pi zk)+\sum_{n=0}^{k-1}\frac{1}{1-zk+n}-\frac{1}{zk} \\
\end{align}
Next, plug $z$ in and simplify:
$$\frac{1}{zk}=\frac1k-\frac{i}{k}$$
$$\frac{1}{n+1-zk}=\frac{2}{2n+2-(1+i)k}=\frac{4n+4-2k}{(2n+2-k)^2+k^2}+i\frac{2k}{(2n+2-k)^2+k^2}$$
$$\cot(\pi zk)=\cot\left(\frac{\pi}{2}(1+i)k\right)=-\frac{\sin(k\pi)}{\cos(k\pi)-\cosh(k\pi)}+i\frac{\sinh(k\pi)}{\cos(k\pi)-\cosh(k\pi)}$$
and finally, putting these values back in $S_1$, together with the fact that $\sin(k\pi)=0$ and $\cos(k\pi)=(-1)^k$, we get
\begin{align}
\color{red}{S_1}&=\frac{i}{2k}(\psi(1+\bar{z}k)-\psi(1+zk)) \\
&=-\frac{\pi}{2k}\frac{\sinh(k\pi)}{(-1)^k-\cosh(k\pi)}-\frac{1+i}{2k^2}+\frac{i}{k}\sum_{n=0}^{k-1}\frac{2n+2-k}{(2n+2-k)^2+k^2}-\sum_{n=0}^{k-1}\frac{1}{(2n+2-k)^2+k^2} \\
\end{align}
And now, since the original series was always real, we can extract the real part and leave the result unchanged! Therefore
$$\boxed{\color{red}{S_1}=-\frac{1}{2k^2}-\frac{\pi}{2k}\frac{\sinh(k\pi)}{(-1)^k-\cosh(k\pi)}-\sum_{n=0}^{k-1}\frac{1}{(2n+2-k)^2+k^2}}$$
Additionally, since the result is real, then the imaginary part is $0$, and this leads to the first side-result:
$$\boxed{\sum_{n=0}^{k-1}\frac{2n+2-k}{(2n+2-k)^2+k^2}=\frac{1}{2k}}$$
Now that we are done with $S_1$, let's attack $S_2$.
PART $2$: $S_2$
Remember that
$$\color{blue}{S_2}=\frac{i}{2(k+1)}(\psi(z-\bar{z}k)-\psi(\bar{z}-zk))$$
As before, we will start by using the fact that $\bar{z}=1-z$:
\begin{align}
\psi(z-\bar{z}k) & =\psi(z-k+zk) \\
& \stackrel{\text{(4)}}{=}\psi(-z+k-zk)+\frac{1}{k-z-zk}+\pi\cot(\pi(k-z-zk)) \\
& \stackrel{\text{(3)}}{=}\psi(-z-zk)+\sum_{n=0}^{k-1}\frac{1}{n-z-zk}+\frac{1}{k-z-zk}+\pi\cot(\pi(k-z-zk)) \\
\end{align}
while for the second digamma:
\begin{align}
\psi(\bar{z}-kz) & =\psi(1-z-zk) \\
& \stackrel{\text{(1)}}{=}\psi(-z-zk)-\frac{1}{z+zk} \\
\end{align}
and therefore
$$\psi(z-\bar{z}k)-\psi(\bar{z}-zk)=\frac{1}{k-z-zk}+\frac{1}{z+zk}+\sum_{n=0}^{k-1}\frac{1}{n-z-zk}+\pi\cot(\pi(k-z-zk))$$
Again, plug $z$ in and simplify:
$$\frac{1}{k-z-zk}=\frac{2}{k-ik-(1+i)}=\frac{k-1}{k^2+1}+i\frac{k+1}{k^2+1}$$
$$\frac{1}{z+zk}=\frac{1}{1+k}-i\frac{1}{1+k}$$
$$\frac{1}{n-z-zk}=-\frac{k-2n+1}{n^2+(n-1)^2-2kn+k^2+2k}+i\frac{k+1}{n^2+(n-1)^2-2kn+k^2+2k} $$
$$\cot(\pi(k-z-zk))=-\frac{\sin(k\pi)}{\cos(k\pi)-\cosh((k+1)\pi)}+i\frac{\sinh((k+1)\pi)}{\cos(k\pi)-\cosh((k+1)\pi)}$$
and now, putting these values back in $S_2$, we end up with
\begin{align}
\color{blue}{S_2} & =\frac{i}{2(k+1)}(\psi(z-\bar{z}k)-\psi(\bar{z}-zk)) \\
& =\frac{i}{2}\frac{k-1}{(k+1)(k^2+1)}-\frac{1}{2(k^2+1)}+\frac{1+i}{2(k+1)^2}-\frac{i}{2(k+1)}\sum_{n=0}^{k-1}\frac{k-2n+1}{(k-n+1)^2+n^2}-\frac12\sum_{n=0}^{k-1}\frac{1}{(k-n+1)^2+n^2}-\frac{\pi}{2(k+1)}\frac{\sinh((k+1)\pi)}{(-1)^k-\cosh((k+1)\pi)}
\end{align}
but again, $S_2$ is real! So extracting the real part leads us to the final value for $S_2$:
$$\boxed{\color{blue}{S_2}=\frac{1}{2(k+1)^2}-\frac{1}{2(k^2+1)}-\frac12\sum_{n=0}^{k-1}\frac{1}{(k-n+1)^2+n^2}-\frac{\pi}{2(k+1)}\frac{\sinh((k+1)\pi)}{(-1)^k-\cosh((k+1)\pi)}}$$
And, for the second side-result, just take the imaginary part and set it equal to $0$:
$$\boxed{\sum_{n=0}^{k-1}\frac{k-2n+1}{(k-n+1)^2+n^2}=\frac{1}{k+1}+\frac{k-1}{k^2+1}}$$
Conclusion:
Finally, remember that we were after $\color{red}{S_1}+\color{blue}{S_2}$, and notice that
\begin{align}
\frac{\sinh(x)}{(-1)^k-\cosh(x)} & =-\coth\left(\frac{x}{2}\right) \hspace{1cm}\text{k even} \\
& =-\tanh\left(\frac{x}{2}\right) \hspace{1cm}\text{k odd} \\
\end{align}
which shows why in the conjectured form there are only $\tanh$ and $\coth$.
Furthermore, if I haven't made any typo, we now have an expression for the ratio $\frac{a_k}{b_k}$:
\begin{align}
\color{green}{\frac{a_k}{b_k}} & =\frac{1}{2(k+1)^2}-\frac{1}{2(k^2+1)}-\frac{1}{2k^2}-\frac12\sum_{n=0}^{k-1}\frac{1}{(k-n+1)^2+n^2}-\sum_{n=0}^{k-1}\frac{1}{(2n+2-k)^2+k^2} \\
\end{align}
and I have checked that this agrees with the first few cases. This ugly representation also explains why the two sequences $a_k$, $b_k$ were so hard to guess. It remains to check if this is consistent with @Ydd's result for $\frac{a_k}{b_k}$ as $k\to\infty$, but honestly I can't see how to do it now.
Out of all of this, what I find most intriguing are the two closed forms for the sums obtained setting $\Im=0$.
//Expandon your Mathematica code. – JimB Jul 23 '23 at 21:23