Intuitively (from the viewpoint of Fourier Analysis),
I would expect the prime zeta function $P(s)=\zeta_P(s)$
to be undefined on the imaginary axis, as argued below.
Wikipedia says that "singularities cluster near all points of"
the line $\Re s=0$, which also bounds the critical strip
of the (full) Riemann zeta function on the left in the complex plane.
Perhaps you should call the partial sums (over $p\leq x$),
which you are treating, the "incomplete" or "truncated"
prime zeta function.
For $s=it$ with $0\neq t\in\mathbb{R}$ and $0\neq k\in\mathbb{Z}$,
$$
p^{-ks}=e^{-ikt\log p}=\cos(kt\log p)-i\sin(kt\log p)
$$
so that for a fixed $k$, the sum over primes $p$
$$
P(ks)
=\sum_{p\text{ prime}}p^{-ks}
=\sum_{p}e^{-ikt\log p}
=\sum_{p}e^{-ik\theta_p}
$$
should behave like a "random walk" with unit (or "quantum") steps
but whose directions $-k\theta_p$ are not discrete but, rather, continuous, "incoherent" or "decoherent" (in the spectral/polarization/quantum physical sense, in the sense of compressed sensing, where its significance connects with the concept of basis functions in Fourier analysis, or in the general sense of meaning uncorrelated), aperiodic and decreasing, since
$$
\theta_p=2\pi
\left(
\frac{t\log p}{2\pi}-\left\lfloor
\frac{t\log p}{2\pi}\right\rfloor
\right)
$$
where the rate of change of direction (modulo $2\pi$ we can ignore the jumps)
$
\frac{\partial}{\partial p}\theta_p=\frac{t}{p}
$
decreases in magnitude (deccelerates) slowly
(and hence the autocorrelation between terms increases)
for $p$ sufficiently large.
Since each $\log p$ is irrational (presumably also transcendental)
and each $\theta_p$ is probably also irrational
(except perhaps for one value of $p$, due to $t$),
the directions never repeat and are therefore dense modulo $2\pi$,
so that the sum or expected value should be undefined.
Compare this with a Fourier Series or the Jacobi Theta function.
I would not expect convergence from the sequence of partial sums,
unless some kind of symmetry were exploited (for example between $k$ & $m$)
or imposed in some manner, for example analogously to
Cauchy principle values for improper integrals,
or perhaps some other sort of regularization.
Because of the irrationality of $\log p$ for each $p$,
I also would not expect any such symmetry trick for $m,n$
as one sees for example when computing that
$\{e^{inx}\}_{n\in\mathbb{Z}}$ form an orthonormal basis
in Hilbert space.
If we were given a function $f\in L^2([-\pi,\pi])$
so that the sequence of Fourier coefficients
$\{c_k=\hat{f}\}$ in $\mathbb{C}$
tended toward zero as $k\rightarrow\pm\infty$,
then we could form a convergent Fourier series
and interchange the order of summations
in the double summation below to obtain
$$
\sum_{k\in\mathbb{Z}}c_kP(ks)=
\sum_{k\in\mathbb{Z}}
\sum_{p}c_ke^{-ik\theta_p}=
\sum_{p}
\sum_{k\in\mathbb{Z}}
c_ke^{-ik\theta_p}=
\sum_{p}f(-\theta_p),
$$
but this would not be allowed for $c_k\equiv1$
since then $f$ would not be integrable.
The truncated sums of $P(ks)$,
except for being deterministic (and the angles $\theta_p$ being aperiodic),
also bear some resemblance to quantum random walks
(see the discussion around equation 44 in the linked introductory paper).
A characteristic difference between classical and quantum random walks
is that the expected value $\mathbb{E}|X_n|$
is $O(\sqrt n)$ for classical, but $O(n)$ for quantum random walks,
i.e. QRWs are ballistic rather than diffusive.
In this context, it would be interesting to see
if your results show ballistic-diffusive "crossover".
