If a (suitably normalised) holomorphic cusp newform has q-expansion $$f(z) = \sum_n \lambda_f(n) e(nz),$$

then we know the Hecke relations for $(mn,q)=1$, $$(\star) \qquad \lambda_f(m)\lambda_f(n) = \sum_{d | (m, n)} \lambda_f\left( \frac{mn}{d^2} \right).$$

In the case of a Maass cusp form, we can write the expansion $$u(z) = \sum_{n \geqslant 1} \rho_u(n) W_s(nz),$$

where $s$ is related in an explicit way to the associated eigenvalue and $W_s$ is a Whittaker function. Similarly, the continuous part of the spectrum made of Eisenstein series at a cusp $a$ admits expansion of the form $$E_a(s,z) = \phi_a y^s + \phi_a(s) y^{1-s} + \sum_{n \geqslant 1} \phi_a(n,s)W_s(nz)$$

I am interested in the following question :

Are there analogous "Hecke relations" for the coefficients of Maass forms, $\rho_u(n)$ and $\phi_a(n,s)$?

I suppose so, but I do not have any good reference for these matters.

The subconvexity problem for Artin $L$-functions(Inventiones 149, 2002), Section 6. In particular, the Hecke relations (or multiplicative commuting operators) are valid for any functions of period one. Since these operators commute with the Laplacian, Bob's your uncle and you can decompose into simultaneous eigenfunctions. Is this not also in one of Iwaniec's books? doi.org/10.1007/s002220200223 $\endgroup$Spectral Methods of Automorphic Forms(AMS Grad Studies in Math 53, 2002). doi.org/10.1090/gsm/053 $\endgroup$