Find Harmonic Numbers for Imaginary and Complex Values The Normal definition of Harmonic numbers with $ n \in \mathbb{N} $ is
$$ H_n = \sum_{k=1}^{n}\frac{1}{k} \tag{1}\label{eq1A} $$
This can be expanded to  $ n \in \mathbb{C} $
$$ H_n = \psi_0(n+1) + \gamma \tag{2}\label{eq2A}$$
Where $\psi_0(n)$ is the $0$th degree Polygamma function, which is defined for complex values of n, and $\gamma$ is the Euler-Mascheroni constant.
So my question is there a general solution to $H_{ji}$ where $j \in \mathbb{N}  $ and $i$ is the imaginary unit?
I considered in using the series formula of the Polygamma which is 
$$ \psi_0(z+1) = -\gamma + \sum_{n=1}^\infty(\frac{1}{n}-\frac{1}{n+z}) \tag{3}\label{eq3A}$$
Simplify the $(2)$ with $(3)$ we get
$$ H_{ji} = \sum_{n=1}^\infty(\frac{1}{n}-\frac{1}{n+ji}) \tag{4}\label{eq4A}$$
But I don't know how to simplify this series.
I tried to use the function above and Wolfram Alpha just simplifies it back to $\psi_0(x+1) + \gamma $ so this type of method seems to just be dead end.
Another method I have considered is the integral representation of $\psi_0(z)$ which is
$$ \psi_0(z) = \int_0^\infty (e^{-t}-\frac{1}{(1+t)^z})\frac{dt}{t} \tag{5}\label{eq5A}$$
Which transforms our $(2)$ into
$$ H_{ji} = \int_0^\infty (e^{-t}-{(1+t)^{-ji+1}})\frac{dt}{t} + \gamma \tag{6}\label{eq6A}$$
Expanding the integral int $(6)$ we get
$$ H_{ji} = \int_0^\infty \frac{e^{-t}}{t}dt - \int_0^\infty \frac{(1+t)^{-ji+1}}{t}dt + \gamma \tag{7}\label{e71A}$$
Which is just
$$ H_{ji} = \Gamma(0) - \int_0^\infty \frac{(1+t)^{-ji+1}}{t}dt +\gamma \tag{8}\label{eq8A}$$
Where $ \Gamma(z) $ is the Gamma function.
Seeing that $ \lim_{z\to 0} \Gamma(z) \rightarrow \infty $ and that the second integral doesn't converge there has to be some type of manipulation for $\Gamma(0)$ and the second integral to get a value of $H_{ji}$.
You can also use the identity of $ \psi(z+1) = \psi(z)+\frac{1}{z} $ to obtain for the previous function as
$$H_{ji} = \Gamma(0) - \int_0^\infty \frac{1}{(1+t)^{ji}t}dt - \frac{i}{j} +\gamma\tag{9}\label{eq9A}$$
Letting $ 1+t = u $ we can see our $(9)$ will change to
$$H_{ji} = \Gamma(0) - \int_1^\infty \frac{1}{u^{ji}(u-1)}du - \frac{i}{j} +\gamma \tag{10}\label{eq10A}$$
Doing partial fraction decomposition of $ \frac{1}{u^{ji}(u-1)} = \frac{A}{u^{ji}}+\frac{B}{u-1} $ we see that
$$ A = -1 \\
 B = 1^{1-ji} \tag{11}\label{eq11A}$$
Expanding our $(10)$ with $(11)$ we get
$$H_{ji} = \Gamma(0) - \int_1^\infty (\frac{1^{1-ji}}{u-1} - \frac{1}{u^{ji}})du - \frac{i}{j} +\gamma \tag{12}\label{eq12A}$$
Constructing $(12)$ into two integrals and simpifying we see find
$$ H_{ji} = \Gamma(0) - \int_1^\infty \frac{1}{u-1} du - \int_1^\infty \frac{1}{u^{ji}} du - \frac{i}{j} + \gamma \tag{13}\label{eq13A}$$
But it is obvious that these two integrals don't converge to any value, so it seems that the partial fraction decomposition is also a dead end.
Thanks to @AliShather, they noticed that the integral in $(8)$ is very closly related to the Beta function, where the beta function is
$$ B(x,y) = \int_0^\infty \frac{t^{x-1}}{(1+t)^{x+y}}dt = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} \Re\{x,y\}>0 \tag{14}\label{eq14A}$$
Comparing this to the integral in $(8)$ we can see that $\Re\{x,y\} \ngtr 0$, but very close!
Is there any better way to solve this?
Thank you for your time and patience!
 A: Let us consider the simplest case with an imaginary argument of the harmonic number, namely $H_i$. This can easily be generalized to the requested case $H_{ji}$.
