Numerical evaluation of the Petersson product of elliptic modular forms

It is known how to compute the Fourier expansion of elliptic modular forms using modular symbols, and it is known how to get numerical evaluations of $L$-functions of various type ; it's possible to get explicit values on those matters with sage already.

It is also known that when an Eisenstein series is involved, it's possible to relate the Petersson scalar product to $L$-functions, and hence to evaluate them.

I have seen that sage bug about various pairings for modular forms, but it looks more like it's about the pairing between modular forms and modular symbols than the Petersson scalar product.

My question is: does there exist general formulas to compute the Petersson scalar product of two elliptic modular forms numerically?

EDIT(2012-12-23): I insist on the numerically: having an expansion with estimates on the order of the error with constants which depends on this or that (I'm thinking about those which can be found in chapter 5 of Iwaniek's "Topics in classical automorphic forms" for example) is very nice from a theoretical point of view, but doesn't help when one wants to actually compute with specific forms and to a given precision. In fact, I want to compute various things with the Petersson scalar product, so this question is to check whether I can directly work on them or if I should write something about the matter before.

• for Jacobi forms, see arxiv.org/abs/1009.3198: "a numerical method to compute the Petersson scalar products of Jacobi Forms is developed and discussed in detail." – Carlo Beenakker Dec 1 '12 at 20:24

There is a "quick and dirty" way to find the inner product of two cusp forms that are not necessarily Hecke eigenforms. I learned this from Akshay Venkatesh.

The formula is that \begin{equation*}\langle f, g \rangle = \lim_{y \rightarrow 0^+} y^k \int_0^{1} f(x+iy) \overline{g(x+iy)} dx, \end{equation*} where $f$ and $g$ are weight $k$ and the inner product is normalized via \begin{equation*} \langle f, g \rangle = \int_{\Gamma \backslash \mathbb{H}} y^k f(z) \overline{g(z) }\frac{1}{V} \frac{dx dy}{ y^2}, \end{equation*} where $V$ is the volume of $\Gamma \backslash \mathbb{H}$. The philosophy behind the proof is that the horocycle $x+iy: 0 \leq x \leq 1$ equidistributes in the fundamental domain as $y \rightarrow 0$. You can prove the formula by spectrally decomposition $f \overline{g}$. The projection onto the constant eigenfunction gives $\langle f, g\rangle$. The projections onto the cusp forms integrate out to zero. The projection onto the Eisenstein series leaves the constant terms which are bounded by $\sqrt{y}$, and hence have limit zero as $y$ tends to $0$.

If $f(z) =\sum_n a(n) e(nz)$ and $g(z) = \sum_n b(n) e(nz)$, then of course \begin{equation*} \int_0^{1} f(x+iy) \overline{g(x+iy)} dx = \sum_{n \geq 1} a(n) \overline{b(n)} \exp(-4 \pi n y). \end{equation*}

• This formula looks nice from a theoretical point of view, but having a limit on an integral isn't a good start for numerical evaluation... I'm looking for something like a series expansion with an explicit bound on the remainder when cut. My question is really about practical numerical computation. – Julien Puydt Dec 2 '12 at 17:35
• @Julien : Actually, my guess is that this kind of limit formula will give a good approximation in practice (because of the sum converging exponentially), but that it will not obvious to bound the error term rigorously. – François Brunault Dec 2 '12 at 21:15
• The last expression is the one you want to use, which does not involve an integral. Instead you suppose you have computed Fourier coefficients up to some bound $n \leq N$, and then choose $y$ so that the sum of the tail is quite small, using whatever available bound you have for the Fourier coefficients. It's not the best algorithm, but it's easy to implement. – Matt Young Dec 3 '12 at 1:50
• @MattYoung I just reread your comment... In the first formula I'm supposed to take the limit as y goes to zero, so I can hardly fix it high enough to get a good series convergence in the second expression ;-) – Julien Puydt Nov 27 '13 at 17:02
• Of course, but for small values of $y$ the formula gives a good approximation to the Petersson inner product. The number of Fourier coefficients you need to compute and sum is going to be proportional to $1/y$. It takes some additional work to explicitly bound the difference between $\langle f, g \rangle$ and $\sum_{n \leq N} a(n) \overline{b(n)} \exp(-4 \pi n y)$. – Matt Young Nov 27 '13 at 17:25

