Let $D \equiv 1, 2 \, (\textrm{mod }4)$ be a positive, squarefree integer. Let \begin{align} r_D(n) = \{(x, y) \in \mathbf{Z}^2 \mid x^2 + Dy^2 = n\} \end{align} for any positive integer $n$. Consider the Dirichlet series \begin{align} R_D(s) := \sum_{n \geq 1} \frac{r_D(n)}{n^s} = \sum_{\substack{(x, y) \in \mathbf{Z}^2 \\ (x, y) \neq (0, 0)}} \frac{1}{|x + y\sqrt{-D}|^2} = D^{-s/2} \sum_{\substack{(x, y) \in \mathbf{Z}^2 \\ (x, y) \neq (0, 0)}} \frac{(\sqrt{D})^s}{|x + y\sqrt{-D}|^2}. \end{align} Therefore the values taken by $R_D(s)$ are (almost) special values of the Eisenstein series \begin{align} \tilde{G}_{2s}(\tau) = \sum_{\substack{(x, y) \in \mathbf{Z}^2 \\ (x, y) \neq (0, 0)}} \frac{\textrm{Im}(\tau)^s}{|x + y\tau|^{2s}} \end{align} obtained by taking $\tau = \sqrt{-D}$. The function $\tilde{G}_{2s}(\tau)$ is weakly modular of weight zero with a simple pole at $\infty$. Any such function is a rational function in $j$, where \begin{align} j(\tau) = \frac{1}{q} + 744 + 196884q + \dots \end{align} is the $j$-invariant function. Hence we can write \begin{align} R_D(s) = D^{-s/2}P_s\Big(j\big(\sqrt{-D}\big)\Big) \end{align} where $P_s(z) \in \mathbf{C}(z)$ are some rational functions, depending on $s$. I want to know what can be said about these $P_s$ functions as $s$ varies.
To motivate: if the coefficients of $P_s$ are rational functions of $s$, I immediately find that the sums \begin{align} \sum_{n \geq 1} \frac{r_D(n)}{n^k} \end{align} are algebraic numbers for all rational $k$ where the sum converges.
Note that the function $\tilde{G}_{2s}(\tau)$ only takes real values, so it cannot be meromorphic unless it is constant (it would be constant on an open set by Cauchy-Riemann, hence constant everywhere by identity theorem). Although it is weakly modular of weight zero, we cannot say that it is a rational function in $j$, i.e., it is not meromorphic so it is not in the function field of the curve $X(1)$.