Introduction
I have taken the approach of using the Abel-Plana formula to evaluate the Lerch Phi as follows:
$$ \Phi(z,s,a) = \frac{1}{2 a^s} + \frac{(-\log(z))^{s-1}}{z^a} \Gamma(1-s, -a \log(z)) - 2 \int_0^{\infty} \frac{\sin(t \log z - s \operatorname{arctan}(t/a)}{(a^2 + t^2)^{s/2} (e^{2 \pi t}-1)} dt $$
According to Wikipedia, this definition is correct iff $\Re(a) > 0$. However, around the branch cut when trying to compute $\Phi(z,1-i,1+i)$, the results as given by the Abel-Plana formula and the results given by my copy of Mathematica 12 differ significantly: $\Phi(-5i, 1-i, 1+i)=-82.99760522684576 + 19.43284862059917i$ according to Abel-Plana, $\Phi(-5i, 1-i, 1+i)=-0.253694504765888 - 0.096558528415123i$ according to Mathematica 12. Curiously, this is not observed for $-4i$. At first I had thought that Wikipedia is once again mistaken and does not list the basics of the domain properly, but the branch cut seems so arbitrary that the $|z|<1$ condition doesn't seem to be the case here. Another speculation was that this happens due to choosing the wrong branch cut in upper incomplete gamma function. I have read about approximating the Lerch Phi using an integral representation with a Hankel contour, but performing contour integration is not particularly convenient in my scenario due to software limitations. The source seems to have noticed the problem with this formula as well.
My second thought was the following integral representation, which according to Wolfram Research converges for $Re(a)>0$, a condition I can easily satisfy: $$ \Phi (z,s,a)=\int_0^\infty (a+t)^{-s} z^t\ dt - 2 \int_0^\infty \frac{\sin\left(t \log\left(z\right)-s \operatorname{arctan}\left(\frac{t}{a}\right)\right)}{\left(a^2+t^2\right)^{0.5s} \left(e^{2\pi i} - 1\right)} dt + \frac{1}{2a^s} $$
However, this integral does not seem to converge under the specified conditions according to Mathematica 12:
Wikipedia specifies a further requirement of $|z|<1$ (so this times Wolfram Research is wrong), but this is not very helpful to me since I am still left without a way to approximate the Lerch Phi function due to problems with incomplete gamma function branch cuts.
For the $z,s,a \in \mathbb{R}$ case I use a method akin to the combined nonlinear-condensation transformation and the Van Wijngaarden transformation to accelerate alternating series as described in this paper (with a few corner-case issues in the original paper fixed), however, this method does not seem to generalise to the complex case.
Question
What technique should I use to approximate the Lerch Phi in the general case (regardless of the bounds of $\Re\ a$) that ideally does not involve contour integration? Or, how could I alter either of the integral formulae to converge to the correct value - perhaps by implementing an approximation to the incomplete gamma function that picks the desired branch cut somehow?
