Let $f(z) = \sum_{n=1}^\infty a(n) e^{2i \pi nz}$ be an eigenform of $S_k(\Gamma_0(N))$. Since the Hecke operator acts by $T_p f = a_p f$ the Riemann hypothesis for $f$'s L-function is $$ \!\!\! \!\!\! \!\!\! \!\!\!(\sum_{n=1}^\infty a_n n^{-s})^{-1} =\prod_p (1-a(p) p^{-s}+p^{k-1-2s}) =\frac{\displaystyle[\prod_p (1-T_p p^{-s}+p^{k-1-2s})] f(z)}{f(z)} \tag{1}$$ converges and is analytic for $\Re(s) > k/2$.
Then, the Riemann hypothesis for all the eigenforms of $S_k(\Gamma_0(N))$ is that for any $f \in S_k(\Gamma_0(N))$, $$[\prod_p (1-T_p p^{-s}+p^{k-1-2s})] f(z) \quad \text{ is analytic for } \Re(s) > k/2\tag{2} $$
Questions :
This suggests to define a Riemann hypothesis for the Hecke operators themselves, and I would like to know if there is a well-known way to think about that.
The RH for the Hecke operators could be the convergence in operator norm of $\lim_{x \to \infty} \prod_{p \le x} (1-T_p p^{-s}+p^{k-1-2s}),$ for $\Re(s) > k/2$, with the norm coming from the Petersson inner product.
But there are other possible norms, for example $\langle f,g\rangle = \int_0^\infty f(ix) \overline{g(ix)} x^{2 \sigma -1}dx$. Indeed, the choice of the norm is a major problem : when defining $\displaystyle T_p f(z) = p^{k-1}\sum_{ad =p, b \bmod p} d^{-k} f(\frac{az+b}{d})$$\displaystyle T_p f(z) = p^{k-1}\sum_{ad =p, b \bmod p} d^{-k} f(\frac{az+b}{d}) \tag{3}$ the statement in $(1)$ works the same way whenever $f(z) = \sum_{n=1}^\infty a_n e^{2i \pi n z}$ and $\sum_{n=1}^\infty a_n n^{-s} = \prod_p (1+a_p p^{-s}+ p^{k-1-2s})^{-1}$, no matter that $f$ is modular or not, so that in general $f$ doesn't have a Riemann hypothesis (nor the $T_p$ operator acting on it). So we really need a norm and a statement specific to modular forms, for example about $\prod_p (1-T_pp^{-s}+\langle p \rangle p^{k-1-2s})$ (which works for $S_k(\Gamma_1(N))$ too).
For short, RH is believed true for Dirichlet series with Euler product and functional equation. $\prod_p (1-T_p p^{-s}+p^{k-1-2s})$ is just an Euler product. So hopefully, we only need to add a reference to the modularity (implying the functional equation) to make a viable RH statement.
Assuming we solved that part, the Hecke operators depend on $N$ only for finitely many $p$, so can we expect the statement to imply the Riemann hypothesis for $\displaystyle\bigcup_N S_k(\Gamma_0(N))$ ? In that case, what's about the dependence on $k$ ?
Also looking at the weight-$\frac{1}{2}$ forms $\sum_{n \ge 1}^\infty \chi(n) e^{2i \pi n^2 z}$ could help.