Computing the Petersson norm of newforms of weight 2 from the symmetric square $L$-function Let $f \in S_2(\Gamma_0(N))$ be a newform with trivial character. I want to compute the Petersson norm $\lVert f\rVert^2$ of $f$, not normalized by $1/[\operatorname{SL}_2(\mathbf{Z}):\Gamma_0(N)]$, as in Gross–Zagier.
From Numerical evaluation of the Petersson product of elliptic modular forms, I came across the formula $$\lVert f\rVert^2 = \frac{(k-1)!}{2^{2k - 1}\pi^{k + 1}}L(\operatorname{Sym}^2(f),2)$$ with $k = 2$. I implemented this in Magma. If $N$ is not square-free, I guess the correct Euler factor for the symmetric square at $p^2 \mid N$ by testing if the functional equation for the symmetric square is satisfied with $1 \pm x$ or $1 \pm px$.
However, comparing with the result of PARI/gp (which is normalized by $1/[\operatorname{SL}_2(\mathbf{Z}):\Gamma_0(N)]$, so removing that normalization), it seems that I have to multiply my result by $N$ if $N$ is square-free, the reason for which I don't understand (maybe it's a convention of the implementation of the symmetric square $L$-function in Magma?). It is even worse for $N$ not square-free, e.g. $f \in S_2(\Gamma_0(125))^+$ or $S_2(\Gamma_0(147))^{w_3,w_{49}}$, where the normalization factors seem to be $125$ and $147 \cdot 7/8$, respectively.
PARI's code is hard to read, and Petersson scalar products are not implemented for $N \neq 1$ in Sage.
Can someone please shed light on this?
 A: Let $f$ be a newform of weight $k$, level $q$, and nebentypus $\chi$, where $\chi$ is a primitive Dirichlet character modulo $q_1 \mid q$, and let $f(z) = \sum_{n = 1}^{\infty} \lambda_f(n) n^{\frac{k - 1}{2}} e(nz)$ be the Fourier expansion of $f$, where the Hecke eigenvalues $\lambda_f(n)$ are normalised such that $|\lambda_f(n)| \leq d(n)$. Let $d\mu(z) = \frac{dx \, dy}{y^2}$ and let $E(z,s) = \sum_{\gamma \in \Gamma_{\infty} \backslash \Gamma_0(q)} \Im(\gamma z)^s$ denote the real analytic Eisenstein series for $\Gamma_0(q)$. Then by unfolding,
$$\int_{\Gamma_0(q) \backslash \mathbb{H}} |y^{k/2} f(z)|^2 E(z,s) \, d\mu(z) = \sum_{n = 1}^{\infty} |\lambda_f(n)|^2 n^{k - 1} \int_{0}^{\infty} y^{s + k - 1} e^{-4\pi ny} \, \frac{dy}{y}.$$
We make the change of variables $y \mapsto y/(4\pi n)$. Since $\Gamma(s) = \int_{0}^{\infty} y^s e^{-s} \, \frac{dy}{y}$, we arrive at the identity
$$\int_{\Gamma_0(q) \backslash \mathbb{H}} |y^{k/2} f(z)|^2 E(z,s) \, d\mu(z) = (4\pi)^{1 - s - k} \Gamma(s + k - 1) \sum_{n = 1}^{\infty} \frac{|\lambda_f(n)|^2}{n^s}.$$
We take the residue of both sides at $s = 1$. Since the residue of $E(z,s)$ at $s = 1$ is
$$\frac{1}{\operatorname{vol}(\Gamma_0(q) \backslash \mathbb{H})} = \frac{1}{[\mathrm{SL}_2(\mathbb{Z}) : \Gamma_0(q)] \operatorname{vol}(\mathrm{SL}_2(\mathbb{Z}) \backslash \mathbb{H})} = \frac{3}{\pi q \prod_{p \mid q} (1 + p^{-1})},$$
we deduce that
$$\int_{\Gamma_0(q) \backslash \mathbb{H}} |y^{k/2} f(z)|^2 \, d\mu(z) = \frac{\pi q \prod_{p \mid q} (1 + p^{-1}) \Gamma(k)}{3 (4\pi)^k} \operatorname*{Res}_{s = 1} \sum_{n = 1}^{\infty} \frac{|\lambda_f(n)|^2}{n^s}.$$
Next, we use the fact that
$$\sum_{n = 1}^{\infty} \frac{|\lambda_f(n)|^2}{n^s} = \frac{\zeta^q(s) L^q(s,\operatorname{ad} f)}{\zeta^q(2s)} \prod_{p \mid q} \sum_{r = 0}^{\infty} \frac{|\lambda_f(p^r)|^2}{p^{rs}},$$
where I write $L^q(s,\pi)$ to denote the $L$-function with the Euler factors dividing $q$ omitted, and $L(s,\operatorname{ad} f)$ denotes the adjoint $L$-function; if the nebentypus of $f$ is trivial, this is the same as the symmetric square $L$-function $L(s,\operatorname{sym}^2 f)$.
Since $\operatorname*{Res}_{s = 1} \zeta^q(s) = \prod_{p \mid q} (1 - p^{-1})$, whereas $\zeta^q(2) = \frac{\pi^2}{6} \prod_{p \mid q} (1 - p^{-2})$, we arrive at the identity
$$\int_{\Gamma_0(q) \backslash \mathbb{H}} |y^{k/2} f(z)|^2 \, d\mu(z) = \frac{2 q\Gamma(k)}{\pi (4\pi)^k} L^q(1,\operatorname{ad} f) \prod_{p \mid q} \sum_{r = 0}^{\infty} \frac{|\lambda_f(p^r)|^2}{p^r}.$$
To simplify this any further, we need to use the local theory of automorphic representations of $\mathrm{GL}_2(\mathbb{A}_{\mathbb{Q}})$ at each prime $p \mid q$.

