What is the limit of this Goldbach Conjecture related Dirichlet series? Background
Let us recollect the Goldbach's conjecture. For any $n \geq 2$ there exists primes $p_i$ and $p_j$ such that:
$$ 2n = p_i + p_j $$
where $p_k$ is the $k$'th prime. Now we define a function $G(2n)$ as the number of "unique" solutions for $2n$. For example:
$$ G(4) = 2 + 2 = p_1 + p_1$$
Since $p_i = p_j$ the number of "unique" solutions is $1$. However for $G(8)$:
$$ G(8) = 3 + 5 = p_2 + p_3$$
Since $p_i \neq p_j$ the number of solutions is $2$. What about sometime like $G(14)$?
$$G(14) = 7+ 7 = 3+ 11 $$
Hence the number of unique solutions is $G(14)=1 + 2 = 3$
Essentially for each prime pairs $p_i$ and $p_j$ whose sum is $2n$ if we add all such pairs (twice if $i \neq j$ and once if $i=j$). If there are no solutions for the number $k$ then $G(k) = 0$ (for example $G(2) = 0$).
Question
Is the below conjecture true or false?
$$ \frac{P(s)^2}{\zeta(s)} \sim 2 \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}} $$
As $s \nearrow 1$. Or:
$$ \lim_{s \to 1} \frac{\frac{P(s)^2}{\zeta(s)} }{ \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}}} = 2$$
Heuristic Reasoning
Let's think of the prime zeta function as Dirchlet series:
$$ P(s) = 2^{-s} + 3^{-s} + \dots $$
Choosing a function $f(x) = e^{-x}$ and this formula:
$$ (e^{-2\epsilon} + e^{-3 \epsilon} + e^{-5 \epsilon} + \dots) \epsilon \sim \frac{P(s)}{\zeta(s)} \times \int_0^\infty f(x) dx $$
as $s \nearrow 1$ and $\epsilon$ tends to $0$.
Removing $e^{-2 \epsilon}$ since it is small
$$ ( e^{-3 \epsilon} + e^{-5 \epsilon} + \dots) \epsilon \sim \frac{P(s)}{\zeta(s)} $$
since $\int_0^\infty f(x) dx = 1$. Squaring both sides and resuming the L.H.S:
$$ ( G(6)e^{-6 \epsilon} + G(8)e^{-8 \epsilon} + \dots) \epsilon^2 \sim (\frac{P(s)}{\zeta(s)})^2 $$
However recognizing the L.H.S as an integral of a function $g(x) = e^{- \sqrt x}$ in it's own right with $\epsilon^2 \to \epsilon' $:
$$ ( G(6)e^{-6 \sqrt \epsilon '} + G(8)e^{-8 \sqrt \epsilon'} + \dots) \epsilon' \sim \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}} \times \frac{1}{\zeta(s)} \times \int_0^\infty g(x) dx $$
Comparing the above $2$ equations:
$$ \frac{P(s)^2}{\zeta(s)} \sim 2 \sum_{r=3}^\infty \frac{G(2r)}{(2r)^{2s}} $$
since $\int_0^\infty g(x) dx = 2$
 A: It is clearly false as $(\frac{P(s)}{\zeta(s)})^2 \to 0$ as $s\to 1^+$ in contrary to your Dirichlet series with non-negative coefficients.
It follows from the PNT that as $x\to 0^+$
$$f(x)=\sum_{p\ge 3} e^{-px}=(1-e^{-x}) \sum_{n\ge 2} (\pi(n)-1)e^{-nx}$$ $$\sim
(1-e^{-x}) \sum_{n\ge 2} \sum_{2\le m\le n} \frac1{\log m} e^{-nx}=
 \sum_{n\ge 2} \frac{e^{-nx}}{\log n} $$
For $s > 2$
$$\sum_{r \ge 1} G(2r)(2r)^{-2s}=\frac1{\Gamma(s)}\int_0^\infty \sum_{r\ge 1} G(2r) e^{-2rx}x^{s-1}dx=\frac1{\Gamma(s)}\int_0^\infty f(x)^2 x^{s-1}dx$$
From this, since $f(x)=O(\frac{e^{-x}}{1-e^{-x}})$, you get that for some $0<A<B$
$$\sum_{r \ge 1} G(2r)(2r)^{-2s}\in (A,B)\int_0^{1/2} f(x)^2 x^{s-1}dx$$
As $x\to 0$, $$f(x)= O(\frac{1}{x\log x}),\qquad \frac{1}{x\log x}=O(f(x))$$ whence you get that for some $0<a<b$
$$\sum_{r \ge 1} G(2r)(2r)^{-2s} \in (a,b)\int_0^{1/2} \frac1{x^2 \log x^2}x^{s-1}dx$$ which converges for $s\ge 2$. The RHS diverges for $s<2$ and so does the LHS.
A: Letting $G_0(x)=\sum_{n=3}^{\infty}G(n)x^n$ and $L(G,s)=\sum_{n=1}^{\infty}\frac{G(2n)}{n^s}$, I suggest that you use the fact that $L(G,s)$ is the Mellin transform of $G_0(e^{-t})$, i.e
$$L(G,s)=\frac{1}{\Gamma(s)}\int_0^\infty x^{s-1}G_0(e^{-x})dx$$
Using the fact that $G_0(x)=P^2_0(x):=\left(\sum_{p\geq 3}x^p\right)^2$ you can get a strong handle on $L(G,s)$ as $s\sim1$, i.e
\begin{align*}
L(G,s)&=\frac{1}{\Gamma(s)}\int_0^\infty x^{s-1}P_0^2(e^{-x})dx\\
&\sim \int_0^1 x^{s-1}P_0^2(e^{-x})dx\tag{Since integral converges for $x\ge1$}\\
&\sim \int_0^1 x^{s-1}\left(\sum_{n}\frac{1}{\log(n)}x^n\right)^2dx \tag{Summation by parts+PNT}
\end{align*}
This leaves you with, in a sense, a totally elementary integral and you should be able to manipulate it to find an asymptotic relationship of your pleasing.
