Can we not proceed in a manner similar to how Riemann did with the Mellin transform?
For example, let me define $\theta(u) = \sum_{n=0}^{\infty}e^{-\pi a_n^2u}$
$$\int_0^{\infty} \theta(y) y^{s/2} dy/y = \int_0^{\infty}\sum_{n=0}^{\infty}e^{-\pi a_n^2y} = \sum_{n=0}^{\infty} \int_0^{\infty}e^{-\pi a_n^2y} y^{s/2} dy/y $$
$$= \pi^{-s/2} \sum_{n \ge 0} \frac{1}{a_n^s} \int_0^{\infty}e^{-u} u^{s/2} du/u$$ (after substituting $\pi n^2 y = u$)
$$= \pi^{-s/2} \Gamma(s/2) F(s)$$
I don't know if Jacobi's identity $\theta(u) = u^{-1/2}\theta(1/u)$ holds or not. It is probable that the identity holds, and then the functional equation will look similar to RFE and then one can probably go on to guess all zeros have real part $1/2$.