Let $f(z) = \sum a(n) e(n z)$ and $g(z) = \sum b(n) e(n z)$ be holomorphic modular forms of weight $k \in 2 \mathbb{N}$ on $\Gamma = \operatorname{SL}_2(\mathbb{Z})$ whose product decays rapidly. Then

$$\int_{\Gamma \backslash \mathbb{H}} y^k \overline{f(z)}g(z) ~ \frac{dx ~ d y}{y^2} = 2 \sum_{n \in \mathbb{N} } \frac{ \overline{a(n)} b(n) }{ n^{k-1} } \sum_{d \in \mathbb{N}} \Phi(4 \pi d \sqrt{n}),$$

where $$\Phi(y) = 2( \frac{y}{8 \pi})^{k-1} (y K_{k-2}(y) - K_{k-1}(y)).$$ Note that $\Phi(y) \asymp_k y^{k-1/2} e^{-y}$ for $y \gg 1$. For a general finite index subgroup $\Gamma$ of $\operatorname{SL}_2(\mathbb{Z})$, a correct formula may be obtained by summing the RHS over the cusps $\mathfrak{a}$ of $\Gamma$, weighted by the width $w$ of $\mathfrak{a}$, and taking for $a(n), b(n)$ ($n \in w^{-1} \mathbb{N}$) the Fourier coefficients at $\mathfrak{a}$.

A reference for a general form of such identities is Theorem 5.6 (p.24) in my paper Evaluating modular forms on Shimura curves; see also Example 5.7, Remark 3.5, and the discussion of Sections 5.3--5.6, which includes a detailed comparison with the other approaches mentioned in this thread that I will summarize briefly here. Pros: no need to Hecke-decompose, unlike the "symmetric square" approach; converges rapidly to the correct value, unlike the vanilla "equidistribution of horocycles" approach; generalizes to non-holomorphic forms lacking a straightforward Hecke decomposition (e.g., certain theta series), although perhaps this feature is not important for your purposes. Cons: requires the Fourier expansion at every cusp, unlike either approach just mentioned (although the "symmetric square" approach is not devoid of such subtlety, since it requires one to compute the conductor and bad Euler factors of the symmetric square of a newform).

Another approach (specific to the holomorphic case) would be to exploit the connection with period polynomials, for which a search just now turned up this article. One variant of that method also requires knowing Fourier expansions at every cusp, and reduces the problem to evaluating a class of incomplete gamma functions some of which reduce to K-Bessel functions as above; another requires only that one be able to compute periods of a cusp form $f$ over split geodesics in $\Gamma \backslash \mathbb{H}$, which can apparently be done using the Fourier expansion at only one cusp. Moreover, one can speed up the computation when $f = g$ is an eigenform by exploiting certain rationality results.

It's easy to reduce to the case of computing the Petersson product of a normalised new eigenform with itself. Here you can use the fact that the product is equal to the value at s=k of the symmetric square L-function of f, and this you can compute using e.g. Tim Dokchitser's algorithms. Here is a thread from the Sage developers mailing list with example code by Martin Raum: https://groups.google.com/forum/m/#!topic/sage-nt/EkBWOogY8yw

For elliptic curves there is also Mark Watkins' Sympow program, which will compute all the symmetric power L-functions.

• Hmmmm... basically decomposing the two forms on an orthogonal basis of normalized eigenforms, thus reducing to an explicit linear combination of squares, which indeed are easy to tackle using existing code? That looks quite promising. – Julien Puydt Dec 2 '12 at 18:40
• Well, you can interpret this in two ways. Firstly, it's saying that the Petersson product of $f$ and $g$ is basically $\lim_{s \to k^+} (s - k)^{-1} \sum_{n \ge 1} a_n \overline{b_n} n^{-s}$. This holds for any $f$, $g$ and is somehow in the same ball-park as Matt Young's answer (but with a less rapidly converging series). But secondly, if you restrict to f=g a new eigenform, then you're computing a value of an $L$-function, and there is a well-developed theory for rigorous numerical computation of L-values things. – David Loeffler Dec 3 '12 at 9:05