Consider the Voronoi summation formula for the sum of squares function $r_2$, in terms of Bessel function $J_0$: $$\sum_{n=0}^\infty r_2(n) \int_0^\infty \pi J_0(2\pi\sqrt{nx}) f(x) \, dx = \sum_{n=0}^\infty r_2(n) f(n) \tag{1}$$

To me, it seems it might be helpful to look at it from the point of view of distributions. Let $K(n,x) = \pi J_0(2\pi\sqrt{nx})$, and $a(n) = r_2(n)$. The Voronoi summation formula becomes: $$\sum_{n=0}^\infty a(n) K(n, x) = \sum_{n=0}^\infty a(n) \delta(x-n) \tag{2}$$ note this is valid for $x \geq 0$ (the summation diverges when $x$ is negative). And there is some Gibbs phenomenon going on. Let's take a picture of $\sum_{n=0}^{500} a(n) K(n, x)$:

And compare it with $a(n)$:

Formula $(2)$ is similar to the case of Poisson summation, where $K(n,x) = e^{-2\pi i n x}$，$a(n) = 1$, and the summation is from $-\infty$ to $\infty$: $$\sum_{n \in \mathbb{Z}} e^{-2\pi i n x} = \sum_{n \in \mathbb{Z}} \delta(x-n) \tag{3}$$

The Poisson summation formula can be made exactly like $(2)$, if we let $a(0) = \frac{1}{2}$, and $a(n) = 1$ otherwise: $$\sum_{n=0}^\infty a(n) e^{-2\pi i n x} = \sum_{n=0}^\infty a(n) \delta(x-n) \tag{4}$$ which is also only valid for $x \geq 0$.

Finally, if we define an integral transform: $$\tilde{f}(n) = \int_0^\infty K(n,x) f(x) \, dx \tag{5}$$ then $(2)$ reduces to an elegant statement, that is, the weighted Dirac comb $\delta_{a} = \sum_{n=0}^\infty a(n) \delta(x-n)$ satisfies (for $x \geq 0$): $$\widetilde{\delta_{a}} = \delta_{a} \tag{6}$$

Here comes the interesting part. The Voronoi summation formula is usually proven using the functional equation of the Dirichlet series $\zeta_{r_2}(s) = 4 \, \zeta(s) \, L(\chi_4, s)$ associated to $r_2$ (by some applications of Mellin transform). And the reverse is true too, i.e., we can prove the functional equation using the summation fomula.

So my idea is, maybe we can prove the functional equation (or, modularity) of the Dirichlet series (or, modular form) associated to $a(n)$ by proving $(2)$ for some integral kernel $K(n,x)$. For example, for the weight $k$ and level $0$ case, we can normalize $a(n)$ to $a(n) n^{-\frac{k-1}{2}}$ (in general, it seems we shall normalize $a(n)$ to be of $O(1)$ growth), and then: $$K(n,x) = 2 \pi (-1)^{\frac{k}{2}} J_{k-1}(4\pi\sqrt{nx})$$ For the level $N$ case, we need some more modifications.

And this might be helpful, because the Poisson summation formula can be proven directly using other methods (besides using the functional equation of $\zeta$ function). So one can hope we can use similar methods to prove $(2)$ for other $a(n)$, and hence the functional equation (or, modularity) of $a(n)$ (the only problem is, whether it will be simpler, or it will ends up in some kind of circular reasoning).

**So my questions is, is the above formulation well-known and well-researched?** If not so, I definitely would like to look into it. Unfortunately I am no longer in the mathematics circle (working in industry now). So I think I shall share what I found at the moment with everyone here. To me it seems there are three methods to attack $(2)$.

Firstly, the Voronoi summation formula $(1)$ can be proved directly using $2$-d Poisson summation: $$\sum_{(x,y) \in \mathbb{Z}^2} \hat{F}(x,y) = \sum_{(x,y) \in \mathbb{Z}^2} F(x,y) $$ and if $F(x,y) = f(\sqrt{x^2+y^2})$, i.e., has radical symmetry, then we have: $$\sum_{n=0}^\infty r_2(n) \breve{f}(\sqrt{n}) = \sum_{n=0}^\infty r_2(n) f(\sqrt{n})$$ where $\breve{f}$ is the Hankel transform: $$\breve{f}(n) = 2\pi \int_0^\infty J_0(2\pi n x) f(x) x \, dx$$ which means: $$\breve{f}(\sqrt{n}) = \pi \int_0^\infty J_0(2\pi\sqrt{nx}) f(\sqrt{x}) \, dx$$ and hence $(1)$ is proved. The problem is, whether this can be done for other $a(n)$.