The harmonic number can be defined for complex $z$ as
$$H_z = \sum_{k=1}^\infty \left( \frac{1}{k} - \frac{1}{k+z} \right)\tag{1}$$
Notice that you must not split the sum into two parts, because both sums
$$\sum_{k=1}^\infty \left( \frac{1}{k} \right),\sum_{k=1}^\infty \left( \frac{1}{k+z} \right) $$
are divergent.
Now for $z=i$ we have
$$H_i = \sum_{k=1}^\infty \left( \frac{1}{k} - \frac{1}{k+i} \right)\tag{2}$$
writing the summand as 
$$ \frac{1}{k} - \frac{1}{k+i} =  \frac{1}{k} - \frac{k-i}{(k+i)(k-i)}
= \frac{1}{k} - \frac{k-i}{k^2+1}=\frac{1}{k} - \frac{k}{k^2+1}+i \frac{1}{k^2+1}\\=\frac{1}{k(k^2+1)} +i \frac{1}{k^2+1}
$$
we get
$$H_i = \sum_{k=1}^\infty \left(\frac{1}{k(k^2+1)} +i \frac{1}{k^2+1}\right)\tag{3}$$
Now splitting is permitted because the two sums are convergent. In fact, the real and imaginary parts of $H_i$, $f$ and $g$, respectively, are
$$f = \sum_{k=1}^\infty \left(\frac{1}{k(k^2+1)}\right)\tag{4}$$
$$g = \sum_{k=1}^\infty \left( \frac{1}{k^2+1}\right)\tag{5}$$
I stop here for a while to let you calculate $f$ and $g$, i.e. find expressions which are not just real and imaginary part of $H_i$.
EDIT
Now, what can be said about $f$ and $g$? 
$g$ has a closed form
$$g = \sum_{k=1}^\infty \left( \frac{1}{k^2+1}\right)= \frac{1}{2}\left(\sum_{k=-\infty}^\infty \left( \frac{1}{k^2+1}\right)-1 
 \right) \\= \frac{1}{2} (\pi  \coth (\pi )-1) = \dfrac{\pi-1}{2}+\dfrac{\pi}{e^{2\pi}-1}\tag{6}$$
which has been given in several places in this forum, for instance here How to prove $\sum_{n=0}^{\infty} \frac{1}{1+n^2} = \frac{\pi+1}{2}+\frac{\pi}{e^{2\pi}-1}$ (notice that there the sum starts at $k=0$) and derived with complex contour integration here How to sum $\sum_{n=1}^{\infty} \frac{1}{n^2 + a^2}$?.
For $f$ I have not found a closed expression other than the information that $f$ is the real part of $H_i$. However, normally this would be considered a closed form as well.
Technically, the deeper reason for the different behaviour of $f$ and $g$ is that whereas $g$ can be written as a symmetric sum from $-\infty$ to $\infty$ which allows complex contour integration with a kernel $\pi \cot(\pi z)$, $f$ is a one-sided sum which has the kernel $H_{-z}$. The latter kernel then just brings us back to where we came from. Usage of contour integrals for infinite sums is described for example in chapter 2 of  https://projecteuclid.org/download/pdf_1/euclid.em/1047674270
A: In case the OP is curious about how to evaluate $f$ and $g$ in  Dr. Wolfgang Hintze's solution above,

For $f$, we use the generalization $$H_a=\sum_{k=1}^\infty\frac{a}{k(k+a)}$$
so
$$ H_a+H_{-a}=\sum_{k=1}^\infty\frac{a}{k(k+a)}+\sum_{k=1}^\infty\frac{-a}{k(k-a)}=\sum_{k=1}^\infty\frac{-2a^2}{k(k^2-a^2)}$$
set $a=i$ we get
$$ H_i+H_{-i}=\sum_{k=1}^\infty\frac{2}{k(k^2+1)}$$
for $g$, we use the classical generalization 
$$\sum_{k=1}^\infty\frac1{k^2+a^2}=\frac1{2a^2}\left(a\pi\coth(a\pi)-1\right)$$
$$\Longrightarrow g=\sum_{k=1}^\infty\frac{1}{k^2+1}=\frac12(\pi\coth(\pi)-1)$$
so
$$H_i=f+ig=\frac12H_i+\frac12H_{-i}+\frac i2(\pi\coth(\pi)-1)$$
or
$$H_i-H_{-i}=i(\pi\coth(\pi)-1)\tag1$$

By the way, we can quickly reach the result in $(1)$ if we use the identity
$$H_{-a}-H_{a-1}=\pi\cot(a\pi)$$
$$\Longrightarrow H_{-i}-H_{i-1}=-i\pi\coth(\pi)$$
plug $H_{i-1}=H_i-\frac1i=H_i+i$
we get
$$H_{-i}-H_i=-i(\pi\coth(\pi)-1)$$

Its obvious that $$\Re\{H_{-i}\}=\Re\{H_{i}\}$$
and
$$\Im\{H_{-i}\}=-\Im\{H_{i}\}$$
so $$\Im\{H_{i}\}-\Im\{H_{-i}\}=2\Im\{H_{i}\}=(\pi\coth(\pi)-1)$$
or
$$\Im\{H_{i}\}=\frac{\pi\coth(\pi)-1}{2}$$
or we know from Dr. Wolfgang Hintze's solution that $g=\Im\{H_{i}\}$ 

It would be interesting if we can find 
$$H_{i}+H_{-i}\tag2$$
because from $(1)$ and $(2)$ we can find $H_{i}$ and $H_{-i}$.
