The balanced polygamma function $\Psi^{(s)}(x)$ for $x=1$ can be expressed as:

$$\Psi^{(s)}(1)=\dfrac{\big(\Psi(-s)+\gamma\big)\,\zeta(s+1)+\zeta'(s+1)}{\Gamma(-s)}$$

Note $\Psi(s)$ is the digamma function and $\gamma$ the Euler–Mascheroni constant.

$\Psi^{(s)}(1)$ has complex zeros and it appears their real parts all reside in a strip within $[-1,-\frac12]$. Their imaginary parts seem correlated and get closer to $\Im(\rho)$ for increased $T$ ($\rho$ is a non-trivial zero of $\zeta(s)$).

However, numerical evidence suggests the complex zeros (real ones exist as well) of

$$\Psi^{(s)}(1) \pm \Psi^{(1-s)}(1)$$

all reside on the critical line $\Re(s)=\frac12$.

Only focusing on the critical strip $0 \lt\Re(s) \lt 1$ allows rewriting $\zeta(s+1)$ and $\zeta'(s+1)$ into:

$$\frac{1}{\Gamma(-s)} \sum_{n=1}^{\infty}\frac{\gamma +\Psi(-s)-\ln(n)}{n^{s+1}} \pm \frac{1}{\Gamma(s-1)} \sum_{n=1}^{\infty}\frac{\gamma +\Psi(s-1)-\ln(n)}{n^{2-s}}$$

Another approach could be to express $\zeta(s+1)$ and $\zeta(2-s)$ as their Euler-products.

Is there a way to simplify this further or to prove that all complex zeros are on the line $\Re(s)=\frac12$?

The graph below illustrates the claim for $\Psi^{(s)}(1) - \Psi^{(1-s)}(1)$. (Note: at a first glance the distribution of the zeros looks deceitfully regular; however, small distortions do seem to occur in the vicinity of the $\rho$s.)