Let $f$ a weight $k$ holomorphic Hecke cusp form with $\|f\|^2=\langle f,f\rangle=1$ with fourier expansion $$f(z)=\sum_{r\geq 1}a_f(r)e(rz)$$ Let $\displaystyle\lambda_f(r)=\frac{a_f(r)r^{(-k+1)/2}}{a_f(1)}$, then $\lambda_f(r)$ are the eigenvalues of the normalized Hecke operators $n^{(1-k)/2}T_n$. We define $$L(s,\operatorname{sym}^2(f))=\zeta(2s)\sum_{n=1}^\infty \frac{\lambda_f(n^2)}{n^s}$$ I am trying to prove the following (this appears in different papers like this one in Luo-Sarnak)
$$|a_f(1)|^2=\frac{(4\pi)^{k-1}}{\Gamma(k)}\frac{2\pi^2}{L(1,\operatorname{sym}^2(f))}.$$
What I was able to prove so far is that $$a_f(1)=\frac{\zeta(2s)}{L(s,\operatorname{sym}^2 (f))}\frac{L(k+s-1,\operatorname{sym}^2 (f))}{\zeta(2(k+s-1))}$$ If I put $s=1$, the problem reduces to find $L(1,\operatorname{sym}^2(f))$ and $L(k,\operatorname{sym}^2(f))$, but I don't think this is a good approach. Any ideas or references would be appreciated.
Thanks to @Peter Humphries in telling me the approach in his comment. Denote $\Gamma=SL_2(\mathbb{Z})$ and consider $$F(s)=\int_{\Gamma\setminus\mathbb{H}}|f(z)|^2y^kE(z,s)d\mu(z)$$ where $E(z,s)=\sum_{\gamma\in\Gamma_\infty\setminus\Gamma}(\operatorname{Im}(\gamma z))^s$ is the usual Eisenstein series.
First, since $\operatorname{Res}_{s=1}(E(z,s))=\frac{3}{\pi}$, we have $$\operatorname{Res}_{s=1}F(s)=\frac{3}{\pi}\int_{\Gamma\setminus\mathbb{H}}|f(z)|^2y^kd\mu(z)=\frac{3}{\pi}\langle f,f\rangle=\frac{3}{\pi}$$ Now, let's compute this residue in a different manner:
By modularity, it's easy to check that $|f(z)|^2y^k=|f(z)|^2(\operatorname{Im}(z))^k$ is $\Gamma$-invariant, i.e. $|f(z)|^2(\operatorname{Im}(z))^k=|f(\gamma z)|^2(\operatorname{Im}(\gamma z))^k$. So, we can unfold the integral obtaining $$F(s)=\int_{\Gamma_\infty\setminus\mathbb{H}}|f(z)|^2y^ky^sd\mu(z)=\int_{0}^\infty\int_{0}^1|f(z)|^2y^ky^sy^{-2}dxdy$$ Using the Fourier expansion of $f(z)$ and integrating with respect to the $x$ variable, we have $$\int_{0}^1|f(z)|^2 dx=\sum_{r}|a_{f}(r)|^2e^{-4\pi ry}$$ So, $$F(s)=\sum_{r}|a_f(r)|^2\int_{0}^\infty y^{k+s-2}e^{-4\pi ry}dy$$ After a change of variable $t=4\pi ry$, we can see that $\int_{0}^\infty y^{k+s-2}e^{-4\pi ry}dy=\frac{\Gamma(k+s-1)}{(4\pi r)^{k+s-1}}$. Then, \begin{align*}F(s)&=\sum_{r}|a_f(r)|^2\frac{\Gamma(k+s-1)}{(4\pi r)^{k+s-1}}\\ &=|a_f(1)|^2\frac{\Gamma(k+s-1)}{(4\pi )^{k+s-1}}\sum_{r}\frac{\lambda_f(r)^2}{r^s}\\ &=|a_f(1)|^2\frac{\Gamma(k+s-1)}{(4\pi )^{k+s-1}}\frac{\zeta(s)}{\zeta(2s)}L(s,\operatorname{sym}^2f)\end{align*} Then $$\operatorname{Res}_{s=1}F(s)=|a_f(1)|^2\frac{\Gamma(k)}{(4\pi )^{k}}\frac{1}{\zeta(2)}L(1,\operatorname{sym}^2f)$$ Equating $$\frac{3}{\pi}=|a_f(1)|^2\frac{\Gamma(k)}{(4\pi )^{k}}\frac{1}{\zeta(2)}L(1,\operatorname{sym}^2f),$$ we get our result.