*

*If $p \parallel q$ but $p \nmid q_1$, then the local component $\pi_p$ at $p$ of the automorphic representation $\pi_f$ associated to $f$ is a twist-minimal special representation (i.e. a special representation associated to an unramified character of $\mathbb{Q}_p^{\times}$), in which case $|\lambda_f(p^r)| = p^{-r/2}$, so that $\sum_{r = 0}^{\infty} |\lambda_f(p^r)|^2 p^{-r} = (1 - p^{-2})^{-1}$. In this case, we have that $L_p(1,\operatorname{ad} f) = (1 - p^{-2})^{-1}$ as well.

*If $p \mid q$ but $p \nmid \frac{q}{q_1}$, so that $p^k \parallel q$ and $p^k \parallel q_1$, then $\pi_p$ is a twist-minimal principal series representation (i.e. $\pi_p = \omega_1 \boxplus \omega_2$ with one of $\omega_1,\omega_2$ unramified), in which case $|\lambda_f(p^r)| = 1$, so that $\sum_{r = 0}^{\infty} |\lambda_f(p^r)|^2 p^{-r} = (1 - p^{-1})^{-1}$. In this case, we have that $L_p(1,\operatorname{ad} f) = (1 - p^{-2})^{-1}$ as well.

*In all remaining cases (so that $\pi_p$ is supercuspidal or a non-twist-minimal principal series representation or a non-twist-minimal special representation), $\lambda_f(p^r) = 0$ for all $r \geq 1$, so that $\sum_{r = 0}^{\infty} |\lambda_f(p^r)|^2 p^{-r} = 1$. However, we do not necessarily have that $L_p(1,\operatorname{ad} f) = 1$ as well, and in fact this only occurs when $\pi_p$ is supercuspidal and not twist-invariant by the unramified quadratic character of $\mathbb{Q}_p^{\times}$.

Here we can determine the local factors at $s = 1$ of $L(s,\operatorname{ad} f) = L(s,f \otimes \widetilde{f})/\zeta(s)$ via work of Gelbart and Jacquet (namely Corollary (1.3) and Proposition (1.4)).
A: My guess is that the formula you are trying to use is only valid for $N=1$, and thus needs correction in general.
Maybe Shimura's paper can help sort this out. https://doi.org/10.1002/cpa.3160290618
In (2.1), which Shimura writes for $\Gamma_1(N)$ but that doesn't matter when the character is trivial, his definition is
$$\langle f,f\rangle={3/\pi\over [SL_2(Z):\Gamma_0(N)]}\int_\Phi |f|^2 dx dy.$$
Then in (2.5) you have
$$\langle f,f\rangle={\Gamma(2)\over (4\pi)^2}\cdot\mathop{\rm res}\limits_{s=2} D(s,f,f)$$
where by the last display of Section 1 he defines
$$D(s,f,f)=\sum_{n=1}^\infty {a_n^2\over n^s}.$$
Now a comparison of Euler products gives that the local factors
of $D(s,f,f)\zeta(2s-2)$ and $L(s,Sym^2 f)\zeta(s-1)$ match, at least away from $p$ that divide $N$ (this discrepancy is the issue that David Loeffler raises). This Euler product comparison is mentioned in another paper of Shimura, see (0.4) of https://doi.org/10.1112/plms/s3-31.1.79
Anyway, this gives the answer up to the bad factors,
namely
$$\int_\Phi |f|^2 dx dy={[SL_2(Z):\Gamma_0(N)]\over 3/\pi}\langle f,f\rangle$$
$$=[SL_2(Z):\Gamma_0(N)]{\pi\over 3}{1\over (4\pi)^2}\mathop{\rm res}\limits_{s=2} D(s,f,f)$$
$$=[SL_2(Z):\Gamma_0(N)]{\pi\over 48\pi^2}{1\over\zeta(2)}L(2,Sym^2 f)\prod_{p|N} C_p$$
$$=N\prod_{p|N}(1+1/p)\cdot{1\over 8\pi^3}L(2,Sym^2 f)\prod_{p|N}C_p$$
Note that this matches your asserted formula when $N=1$ and $k=2$.
In the more general case, considering a bad prime $p|N$, the Euler factor from $\zeta(s-1)/\zeta(2s-2)$ evaluated at $s=2$ exactly cancels out factor of $(1+1/p)$ in the index formula. Meanwhile, the Euler factor of $L(s,Sym^2f)$ when $p$ exactly divides $N$ is $(1-1/p^s)^{-1}$, as is the Euler factor of $D(s,f,f)$ in this case (since $a_p^2=1$).
Finally, when $p^2|N$, the Euler factor of $D(s,f,f)$ is trivial since $a_p^2=0$, while that of $L(s,Sym^2f)$ can be known either by theory or trial-and-error computation. For the theoretical side, one can presuably work with the $p$-minimal twist of $f$ where this minimality allows twists with nontrivial Nebentypus - see 2.1 of Coates and Schmidt, particularly (2.12). https://doi.org/10.1515/crll.1987.375-376.104
I think one aspect is that if $v_p(N)$ is odd then the Euler factor of $L(s,Sym^2f)$ is trivial; while if $v_p(N)$ is even and $f$ is itself $p$-minimal then the factor is $(1+p/p^s)^{-1}$; and otherwise the Euler factor comes from that of the $p$-minimal twist (though perhaps not completely transparently, again with this $(1+p/p^s)^{-1}$ possibly appearing).
