I'm writing a paper with Eran Assaf, Chan Ieong Kuan, and Alex Walker on the Fibonacci zeta function \begin{equation*} \sum_{n \geq 1} \frac{1}{F(n)^s}, \end{equation*} where $F(n) \DeclareMathOperator{\SL}{SL} \DeclareMathOperator{\GL}{GL} \renewcommand{\Im}{\operatorname{Im}} \renewcommand{\Re}{\operatorname{Re}} \newcommand{\ZZ}{\mathbb{Z}} \DeclareMathOperator*{\Res}{Res}$ is the $n$th Fibonacci number. We show that it has a meromorphic continuation through modular forms. (It's not surprising that it has a meromorphic continuation, but it's not at all obvious why modular forms have anything to do with it).

While generalizing our results, we struggled to find a computational error. We could check this if we could estimate certain Petersson inner products sufficiently closely. Magma can compute some Petersson inner products, but not between the objects we were studying.

In this note,^{1 1}This is also available as pdf note. I show that if we can compute enough Fourier coefficients, we can easily approximate the inner product. I record some observations and simplified experiments.

And I'll note that this was enough to help us resolve our computational error by showing that our theoretical inner product computations agreed with the experimental. The error was elsewhere. I hope that we'll have a preprint out shortly.

Computing Petersson Inner Products via Hidden Unfolding

Suppose we are given two weight $k$ holomorphic cusp forms $f$ and $g$ on $\Gamma_0(N) \subset \SL(2, \ZZ)$, and we want to compute the Petersson inner product \begin{equation}\label{eq:desired} \langle fy^{k/2}, gy^{k/2} \rangle = \iint_{\Gamma_0(N) \backslash \mathcal{H}} f(z) \overline{g(z)} y^k d\mu(z). \end{equation} Numerically computing the Petersson inner product in general is hard.

One approach is to insert the weight $0$ Eisenstein series \begin{equation} E(z, s) := \sum_{\gamma \in \Gamma_\infty \backslash \Gamma_0(N)} \Im(\gamma z)^s. \end{equation} Consider the inner product \begin{equation}\label{eq:instead} \langle f \overline{g} y^k, E(z, \overline{s}) \rangle. \end{equation} Writing $f = \sum a(n) q^n$ and $g = \sum b(n) q^n$, the standard unfolding argument with Eisenstein series shows that \begin{equation} \langle f \overline{g} y^k, E(z, \overline{s}) \rangle = \sum_{n \geq 1} a(n) \overline{b(n)} \int_0^\infty e^{-4\pi n y} y^{s + k - 1} \frac{dy}{y} = \frac{\Gamma(s + k - 1)}{(4\pi)^{s + k - 1}} \sum_{n \geq 1} \frac{a(n) \overline{b(n)}}{n^{s + k - 1}}. \end{equation} The point is that the Eisenstein series has a pole at $s = 1$ with constant residue $c$. Taking the residue shows that \begin{equation} c\langle f y^{k/2}, \overline{g} y^{k/2} \rangle = \Res_{s = 1} \langle f \overline{g} y^k, E(z, \overline{s}) \rangle = \frac{\Gamma(k)}{(4\pi)^{k}} \Res_{s = 1} \sum_{n \geq 1} \frac{a(n) \overline{b(n)}}{n^{s + k - 1}}. \end{equation}

Stated differently, computing \eqref{eq:desired} can be transformed into a problem about computing a residue of, essentially, a Rankin–Selberg $L$-function.

I also remark that if $f$ and $g$ are merely cusp forms (and not holomorphic cusp forms), analogous work applies and the result is the same as long as the resulting integral over $y$ is understandable. For example, if $f$ and $g$ are Maass forms, similar results hold.

If $f$ and $g$ have different eigenvalues, then $\langle f, g \rangle = 0$. If they have the same eigenvalue, then the resulting $y$ integral is understandable (in the sense that one can look it up in a table of integrals or ask Mathematica). (Our motivating problem occurred when $f$ and $g$ are different Maass forms with the same eigenvalue).

Numerically estimating residues of $L$-functions

We now look at the reduced problem. Given a Dirichlet series \begin{equation} D(s) = \sum_{n \geq 1} \frac{a(n)}{n^s} \end{equation} that is presumed to have meromorphic continuation to $\Re s > 1 - \delta$ for some $\delta > 0$ and a simple pole at $s = 1$, how can we numerically estimate the residue?

Direct estimation

We begin very naively.

In a neighborhood of $1$, we have \begin{equation}\label{eq:laurent} D(s) = \frac{r}{s - 1} + f(s) \end{equation} for a holomorphic function $f$. Thus \begin{equation} \lim_{s \to 1} (s - 1) D(s) = r, \end{equation} the residue. If we can compute $D(s)$ along a sequence of points tending to $1$, then we might just try to estimate this limit.

In practice, $D(s)$ often converges absolutely for $\Re s > 1$, so we choose some small $\epsilon > 0$ and try looking at \begin{equation} (1 + \epsilon - 1) D(1 + \epsilon) = \epsilon \cdot D(1 + \epsilon) \approx r. \end{equation}