Indeed this can be done. As the most simple example, let: $$a(n) = \#\{ x^2 + k y^2 = n \,|\, x, y \in \mathbb{Z}\}$$ for some $k>0$, $k \in \mathbb{Z}$. Then let our $F(x,y) = f(\sqrt{x^2+k y^2})$, and after working out the new "Hankel transform" (a fun exercise of polar coordinates), we'll see our $K(n,x)$ is $\pi J_0(2\pi\sqrt{\frac{nx}{k}})$. Of course, this is usually proven using theta function, but the above method is still cool (personally, I think the nicest part is, we can simply plot $(2)$ in Mathematica [try it if you haven't] to verify whether our computation is correct).

In fact, it seems that, for any $\rho: \mathbb{R}^n \to \mathbb{R}_{\geq 0}$ satisfying:

(1) $\rho$ is a Morse function with one single critical point at $0$. Although it seems something can still be done when there are multiple critical points, let's keep it simple.

(2) The image of $\mathbb{Z}^n \subset \mathbb{R}^n$ is still in $\mathbb{Z}_{\geq 0}$.

(3) Notice $\rho$ is invariant under the Galois action of the permutation of roots of $\rho(v) = n$.

Then we can have a summation formula, where $a(n) = \#\{ \rho(v) = n \,|\, v \in \mathbb{Z}^n\}$. Note this may or may out descend to a functional equation of the L-function in the usual sense, probably because the summation formula is in some sense more general than the functional equation.

At this point, the natural thing to do, is to see whether this can work for $a(n)$ coming from point counting of varieties in $\mathbb{F}_p$, for instance, point counting of a polynomial $f(x)$ with a $S_3$ Galois group. I played around with it a bit, and it seems it might work if we consider point counting of $f(x) - n p = 0$ in $\mathbb{Z}^3$ (where $n$ and $p$ are new dummy variables, in some sense a bit like blowing up $f(x)$. reminds us of the folklore that $\mathbb{Z}$ is $3$-d...) instead of $f(x) = 0$ (p.s. I am still working on it, so is it a good idea to talk about it? As I am writing this, I realize there might indeed be some chance that this gives us the twists required for the level N summation formula when we extract the $\mathbb{F}_p$ information from the solutions of $f(x) - n p = 0$, because the information will be periodic, and hence we shall hit it with characters). The nice thing is the Galois group of $f$ put constraints on $\rho$, so we are using the Galois group.

[UPDATE: It seems we shall first see whether this method can prove Kronecker-Weber, as a sanity check. I am working on it now...]

[**UPDATE**: Probably the idea shall be put in this way: The modularity of theta functions, or, a(n) coming from solution counting of quadratic forms, are easy to prove, because it's related to point counting in lattice. But the modularity of other arithmetic a(n) coming from Galois representation, and in the simplest case, from $a(p) = \#\{f(x) = 0 $ in $ F_p\} - 1$ (and other a(n) calculated from a(p) by considering Hecke operator in the usual way), is difficult in general. So it will be lucky if we can still transfer the problem to some kind of point counting in lattice.

It seems that it might help if we consider $f(x) - n p = 0$ in $\mathbb{Z}^3$, for example. If f(x) has a solution in $F_p$, i.e. there is a grid point (x, n, p) for some n, then (x+kp, n+stuff, p) is also a grid point for any k. So the grid point density in $x$-direction is related to $a(p)$, and it's also periodic. Moreover, if the Galois group of $f(x)$ is Abelian, then due to Kronecker-Weber we will see periodicity in $p$-direction as well.

Because we hope to duplicate the proof of the $a(n) = r_2(n)$ case, if we write down all the details in the proof, we will notice a few things. Firstly, $F(x,y)$ shall be invariant in each level set of $\rho$, such that it descends to a function $f(n)$ on level sets. This is similar to the fact that $\sqrt{x^2 + y^2}$ is invariant for all solutions of $x^2 + y^2 = n$, such that we can collect the terms and write $r_2(n) \cdot f(\sqrt{n})$. So it's more like saying the level sets contain some kind of Galois orbits, although the Galois action usually cannot descend to $\mathbb{R}^n$.

Secondly, the proof of the $a(n) = r_2(n)$ case, relies on the $(r, \theta)$ polar coordinate, where $\theta$ is actually the coordinate in each level set, so it will be easier if our level sets are connected, hence the requirement on critical points.]

