Is it possible to justify these approximations about prime numbers? A recently closed question asked for a possible closed form of the infinite summation
$$f(a)=\sum _{i=1}^{\infty } a^{-p_i}$$ for which I already proposed a first simple but totally empirical approximation.
Since we quickly face very small numbers, I tried to find approximations of
$$g(a)=\Big[\sum _{i=1}^{\infty } a^{-p_i}\Big]^{-1} \qquad \text{and} \qquad h(a)=\Big[\sum _{i=1}^{\infty } (-1)^{i-1} a^{-p_i}\Big]^{-1}$$
All calculations where done with integer values of $a$ for the range $2 \leq a \leq 1000$.
What I obtained is
$$\color{blue}{g(a)\sim\frac{(a-1) (2a^3+2a-1)}{2 a^2}}\qquad \text{and} \qquad \color{blue}{h(a)\sim\frac{(a-1) \left(a^3+2 a^2+3 a+4\right)}{a^2}}$$
If the corresponding curve fits were done, in both cases we should have $R^2 > 0.999999999$.
For the investigated values of $a$,
$$\text{Round}\left[\frac{(a-1) (2a^3+2a-1)}{2 a^2}-{g(a)}\right]=0$$
$$\text{Round}\left[\frac{(a-1) \left(a^3+2 a^2+3 a+4\right)}{a^2}-{h(a)}\right]=0$$

Not being very used to work with prime numbers, is there any way to justify, even partly, these approximations ?


These estimates are correct within a reasonable degree of accuracy. Below is the explanation for $f(a)$; the case for $h(a)$ can be dealt similarly. We have
$$
f(a) = \frac{1}{a^2} + \frac{1}{a^3} + \frac{1}{a^5} + O\bigg(\frac{1}{a^7}\bigg)
$$
whereas
$$
\frac{2a^2}{(a-1)(2a^3 + 2a - 1)} = \frac{1}{a^2} + \frac{1}{a^3} + \frac{1}{2a^5} + \frac{3}{2a^6} + O\bigg(\frac{1}{a^7}\bigg).
$$
Hence,
$$
f(a) = \frac{2a^2}{(a-1)(2a^3 + 2a - 1)} + O\bigg(\frac{1}{a^5}\bigg)
$$
For large values of $a$ the error would obviously be negligible, since it grows no faster than a constant times $a^{-5}$. So this may or may not be a good estimate depending upon weather you are satisfied with the magnitude of the error term $O(a^{-5})$.
The best possible estimate of the form $\dfrac{Ax^2}{(x-1)(Bx^3 + Cx^2 + Dx + E)}$ is obtained by the Laurent series expansion of about the point $x = \infty$ and equating the coefficient of smallest non prime powers to zero which gives $A = B = D = 1, C = 0,E = -1$.
Hence we have,
$$
f(a) = \frac{a^2}{(a-1)(a^3 + a - 1)} + O\bigg(\frac{1}{a^6}\bigg)
$$
which reduces the error by a factor of $a$.
Update 21-Jul-2020: However, using basic properties of primes we can get remarkably sharper estimates. Since every primes $\ge 5$ are of the form $6k \pm 1$, by summing up the geometric sequences $a^{-6k-1} + a^{-6k+1}$ for $k = 1,2,\ldots, \infty$ and adding $a^{-2} + a^{-3}$, and taking advantage of the fact the the density of primes among the first few numbers of these form is high we get
$$
f(a) = \frac{a^7 + a^6 + a^4 + a^2 -a - 1}{a^3(a^6 - 1)} + O\bigg(\frac{1}{a^{25}}\bigg)
$$