In practice, this behaves very badly. Let's look at $\zeta(s)$ and use the obvious approximation \begin{equation} \zeta(1 + \epsilon) = \sum_{n \geq 1} \frac{1}{n^{1 + \epsilon}} \approx \sum_{n \leq N} \frac{1}{n^{1 + \epsilon}}. \end{equation}

def zeta_approx(N, s):
    return sum(1/(m**n(s)) for m in range(1, N + 1))

for epspower in (-2, -3, -4, -5, -6, -7):
    eps = 10**epspower
    s = 1 + eps
    for npower in (3, 4, 5, 6):
        N = 10**npower
        print(f"N {N}, eps {eps}: {n(eps * zeta_approx(N, s))}")

This produces the following output:

N 1000, eps 1/100: 0.0725297980739927
N 10000, eps 1/100: 0.0937690500268027
N 100000, eps 1/100: 0.114528539813696
N 1000000, eps 1/100: 0.134815847783707

N 1000, eps 1/1000: 0.00746173653014950
N 10000, eps 1/1000: 0.00974539312390463
N 100000, eps 1/1000: 0.0120241987617199
N 1000000, eps 1/1000: 0.0142978033378315

N 1000, eps 1/10000: 0.000748309249198581
N 10000, eps 1/10000: 0.000978337305495686
N 100000, eps 1/10000: 0.00120835285756249
N 1000000, eps 1/10000: 0.00143831949805017

N 1000, eps 1/100000: 0.0000748523297425421
N 10000, eps 1/100000: 0.0000978718262077023
N 100000, eps 1/100000: 0.000120894841456149
N 1000000, eps 1/100000: 0.000143917731532429

N 1000, eps 1/1000000: 7.48544707142625e-6
N 10000, eps 1/1000000: 9.78756369334501e-6
N 100000, eps 1/1000000: 0.0000120900799291495
N 1000000, eps 1/1000000: 0.0000143926313619477

N 1000, eps 1/10000000: 7.48546848163299e-7
N 10000, eps 1/10000000: 9.78760180176266e-7
N 100000, eps 1/10000000: 1.20901395097692e-6
N 1000000, eps 1/10000000: 1.43927171867352e-6

We cannot usually expect better behavior than from $\zeta(s)$, and this looks abysmal! Thinking more about it, we can understand some of this behavior. Clearly the coefficients are all $1$, and adding up $N$ terms $1/n^{1 + \epsilon}$ will be at most $N$. It doesn't make sense to choose $N$ and $\epsilon$ such that $N \epsilon \ll 1$: the resulting residue approximation will trivially tend to $0$. This explains why the estimated residues are getting smaller as $1 + \epsilon \to 1$.

Of course, we know how to handle $\zeta(s)$ very well. But this illustrates the problem: if understanding the behavior of the pole is hard, then understanding behavior \underline{near} the pole might also be hard. Similarly, if we don't know how to understand the pole purely abstractly, we presumably don't have particularly good ways to understand the meromorphic continuation. We're then trying to study the continuation in its worst region of convergence.

This isn't the way.

Partial sums

Let's use a different strategy that avoids estimating $D(s)$ directly. Perron's formula says that \begin{equation} \sum_{n \leq X} a(n) = \frac{1}{2 \pi i} \int_{(\sigma)} D(s) \frac{X^s}{s} ds, \end{equation} where $\sigma$ is large enough to be in the region of absolute convergence. We can be more precise and give a truncated Perron formula (see for example section 5.3 of Montgomery-Vaughan) of the form \begin{equation} \sum_{n \leq X} a(n) = \frac{1}{2 \pi i} \int_{\sigma - iT}^{\sigma + iT} D(s) \frac{X^s}{s} ds + R(X, t) \end{equation} for a boundable remainder $R(X, t)$. If $\lvert a(n) \rvert \ll 1$ (or indeed, even if this holds \emph{on average}), the remainder is bounded by $X \log X / T$.

Using handwavy analysis, we should expect that we can shift $\sigma$ slightly left of $1$, picking up a pole at $s = 1$ from $D(s)$. The top and bottom contours are approximately bounded by $X \log X / T$. The left contour is annoying and really depends on our understanding of the meromorphic properties of $D(s)$. \emph{Usually}, it has some reasonable growth property that means for some optimal choice of $T$, the residue dominates the integral and remainder.

Then we (usually) have a relationship of the form \begin{equation} \sum_{n \leq X} a(n) = \Res_{s = 1} D(s) \frac{X^s}{s} + o(X) = X r, \end{equation} using $r$ from the Laurent expansion \eqref{eq:laurent}. Thus \begin{equation} r \approx \frac{1}{X} \sum_{n \leq X} a(n), \end{equation} with convergence depending on the actual bounds for the shifted integrals in Perron's analysis.

For $\zeta(s)$, this is of course perfect. The sum of the integers up to $X$ is exactly $X$, so this recovers the residue exactly. Nonetheless, sometimes convergence is slow — and proving convergence requires some finicky analysis (very reminiscent of some proofs of the prime number theorem).