Secondly, there are two methods to prove the Poisson summation formula directly (which may leads to two other methods to prove $(2)$). Firstly, by harmonic analysis on $\mathbb{R}$ and $\mathbb{Z}$ (Fourier series). So this leads to representation theory / non-abelian harmonic analysis / trace formula methods. For example, will our $K(n,x)$ always be the character / matrix coefficient of some representation? This seems likely, because firstly Bessel functions are indeed from representation theory. (I remember seeing somewhere that trace formulas can also be seen from analyze of Dirac distributions, but I forget the reference.)

And the second method is by brute-force. To me, it seems the Poisson summation formula actually comes from this "Euler-Maclaurin" statement: $$\sum_{n = -N}^{N} e^{-2\pi i n x} = \frac{\sin((2N+1)\pi x)}{\sin(\pi x)} = \frac{\pi x}{\sin(\pi x)} \int_{-N-\frac{1}{2}}^{N+\frac{1}{2}} e^{-2\pi i nx} \, dn \tag{7}$$ which is valid for $|x| < 1$, such that $\frac{\pi x}{\sin(\pi x)}$ has no poles. (side remark: $\sum_{n = -N}^{N} e^{-2\pi i n x}$ is the character of a $SU(2)$ representation. I wonder whether this is a coincidence.)

Let $N \to \infty$ and we have (also valid for $|x| < 1$) the "orthogonality" statement: $$\sum_{n \in \mathbb{Z}} e^{-2\pi i n x} = \frac{\pi x}{\sin(\pi x)} \int_{-\infty}^{\infty} e^{-2\pi i nx} \, dn = \frac{\pi x}{\sin(\pi x)} \delta(x) = \delta(x) \tag{8}$$

Finally note $\sum_{n \in \mathbb{Z}} e^{-2\pi i n x}$ is invariant under $x \to x+1$ ("modular symmetry"), so we have $\sum_{n \in \mathbb{Z}} e^{-2\pi i n x} = \sum_{n \in \mathbb{Z}} \delta(x-n)$.

Unfortunately, moving this to $K(n,x)$, such as $\pi J_0(2\pi\sqrt{nx})$, will run into problems (in some sense this is good because this means something interesting is going on). Let's see what we get. Let's try to prove (when $|x| < 1$): $$\sum_{n=0}^\infty a(n) K(n, x) = a(0) \delta(x) \tag{9}$$

"Euler-Maclaurin": $$ \sum_{n = 0}^{N} a(n) K(n,x) = R(N, x) \int_{0}^{N} a(n) K(n,x) \, dn \tag{10}$$ This seems fine, because we just need $Q(x) = \lim_{N \to \infty} R(N, x)$ to be pole-free when $|x| < 1$, and $\lim_{x \to 0} Q(x) = 1$. Here we can take $a(n) = a([n])$ where $[n]$ is round to the nearest integer.

"Orthogonality": $$\int_{0}^{\infty} a(n) K(n,x) \, dn = a(0) \delta(x) \tag{11}$$ This seems tough. Let: $$\Delta(N, x) = \int_{0}^{N} a(n) K(n,x) \, dn \tag{12}$$ Now for some $x \neq 0$, $|x| < 1$, we'll find $\Delta(N, x)$ oscillating around $0$ as $N \to \infty$. For example, take $K(n,x) = \pi J_0(2\pi\sqrt{nx})$, $x=\frac{1}{2}$, then here is $\Delta(N, x)$ for $N$ from $0$ to $500$:

A wild guess is it is bounded by $O(N^{\frac{1}{2} + \epsilon})$, although this might ends up as difficult as RH, or, just the same as proving the modularity of $a(n)$. Will some stationary phase etc. method work?

For the Poisson summation case, the $\Delta(N, x)$ is trivially bounded, which is why it is so simple.

- Finally, we no longer have the simple $x \to x+1$ symmetry for $\delta_{a}$, so we can't directly go from $|x| < 1$ to other $x$.

On the other hand, $r_2$ and $J_0$ have lots of symmetries, such as $r_2$'s connection with theta function (now I am afraid we'll run into circular reasoning again) and the PDEs satisfied by $J_k$ and its various series and integral representations. It remains to see how this can help us to solve 2. and 3.

Well, this is a very long post. And I am sure there'll be mistakes and loopholes. Thank you for reading and let me know what you think.