Smoothed partial sums

Alternately, we can use more forgiving approximations than Perron. For example, we can use the Riesz mean (sometimes called the Ces`{a}ro integral transform, and I don't know the actual history), \begin{equation} \sum_{n \leq X} a(n) (1 - n/X)^A = \frac{1}{2 \pi i} \int_{\sigma} D(s) X^s \frac{\Gamma(1 + A) \Gamma(s)}{\Gamma(1 + A + s)} ds. \end{equation} The left hand side remains easy to compute: it's a finite sum. The point is that if $D(s)$ is polynomially bounded in vertical strips, then one can choose $A \gg 1$ large enough that the integral converges beautifully even after shifting to the left.

The residue theorem then gives \begin{align} \sum_{n \leq X} a(n) (1 - n/X)^A &= \Res_{s = 1} D(s) X^s \frac{\Gamma(1 + A) \Gamma(s)}{\Gamma(1 + A + s)} + \int_{(\sigma')} D(s) X^s \frac{\Gamma(1 + A) \Gamma(s)}{\Gamma(1 + A + s)} ds \\ &= \frac{rX}{1 + A} + O_A(X^{\sigma'}), \end{align} where we assume $A$ is large enough to guarantee absolute convergence of the shifted integral. (In practice, we take $A > B + 1$, where $\lvert D(\sigma' + it) \rvert \ll (1 + \lvert t \rvert)^B$ describes the polynomial bound for $D$ and $0 < \sigma' < 1$).

For example, consider $\zeta(s)$ for convenience, we have \begin{equation} \sum_{n \leq X} (1 - n/X)^2 = \frac{X}{3} + O(X^{\epsilon}), \end{equation} as we know the trivial convexity estimate \begin{equation} \lvert \zeta(\epsilon + it) \rvert \ll (1 + \lvert t \rvert)^{\frac{1}{2}}. \end{equation}

Explicitly computing, we find

tot = 0
X = 100000
for m in range(1, X + 1):
    tot += n(1 - m/X)**2
residue = 3 / X * tot
print(residue)
# 0.999985000050005

That is, \begin{equation} \frac{3}{X} \sum_{n \leq X} (1 - n/X)^2 = 0.999985\ldots, \end{equation} and so the residue of $\zeta(s)$ at $1$ is approximately $1$.

After choosing $A$ large enough for initial convergence, choosing larger smoothing powers doesn't necessarily correspond to more rapid convergence. Tail behavior is determined by the growth of the shifted integral. While there is a difference between "quickly growing" and "decaying", there is much less difference between "decaying" and "decaying more".

More technically, the function \begin{equation*} D(s) X^s \frac{\Gamma(1 + A)\Gamma(s)}{\Gamma(1 + A + s)} \end{equation*} has a pole at $s = 0$ from the gamma functions (and depending on the Dirichlet series $D(s)$, possibly also a pole from $D(s)$ at $s = 0$). If $A$ is sufficiently large and $D(s)$ sufficiently well-behaved that we can shift the line of integration to the left of $0$, then the secondary contribution coming from this second pole will control convergence — and this behavior is essentially independent of $A$.

Rigorous estimates

To make these estimately rigorous, we would need an absolute bound or understanding of the constants from the $O(\cdot)$ bounds above.

In the case of automorphic $L$-functions, it would be possible to convert the Riesz mean approach described above into rigorous estimates for the residues (and thus rigorous estimates for the Petersson inner products). To make this fully rigorous, we would either need to have a proved bound for the coefficients $\lvert a(n) \rvert$ or accept a Ramanujan-Petersson type conjecture.

We can do any $\GL(2)$ type Hecke modular form, for example: using either Deligne's result on $\lvert a(p) \rvert$ or the Kim–Sarnak bounds.

Coefficient bounds are also usually available for $L$-functions coming from varieties; the Hecke bound from elliptic curves, for example. (Of course we expect these to all be automorphic in some sense, but what we expect and what we can prove differ).

The structure of making this rigorous would look like

1. Give explicit bounds for $\lvert a(n) \rvert$.
2. Use this to bound $\lvert L(2 + it)\rvert$.
3. Produce the corresponding (weak exponent) convexity estimate on the interval $-1 \leq \Re s \leq 2$. Make the constant dependence explicit!
4. Choose $A \gg 1$ large enough to compensate for the weak convexity estimate growth.
5. Bound the resulting shifted integral along the line $\Re s = \epsilon$.

I note that I use a wider interval for the convexity argument than the typical $[0, 1]$ or $[-\epsilon, 1 + \epsilon]$ argument. This is because we are not optimizing for the best-possible polynomial growth in the strip: we are instead optimizing for easy-to-understand bounds for $L(2 + it)$ and simple, explicit coefficient growth.

A different way to produce rigorous estimates is to study instead contributions coming from a particularly nice basis of forms, and then to express everything in terms of that basis. This is approximately the strategy used in PARI/GP's implementation of Petersson inner products.

Leave a comment

Comment via email