# An asymptotic expansion for a Lambert series associated to Siegel cusp forms

Babita  Babita
Department of Mathematical Sciences
Indian Institute of Technology(Banaras Hindu University), Varanasi
221005, Uttar Pradesh, India<EMAIL_ADDRESS>,  Abhash Kumar Jha
Abhash Kumar Jha
Department of Mathematical Sciences
Indian Institute of Technology(Banaras Hindu University), Varanasi
221005, Uttar Pradesh, India<EMAIL_ADDRESS>,  Abhishek Juyal
Abhishek Juyal
Department of Mathematics,
HNB Garhwal University, BCC Campus,
Srinagar, Uttarakhand, 246174<EMAIL_ADDRESS>and  Bibekananda Maji
Bibekananda Maji
Department of Mathematics
Indian Institute of Technology Indore
Indore, Simrol, Madhya Pradesh 453552, India<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

In 2000, Hafner and Stopple proved a conjecture of Zagier which states that
the constant term of the automorphic function $|\Delta(x+iy)|^{2}$ i.e., the
Lambert series $\sum_{n=1}^{\infty}\tau(n)^{2}e^{-4\pi ny}$ can be expressed
in terms of the non-trivial zeros of the Riemann zeta function. In this
article, we study a certain Lambert series associated to Siegel cusp forms and
observe a similar phenomenon.

2010 Mathematics Subject Classification. Primary 11M06, 11M26; Secondary
11N37.
Keywords and phrases. Riemann zeta function, non-trivial zeros, Lambert
series, Rankin-Selberg $L$-function, Siegel cusp forms

##  1\. Introduction

Let $\zeta(s)$ and $\mu(n)$ denote the Riemann zeta function and the Möbius
function, respectively. It is well-known that the reciprocal of the Riemann
zeta function is the Dirichlet series associated to $\mu(n).$ Ramanujan, in
his second notebook [34, p. 312], [5, eq. (37.3), p. 470] mentioned the
following identity associated to $\mu(n)$:

For a positive real number $x$,

$\displaystyle\sum_{n=1}^{\infty}\frac{\mu(n)}{n}\exp\left({-\frac{x}{n^{2}}}\right)=\sqrt{\frac{\pi}{x}}\sum_{n=1}^{\infty}\frac{\mu(n)}{n}\exp\left(-\frac{\pi^{2}}{n^{2}x}\right).$

However, the above identity is not correct. The corrected version of the above
identity was obtained by Hardy and Littlewood [18, p. 156]. Assuming all the
non-trivial zeros of $\zeta(s)$ are simple, they proved that the following
identity holds:

$\displaystyle\sqrt{\alpha}\sum_{n=1}^{\infty}\frac{\mu(n)}{n}\exp\left({-\frac{1}{\pi}\left(\frac{\alpha}{n}\right)^{2}}\right)$
$\displaystyle-\sqrt{\beta}\sum_{n=1}^{\infty}\frac{\mu(n)}{n}\exp\left({-\frac{1}{\pi}\left(\frac{\beta}{n}\right)^{2}}\right)$
$\displaystyle=-\frac{\sqrt{\pi}}{2\sqrt{\beta}}\sum_{\rho}\frac{\Gamma\left(\frac{1-\rho}{2}\right)}{\zeta^{\prime}(\rho)}\left(\frac{\beta}{\sqrt{\pi}}\right)^{\rho},$
(1.1)

where $\alpha$ and $\beta$ are two positive real numbers with $\alpha\beta=1$
and the sum over $\rho$ runs through all the non-trivial zeros of $\zeta(s)$.
The convergence of the series (1) is quite intricate. In the same paper, Hardy
and Littlewood mentioned that they were unable to prove the convergence of
this series even after assuming the Riemann hypothesis. However, in an
indirect way, they showed the convergence of this series under the assumption
of the following bracketing condition on the non-trivial zeros of $\zeta(s)$,
namely, the non-trivial zeros $\rho_{1}$ and $\rho_{2}$ of $\zeta(s)$ are
included in the same bracket if they satisfy

$\displaystyle|\operatorname{Im}(\rho_{1})-\operatorname{Im}(\rho_{2})|<\exp\left(-\frac{A_{0}\operatorname{Im}(\rho_{1})}{\log(\operatorname{Im}(\rho_{1}))}\right)+\exp\left(-\frac{A_{0}\operatorname{Im}(\rho_{2})}{\log(\operatorname{Im}(\rho_{2}))}\right),$
(1.2)

where $A_{0}$ is some positive constant. The identity (1) also inspired Hardy
and Littlewood to obtain the following equivalent criterion for the Riemann
hypothesis (RH):

$\displaystyle\sum_{n=1}^{\infty}\frac{\mu(n)}{n}\exp\left({-\frac{x}{n^{2}}}\right)=O\left(x^{-\frac{1}{4}+\epsilon}\right),\quad\mathrm{as}\,\,x\rightarrow\infty,$
(1.3)

Over the time, the above equivalent criterion for RH has inspired
mathematicians to find different generalization of (1) and (1.3) for
$L$-functions associated to several objects. For more details on this topic we
refer to [1, 4, 11, 12, 13, 15, 16]. We would like to mention here that, we
also encounter an infinite series similar to the series (1) in the main result
of this paper.

Let $\Delta(z)$ be the Ramanujan delta cusp form. The Fourier series expansion
of $\Delta(z)$ is given by

$\Delta(z)=\sum_{n=1}^{\infty}\tau(n)\exp(2n\pi iz),~{}~{}z\in\mathcal{H}.$

The Fourier coefficients $\tau(n)$ is known as Ramanujan tau function. In
1981, Zagier [36, p. 417], [37, p. 271] speculated that the constant term of
the automorphic function $\mathcal{F}(z):=y^{12}|\Delta(x+iy)|^{2}$, that is,
the following Lambert series

$c_{0}(y):=y^{12}\sum_{n=1}^{\infty}\tau^{2}(n)\exp({-4\pi ny}),$ (1.4)

has an oscillatory behaviour when $y\rightarrow 0^{+}$. Moreover, Zagier
conjectured that its asymptotic expansion can be written in terms of the non-
trivial zeros of $\zeta(s)$. Another interesting point he observed that one
can numerically obtain non-trivial zeros of $\zeta(s)$ using the asymptotic
expansion of (1.4) and the Ramanujan tau function. Further, he claimed that
the constant term $c_{0}(y)$ satisfies the following asymptotic expansion:

$\displaystyle
c_{0}(y)\sim\mathcal{A}+\sum_{\rho}y^{1-\frac{\rho}{2}}A_{\rho},$

where the sum over $\rho$ runs through the non-trivial zeros of $\zeta(s)$,
$\mathcal{A}=\frac{3}{\pi}\langle\Delta,\Delta\rangle$, and the constant
$A_{\rho}$ depends on $\rho$. Assuming the Riemann hypothesis, the above
asymptotic expansion and the oscillatory behaviour of $c_{0}(y)$ have been
proved by Hafner and Stopple [17]. A similar phenomenon for Hecke eigenforms
for the full modular group as well as for the congruence subgroups have been
observed in [9] and [10]. Subsequently, an asymptotic expansion of an infinite
series associated to a Hecke-Maass eigenforms has been studied by Banerjee and
Chakraborty [6]. Recently, this problem in the case of Hilbert modular forms
has been examined by Agnihotri [2]. To know more about similar problems
inspired from the above conjecture of Zagier, we refer to [21, 22, 29, 30].

Siegel modular forms are generalization of elliptic modular forms to higher
dimension. They were first introduced by Siegel to study quadratic forms. A
Siegel modular forms admits a Fourier-Jacobi expansion and the Fourier-Jacobi
coefficients are Jacobi forms. Motivated from the above mentioned conjecture
of Zagier, in this paper, we study certain Lambert series associated with the
Fourier-Jacobi coefficients of Siegel modular forms. Interestingly, a similar
phenomenon has been observed as in the case of elliptic modular forms. We use
the analytic properties of certain Dirichlet series associated with Siegel
modular forms (studied by Kohnen and Skoruppa [23]), bounds of the Whittaker
function, and a special evaluation of the Meijer $G$-function to prove our
result.

In the next section, we provide necessary backgrounds which are essential to
define the main problem of our article.

##  2\. Preliminaries

Let $\mathbb{C}$ and $\mathcal{H}$ denote the complex plane and complex upper
half-plane, respectively. For two matrices $A$ and $B$ of appropriate sizes,
we define $A[B]:={B}^{t}AB,$ where $B^{t}$ is the transpose of the matrix $B.$
We define $J:=\left(\begin{array}[]{cc}0&{I}\\\ {-I}&0\end{array}\right),$
where $I$ denotes the $2\times 2$ identity matrix. The Siegel upper half plane
$\mathcal{H}_{2}$ of degree $2$ is defined as

$\mathcal{H}_{2}:=\\{Z\in M_{2\times
2}(\mathbb{C})~{}|~{}Z^{t}=Z,\operatorname{Im}(Z)~{}\text{is positive
definite}\\}.$

The full Siegel modular group $\Gamma_{2}:=Sp_{4}(\mathbb{Z})$ of genus $2$ is
defined as

$\Gamma_{2}:=\left\\{M\in M_{4\times 4}(\mathbb{Z}),J[M]=J\right\\}.$

If we write the matrix $M\in\Gamma_{2}$ as a block decomposition
$M=\begin{pmatrix}A&B\\\ C&D\\\ \end{pmatrix}$, where $A,B,C,D\in M_{2\times
2}(\mathbb{Z}),$ then we have

$\Gamma_{2}:=\left\\{\begin{pmatrix}A&B\\\ C&D\\\
\end{pmatrix}\bigg{|}A,B,C,D\in M_{2\times
2}(\mathbb{Z}),A^{t}C=C^{t}A,B^{t}D=D^{t}B,A^{t}D-B^{t}C=I_{2}\right\\}.$

The group $\Gamma_{2}$ acts on $\mathcal{H}_{2}$ as follows:

$\begin{pmatrix}A&B\\\ C&D\\\ \end{pmatrix}\cdot Z=(AZ+B)(CZ+D)^{-1}.$

######  Definition 2.1.

A complex-valued holomorphic function $F:\mathcal{H}_{2}\rightarrow\mathbb{C}$
is said to be a Siegel modular form of weight $k\in\mathbb{N}$ and degree two
if it satisfies

$F(M\cdot Z)={\det}(CZ+D)^{k}F(Z),$

for all $M=\begin{pmatrix}A&B\\\ C&D\\\ \end{pmatrix}\in\Gamma_{2}$,
$Z\in\mathcal{H}_{2},$ and having a Fourier expansion of the form

$F(Z)=\sum_{\begin{subarray}{c}T\geq 0\end{subarray}}A_{F}(T)e^{2\pi
i(tr(TZ))},$ (2.1)

where the summation runs over positive semidefinite half-integral $2\times 2$
matrices $T=(t_{ij})$ ($i.e.,2t_{ij},\;t_{ii}\in\mathbb{Z}$) and $tr(\cdot)$
denotes the trace function.

Further, we say $F$ is a Siegel cusp form if and only if the summation in
(2.1) runs over positive definite half-integral matrices $T.$

We denote the space of Siegel modular forms and Siegel cusp forms of weight
$k$ and degree $2$ on $\Gamma_{2}$ by $M_{k}(\Gamma_{2})$ and
$S_{k}(\Gamma_{2}),$ respectively. The Petersson scalar product on
$S_{k}(\Gamma_{2})$ is defined as

$\langle
F,G\rangle:=\int_{\Gamma_{2}\setminus\mathcal{H}_{2}}F(Z)\overline{G(Z)}(det~{}Y)^{k}dZ,$

where $Z=X+iY$ and $dZ=(det~{}Y)^{-3}dXdY$ is an invariant measure under the
action of $\Gamma_{2}$ on $\mathcal{H}_{2}.$ The space
$(S_{k}(\Gamma_{2}),\langle,\rangle)$ is a finite dimensional Hilbert space.
For more details on the theory on Siegel modular forms, we refer to [8, 32].

###  2.1. Fourier-Jacobi expansion

Let $F(Z)\in M_{k}(\Gamma_{2})$ be a Siegel modular form with the Fourier
series expansion as in (2.1).

Write $Z=\begin{pmatrix}\tau&z\\\ z&\tau^{\prime}\\\
\end{pmatrix}\in\mathcal{H}_{2}$ where
$\tau,\tau^{\prime}\in\mathcal{H},~{}z\in\mathbb{C}.$ A matrix $T$ in the
summation (2.1) can be written as $T=\begin{pmatrix}m&r/2\\\ r/2&n\\\
\end{pmatrix}$ with $m,n\in\mathbb{N},r\in\mathbb{Z}$ and $4mn-r^{2}\geq 0.$
Then the Fourier series expansion (2.1) becomes

$\displaystyle F(Z)=\sum_{T\geq 0}A_{F}(T)\exp({2\pi i(tr(TZ))})$
$\displaystyle=$ $\displaystyle\sum\limits_{\begin{subarray}{c}m,n,r\\\
4mn-r^{2}\geq 0\end{subarray}}A_{F}\left(\begin{pmatrix}m&r/2\\\ r/2&n\\\
\end{pmatrix}\right)\exp({2\pi i(m\tau+rz+n\tau^{\prime})})$ $\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}\phi_{n}(\tau,z)\exp({2\pi
in\tau^{\prime}}),$

where

$\displaystyle\phi_{n}(\tau,z)=\sum_{\begin{subarray}{c}m,r\\\ 4mn-r^{2}\geq
0\end{subarray}}A_{F}\left(\begin{pmatrix}m&r/2\\\ r/2&n\\\
\end{pmatrix}\right)\exp({2\pi i(m\tau+rz)}).$

The coefficient $\phi_{n}(\tau,z)$ is called the $n$th Fourier-Jacobi
coefficients of $F$ and is a Jacobi form of weight $k$ and index $n$. Jacobi
forms are natural generalization of modular forms. Jacobi forms were first
systematically studied by Eichler and Zagier and they played a key role in the
proof of Saito-Kurokawa conjecture. For more details on the theory of Jacobi
forms, we refer to [14]. It is easy to check that if $F\in S_{k}(\Gamma_{2}),$
then its Fourier-Jacobi expansion is given by
$F(Z)=\sum\limits_{n=1}^{\infty}\phi_{n}(\tau,z)\exp({2\pi in\tau^{\prime}})$
and each of the Fourier-Jacobi coefficients are Jacobi cusp forms.

###  2.2. Dirichlet series associated to Siegel cusp forms

We now define the Rankin-Selberg type Dirichlet series associated to Siegel
cusp forms $F_{1},F_{2}\in S_{k}(\Gamma_{2})$ with the Fourier-Jacobi
expansion as follows:

$F_{1}(Z)=\sum\limits_{n=1}^{\infty}\phi_{n}(\tau,z)\exp({2\pi
in\tau^{\prime}}),~{}~{}~{}~{}~{}F_{2}(Z)=\sum\limits_{n=1}^{\infty}\psi_{n}(\tau,z)\exp({2\pi
in\tau^{\prime}})$ (2.2) $\displaystyle
D_{F_{1},F_{2}}(s):=\zeta(2s-2k+4)\sum_{n=1}^{\infty}\frac{\langle\phi_{n},\psi_{n}\rangle}{n^{s}},$
(2.3)

where $\langle\phi_{n},\psi_{n}\rangle$ denotes the Petersson scalar product
(defined in [14, p. 27, eq. (13)]) of Jacobi cusp forms $\phi_{n}$ and
$\psi_{n}.$ It is well-known [23, p. 544, Lemma 1] that the Petersson scalar
product of the Fourier-Jacobi coefficients $\phi_{n}$ and $\psi_{n}$ satisfies

$\langle\phi_{n},\psi_{n}\rangle=O(n^{k}).$

Further, Kohnen [25, p. 718], [28, p. 134] made the following Ramanujan-
Petersson conjecture: for any $\epsilon>0$,

$\displaystyle\langle\phi_{n},\phi_{n}\rangle=O_{F_{1},\epsilon}\left(n^{k-1+\epsilon}\right).$
(2.4)

Kohnen and Sengupta [27] proved the above conjecture for Hecke eigenform $F\in
S_{k}(\Gamma_{2})$ which is a Saito-Kurokawa lift. Recently, Kumar and Paul
[31] proved that the above conjecture is true on average for arbitrary degree.
Assuming (2.4), one can clearly see that the Dirichlet series
$D_{F_{1},F_{2}}(s)$ is absolutely convergent for $\operatorname{Re}(s)>k$.

In the next section we state the main results of our paper.

##  3\. Main Results

As already mentioned in the introduction, we are interested to study certain
Lambert series associated to Siegel cusp forms. More precisely, for
$\alpha>0$, we study the following Lambert series:

$\displaystyle\sum_{n=1}^{\infty}\langle\phi_{n},\psi_{n}\rangle\exp(-4\pi
n\alpha),$ (3.1)

where $\phi_{n}$ and $\psi_{n}$ are defined as in (2.2). Before stating the
main result of this article, we define an arithmetic function
$a_{F_{1},F_{2}}(n)$ that satisfy the following generating function:

$\displaystyle\sum_{n=1}^{\infty}\frac{a_{F_{1},F_{2}}(n)}{n^{s}}=\frac{D_{F_{1},F_{2}}(s)}{\zeta(2s+1-2k)},$
(3.2)

which is absolutely convergent for $\operatorname{Re}(s)>k$. Let
$W_{\kappa,\mu}(z)$ be the well-known Whittakar function which is one of the
solutions of the following second order differential equation:

$\displaystyle\frac{d^{2}w}{dz^{2}}+\left(\frac{1-4\mu^{2}+4\kappa
z-z^{2}}{4z^{2}}\right)w=0.$

With the above notation and definitions in our hand, we are ready to state the
main result of this article.

######  Theorem 3.1.

Let $F_{1}$ and $F_{2}$ be two Siegel cusp forms of weight $k$ and degree $2$
with the Fourier-Jacobi expansion given as in (2.2) and satisfying (2.4). Let
$\alpha$ and $\beta$ be two positive real numbers such that $\alpha\beta=1$.
Under the assumption of the simplicity hypothesis of the non-trivial zeros of
$\zeta(s)$, we have

$\displaystyle\sum_{n=1}^{\infty}\langle\phi_{n},\psi_{n}\rangle\exp(-4\pi
n\alpha)$
$\displaystyle=\frac{\beta^{2k-2}}{\pi^{3/2}}\sum_{n=1}^{\infty}\frac{a_{F_{1},F_{2}}(n)}{(4\pi
n\beta)^{\frac{k-1}{2}}}\exp(-2\pi n\beta)W_{\frac{k}{2},\frac{k}{2}-1}(4\pi
n\beta)$
$\displaystyle+R_{k}+\sum_{\rho}\frac{\Gamma\left(\frac{\rho}{2}+k-2\right)D_{F_{1},F_{2}}\left(\frac{\rho}{2}+k-2\right)(4\pi\alpha)^{2-k-\frac{\rho}{2}}}{\zeta^{\prime}(\rho)},$

where the sum over $\rho$ runs through non-trivial zeros of $\zeta(s)$ and
satisfy the bracketing condition (1.2) and the term $R_{k}$ is given by

$\displaystyle R_{k}=\begin{cases}0,&\text{ if }\langle
F_{1},F_{2}\rangle=0,\\\ \frac{90\langle
F_{1},F_{2}\rangle}{\pi^{2}\alpha^{k}},&\text{ if }\langle
F_{1},F_{2}\rangle\neq 0.\end{cases}$ (3.3)

As an application of the above theorem, we have the following asymptotic
expansion of the Lambert series defined in (3.1).

######  Corollary 3.2.

Let $F_{1}$ and $F_{2}$ be as in Theorem 3.1. If $\langle
F_{1},F_{2}\rangle\neq 0$, then for $\alpha\rightarrow 0^{+}$, we have

$\displaystyle\alpha^{k}\sum_{n=1}^{\infty}\langle\phi_{n},\psi_{n}\rangle\exp(-4\pi
n\alpha)\sim\frac{90\langle F_{1},F_{2}\rangle}{\pi^{2}}.$ (3.4)

In the next section, we collect a few results which will be essential for the
proof of our main results.

##  4\. Some Well-known Results

We complete the Dirichlet series $D_{F_{1},F_{2}}(s)$ associated to the Siegel
cusp forms $F_{1}$ and $F_{2}$ of weight $k$ and degree $2$ as follows:

$D^{*}_{F_{1},F_{2}}(s):=(2\pi)^{-2s}\Gamma(s)\Gamma(s-k+2)D_{F_{1},F_{2}}(s).$

The analytic properties of the above completed Dirichlet series is given
below.

######  Theorem 4.1.

[23] The Dirichlet series $D_{F_{1},F_{2}}(s)$ can be analytically continued
to an entire function if $\langle F_{1},F_{2}\rangle=0$. However, if $\langle
F_{1},F_{2}\rangle\neq 0$, then $D_{F_{1},F_{2}}(s)$ has a simple pole at
$s=k$ with residue $\frac{4^{k}\pi^{k+2}}{(k-1)!}\langle F_{1},F_{2}\rangle$.
The functional equation of $D_{F_{1},F_{2}}(s)$ is given by

$D^{*}_{F_{1},F_{2}}(s)=D_{F_{1},F_{2}}^{*}(2k-2-s).$ (4.1)

The above result will play a crucial role in the proof of the Theorem 3.1.
Next, we define a well-known special function known as Meijer $G$-function.

Let $m,n,p,q$ be integers with $0\leq m\leq q$, $0\leq n\leq p$. Let
$a_{1},\cdots,a_{p}$ and $b_{1},\cdots,b_{q}$ be $p+q$ complex numbers such
that $a_{i}-b_{j}\not\in\mathbb{N}$, for $i\in[1,n]$, $j\in[1,m]$. The Meijer
$G$-function [33, p. 415] is defined by the following line integral:

$\displaystyle G_{p,q}^{\,m,n}\\!\left(\,\begin{matrix}a_{1},\cdots,a_{p}\\\
b_{1},\cdots,b_{q}\end{matrix}\;\Big{|}z\right)=\frac{1}{2\pi
i}\int_{L}\frac{\prod_{j=1}^{m}\Gamma(b_{j}+w)\prod_{j=1}^{n}\Gamma(1-a_{j}-w)z^{-w}}{\prod_{j=m+1}^{q}\Gamma(1-b_{j}-w)\prod_{j=n+1}^{p}\Gamma(a_{j}+w)}\mathrm{d}w.$
(4.2)

Here, we assume that the line of integration $L$ separates the poles of the
factors $\Gamma(b_{j}+w)$ from the poles of the factors $\Gamma(1-a_{j}-w)$.
The above integral converges absolutely if $2(m+n)>p+q$ and
$|\arg(z)|<(2m+2n-p-q)\frac{\pi}{2}$. We now state the following special
evaluation of the Meijer $G$-function.

######  Lemma 4.2.

[20, p. 58] For $|\arg(z)|<\frac{\pi}{2}$, we have

$\displaystyle G_{1,2}^{\,2,0}\\!\left(\,\begin{matrix}a_{1}\\\
b_{1},b_{2}\end{matrix}\;\Big{|}z\right)=z^{\frac{b_{1}+b_{2}-1}{2}}e^{-z/2}W_{\kappa,\mu}(z),$
(4.3)

where $W_{\kappa,\mu}(z)$ is the Whittaker function with the variables
$\kappa=\frac{1}{2}(b_{1}+b_{2}-1)-a_{1}$ and $\mu=\frac{1}{2}(b_{1}-b_{2})$.

The next lemma tells us about the behaviour of the Whittaker function for
large argument.

######  Lemma 4.3.

[33, p. 69, eq. (2.11.29)] As $z\rightarrow\infty$ with
$|\arg(z)|<3\frac{\pi}{2}$, one has

$\displaystyle W_{\kappa,\mu}(z)\sim\exp(-z/2)z^{\kappa}.$

######  Lemma 4.4 (Stirling’s bound for $\Gamma(z)$).

[19, p. 151, A.4.] Let $z=\sigma+i{T}$ with $p\leq\sigma\leq q.$ Then we have

$|\Gamma(\sigma+i{T})|\sim\sqrt{2\pi}|{T}|^{\sigma-1/2}e^{-\frac{1}{2}\pi|{T}|},\quad{\rm
as}\quad|{T}|\rightarrow\infty.$ (4.4)

######  Lemma 4.5.

[35, p. 219] Let $T_{n}$ be a sequence of positive real numbers such that
$|T_{n}-\operatorname{Im}(\rho)|>\exp\left(-C\frac{\operatorname{Im}(\rho)}{\log(\operatorname{Im}(\rho))}\right)$
for any non-trivial zeros $\rho$ of $\zeta(s)$, where $C$ is some positive
constant. Then

$\displaystyle\frac{1}{|\zeta(x+iT_{n})|}<\exp(C_{1}T_{n}),$

where $0<C_{1}<\pi/4$.

In the next section, we present a proof of the Theorem 3.1.

##  5\. Proof of main results

  * Proof of Theorem 3.1.

Using the inverse Mellin transform of the gamma function $\Gamma(s)$, we see
that the inverse Mellin transform of the function
$\frac{\Gamma(s)D_{F_{1},F_{2}}(s)(4\pi)^{-s}}{\zeta(2s-2k+4)}$ is the
infinite series
$\sum\limits_{n=1}^{\infty}\langle\phi_{n},\psi_{n}\rangle\exp(-4\pi
n\alpha).$ Therefore, for some large positive constant $c>k$, we have

$\displaystyle\sum_{n=1}^{\infty}\langle\phi_{n},\psi_{n}\rangle\exp(-4\pi
n\alpha)=\frac{1}{2\pi
i}\int_{(c)}\frac{\Gamma(s)D_{F_{1},F_{2}}(s)(4\pi\alpha)^{-s}}{\zeta(2s-2k+4)}\mathrm{d}s.$
(5.1)

Here and throughout the article the symbol $(c)$ denotes the line integral
$c-i\infty$ to $c+i\infty$. Now we recall from the Theorem 4.1 that the
Rankin-Selberg Dirchlet series $D_{F_{1},F_{2}}(s)$ can be extended to an
entire function if $\langle F_{1},F_{2}\rangle=0$. However, if $\langle
F_{1},F_{2}\rangle\neq 0$, then $D_{F_{1},F_{2}}(s)$ has a simple pole at
$s=k$ with residue $\frac{4^{k}\pi^{k+2}}{(k-1)!}\langle F_{1},F_{2}\rangle$.
Moreover, from the functional equation of $D_{F_{1},F_{2}}(s)$ it is clear
that the poles of $\Gamma(s)$ are neutralized by $D_{F_{1},F_{2}}(s)$.
Nevertheless, the non-trivial of zeros of $\zeta(2s-2k+4)$ will give us
infinitely many poles of the integrand function in the strip
$k-2<\operatorname{Re}(s)<k-\frac{3}{2}$, and the trivial zeros of
$\zeta(2s-2k+4)$ are at $k-n$ for $n\geq 3$. Thus, we construct a closed
rectangular contour $\mathcal{C}$ consisting of the end points
$c-iT,c+iT,c_{1}+iT,c_{1}-iT$, where $T$ is some large positive constant.
Here, we consider $c_{1}\in(k-3,k-2)$ so that the non-trivial zeros of
$\zeta(2s-2k+4)$ lie inside the contour $\mathcal{C}$ and the trivial zeros of
$\zeta(2s-2k+4)$ lie outside the contour $\mathcal{C}$. Now, applying the
Cauchy’s residue theorem, we obtain

$\frac{1}{2\pi
i}\int_{\mathcal{C}}\frac{\Gamma{(s)}D_{F_{1},F_{2}}(s)(4\pi\alpha)^{-s}}{\zeta(2s-2k+4)}\mathrm{d}s=R_{k}+\sum_{|\operatorname{Im}(\rho)|<T}R_{\rho},$
(5.2)

where the term $R_{k}$ denotes the residual term at $s=k$ and $R_{\rho}$
denotes the residual term at the non-trivial zero $\rho$ of $\zeta(s)$. The
residual term $R_{k}$ is given by

$\displaystyle R_{k}=\begin{cases}0,&\text{ if }\langle
F_{1},F_{2}\rangle=0,\\\ \frac{90\langle
F_{1},F_{2}\rangle}{\pi^{2}\alpha^{k}},&\text{ if }\langle
F_{1},F_{2}\rangle\neq 0.\end{cases}$ (5.3)

Here, we have used the fact that $\zeta(4)=\pi^{4}/90$. Now we assume the
grand simplicity hypothesis, which states that the non-trivial zeros of the
Riemann zeta function $\zeta(s)$ are simple. This gives us the following
evaluation of $R_{\rho}$,

$\displaystyle
R_{\rho}=\frac{\Gamma\left(\frac{\rho}{2}+k-2\right)D_{F_{1},F_{2}}\left(\frac{\rho}{2}+k-2\right)(4\pi\alpha)^{2-k-\frac{\rho}{2}}}{\zeta^{\prime}(\rho)}.$
(5.4)

Letting $T\rightarrow\infty$ in (5.2) and utilizing Stirling’s bound (4.4) for
$\Gamma(s)$ and Lemma 4.5, one can show that both the horizontal integrals
will vanish. Finally, using (5.1), we obtain

$\displaystyle\sum_{n=1}^{\infty}\langle\phi_{n},\psi_{n}\rangle\exp(-4\pi
n\alpha)=\frac{1}{2\pi
i}\int_{(c_{1})}\frac{\Gamma(s)D_{F_{1},F_{2}}(s)(4\pi\alpha)^{-s}}{\zeta(2s-2k+4)}\mathrm{d}s+R_{k}+\sum_{\rho}R_{\rho},$
(5.5)

where the sum over $\rho$ runs through the non-trivial zeros of $\zeta(s)$.
Here, we emphasize that the sum over $\rho$ is an infinite sum as we know the
existence of infinitely many zeros of $\zeta(s)$ from Hardy’s theorem [35, p.
257]. However, the convergence of this series is quite delicate as we do not
know much about the lower bound for $\zeta^{\prime}(s)$. Thus, under the
assumption of convergence of the sum in (5.5) one can proceed further. Now our
main aim is to simplify the following left vertical integral:

$\displaystyle V_{k}(\alpha):=\frac{1}{2\pi
i}\int_{(c_{1})}\frac{\Gamma(s)D_{F_{1},F_{2}}(s)(4\pi\alpha)^{-s}}{\zeta(2s-2k+4)}\mathrm{d}s.$
(5.6)

From the functional equation of $D_{F_{1},F_{2}}(s)$, i.e., Theorem 4.1, one
obtains

$\displaystyle\Gamma(s)D_{F_{1},F_{2}}(s)=\frac{(2\pi)^{4s-4k+4}\Gamma(2k-s-2)\Gamma(k-s)D_{F_{1},F_{2}}(2k-2-s)}{\Gamma(s-k+2)}.$
(5.7)

Also, the symmetric form of the functional equation of the Riemann zeta
function implies that

$\displaystyle\zeta(2s-2k+4)=\frac{\pi^{2s-2k+4}\Gamma\left(k-s-\frac{3}{2}\right)\zeta(2k-2s-3)}{\sqrt{\pi}\Gamma(s-k+2)}.$
(5.8)

Now substituting (5.7) and (5.8) in (5.6), we obtain

$\displaystyle V_{k}(\alpha)=\frac{\sqrt{\pi}}{16^{k-1}\pi^{2k}}\frac{1}{2\pi
i}\int_{(c_{1})}\frac{\Gamma(2k-2-s)\Gamma(k-s)D_{F_{1},F_{2}}(2k-2-s)}{\Gamma\left(k-s-\frac{3}{2}\right)\zeta(2k-2s-3)}(4\pi\beta)^{s}\mathrm{d}s,$
(5.9)

where $\alpha\beta=1$. To simplify further, we make a change of variable
$w=2k-2-s$ and upon simplification, we obtain

$\displaystyle V_{k}(\alpha)=\frac{\beta^{2k-2}}{\pi^{3/2}}\frac{1}{2\pi
i}\int_{(d_{1})}\frac{\Gamma(w)\Gamma(w-k+2)D_{F_{1},F_{2}}(w)(4\pi\beta)^{-w}}{\Gamma\left(w-k+\frac{1}{2}\right)\zeta(2w+1-2k)}\mathrm{d}w.$
(5.10)

One can easily verify that $k<\operatorname{Re}(w)=d_{1}<k+1$ as we have
assumed $k-3<\operatorname{Re}(s)=c_{1}<k-2$ . This allows us to write the
infinite series expansions for $D_{F_{1},F_{2}}(w)$ and $1/\zeta(2w+1-2k)$ as
both the series are absolutely and uniformly convergent in the region
$\operatorname{Re}(w)>k$. Let us write

$\displaystyle\frac{D_{F_{1},F_{2}}(w)}{\zeta(2w+1-2k)}=\sum_{n=1}^{\infty}\frac{a_{F_{1},F_{2}}(n)}{n^{w}},~{}~{}\operatorname{Re}(w)>k.$
(5.11)

Now plugging the above series expansion (5.11) in (5.10) and interchanging the
summation and integration, we obtain

$\displaystyle
V_{k}(\alpha)=\frac{\beta^{2k-2}}{\pi^{3/2}}\sum_{n=1}^{\infty}a_{F_{1},F_{2}}(n)I_{k}(n,\beta),$
(5.12)

where

$\displaystyle I_{k}(n,\beta):=\frac{1}{2\pi
i}\int_{(d_{1})}\frac{\Gamma(w)\Gamma(w-k+2)(4\pi
n\beta)^{-w}}{\Gamma\left(w-k+\frac{1}{2}\right)}\mathrm{d}w.$ (5.13)

Now invoking the definition (4.2) of the Meijer $G$-function with
$m=2,n=0,p=1,q=2$ and $a_{1}=\frac{1}{2}-k,b_{1}=0,b_{2}=2-k$, we see that the
above integral is the following Meijer $G$-function:

$\displaystyle
I_{k}(n,\beta)=G_{1,2}^{2,0}\left(\begin{array}[]{c}\frac{1}{2}-k\\\
0,2-k\end{array}\Big{|}4\pi n\beta\right).$ (5.16)

One can easily verify that the above Meijer $G$-function does satisfy all the
convergence criteria. Now employing Lemma 4.2, we obtain

$\displaystyle I_{k}(n,\beta)=(4\pi n\beta)^{\frac{1-k}{2}}W_{k/2,k/2-1}(4\pi
n\beta)\exp(-2\pi n\beta).$ (5.17)

Substituting (5.17) in (5.12), the final expression of the left vertical
integral becomes

$\displaystyle V_{k}(\alpha)$
$\displaystyle=\frac{\beta^{2k-2}}{\pi^{3/2}}\sum_{n=1}^{\infty}a_{F_{1},F_{2}}(n)(4\pi
n\beta)^{\frac{1-k}{2}}W_{\frac{k}{2},\frac{k}{2}-1}(4\pi n\beta)\exp(-2\pi
n\beta).$ (5.18)

Now applying the asymptotic expansion of the Whittaker function (Lemma 4.3),
one can easily verify the convergence of the above infinite series since the
arithmetic function $a_{F_{1},F_{2}}(n)$ has polynomial growth and the
Whittaker function decays exponentially. Finally, considering the expression
(5.18) of the left vertical integral $V_{k}(\alpha)$ in (5.5), together with
residual terms (5.3) and (5.4), we complete the proof of Theorem 3.1. ∎

  * Proof of Corollary 3.2.

In Theorem 3.1, we have started with the assumption on $\alpha$, $\beta$ that
$\alpha\beta=1$. Therefore, $\alpha\rightarrow 0^{+}$ implies that
$\beta\rightarrow\infty$. Making use of Lemma 4.3 and multiplying by
$\alpha^{k}$, it is easy to see that

$\displaystyle\frac{\alpha^{k}\beta^{2k-2}}{\pi^{3/2}}\sum_{n=1}^{\infty}\frac{a_{F_{1},F_{2}}(n)}{(4\pi
n\beta)^{\frac{k-1}{2}}}\exp(-2\pi n\beta)W_{\frac{k}{2},\frac{k}{2}-1}(4\pi
n\beta)$
$\displaystyle\ll\beta^{k-2}\sum_{n=1}^{\infty}\frac{a_{F_{1},F_{2}}(n)}{(4\pi
n\beta)^{\frac{k-1}{2}}}\exp(-4\pi n\beta)(4\pi
n\beta)^{\frac{k}{2}}\ll\frac{1}{\beta^{M}},$ (5.19)

for some large positive number $M$. This shows that the above infinite series
vanishes as $\beta\rightarrow\infty$. Again, multiplying by $\alpha^{k}$ in
the infinite residual term in Theorem 3.1, which involves non-trivial zeros of
$\zeta(s)$, we arrive at

$\displaystyle\sum_{\rho}\frac{\Gamma\left(\frac{\rho}{2}+k-2\right)D_{F_{1},F_{2}}\left(\frac{\rho}{2}+k-2\right)(4\pi)^{2-k-\frac{\rho}{2}}\alpha^{2-\frac{\rho}{2}}}{\zeta^{\prime}(\rho)}.$
(5.20)

It is well-known that the non-trivial zeros of $\zeta(s)$ lie in the critical
strip $0<\operatorname{Re}(s)<1$. Therefore, we have
$3/2<\operatorname{Re}\left(2-\frac{\rho}{2}\right)<2$, and this indicates
that the above infinite sum (5.20) vanishes as $\alpha\rightarrow 0^{+}$.
Thus, applying above two observations in Theorem 3.1 and looking at the
residual term (3.3), we complete the proof of (3.4).

∎

##  6\. Acknowledgement

The second author is partially supported by SERB Start-up Research Grant (File
No. SRG/2022/000487) and MATRICS grant (File No. MTR/2022/000659). The last
author’s research supported by SERB MATRICS grant (File No. MTR/2022/000545).
Both the authors sincerely thank SERB for the support.

## References

  * [1] A. Agarwal, M. Garg, and B. Maji, _Riesz-type criteria for the Riemann hypothesis_ , Proc. Amer. Math. Soc., 150 (2022), 5151–5163.
  * [2] R. Agnihotri, _Lambert series associated to Hilbert modular form_ , Int. J. Number Theory, 18 (2022), 1335–1349.
  * [3] A. N. Andrianov, Quadratic Forms and Hecke Operators. Grundlehren der Math. Wissenschaften 286, Springer-Verlag (1987).
  * [4] S. Banerjee and R. Kumar, _Equivalent criterion for the grand Riemann hypothesis associated to Maass cusp forms_ , submitted for publication, 2021. arXiv:2112.08143
  * [5] B. C. Berndt, Ramanujan’s Notebooks, Part V, Springer-Verlag, New York, 1998.
  * [6] S. Banerjee and K. Chakraborty, _Asymptotic behaviour of a Lambert series á la Zagier: Maass case_ , Ramanujan J., 48 (2019), 567–575.
  * [7] S. Böcherer and W. Kohnen, _Estimates for Fourier coefficients of Siegel cusp forms_ , Math. Ann., 297 (1993), no. 3, 499–517.
  * [8] J. H. Bruinier, G. van der Geer, G. Harder, and D. Zagier, The 1-2-3 of modular forms. Lectures from the Summer School on Modular Forms and their Applications held in Nordfjordeid, June 2004. Edited by Kristian Ranestad. Universitext. Springer-Verlag, Berlin, 2008.
  * [9] K. Chakraborty, S. Kanemitsu, and B. Maji, _Modular-type relations associated to the Rankin-Selberg $L$-function_, Ramanujan J., 42 (2017), 285–299.
  * [10] K. Chakraborty, A. Juyal, S. D. Kumar, and B. Maji, _An asymptotic expansion of a Lambert series associated to cusp forms_ , Int. J. Number Theory, 14 (2018), 289–299.
  * [11] A. Dixit, S. Gupta, and A. Vatwani, _A modular relation involving non-trivial zeros of the Dedekind zeta function, and the generalized Riemann hypothesis_ , J. Math. Anal. Appl., 515 (2022), no. 2, 126435.
  * [12] A. Dixit, A. Roy, and A. Zaharescu, _Ramanujan-Hardy-Littlewood-Riesz phenomena for Hecke forms_ , J. Math. Anal. Appl., 426 (2015), 594–611.
  * [13] A. Dixit, A. Roy, and A. Zaharescu, _Riesz-type criteria and theta transformation analogues_ , J. Number Theory, 160 (2016), 385–408.
  * [14] M. Eichler and D. Zagier, The Theory of Jacobi Forms, Progr. in Math. vol. 55, Birkhauser, Boston(1985).
  * [15] M. Garg and B. Maji, _Hardy-Littlewood-Riesz type equivalent criteria for the generalized Riemann hypothesis_ , Monatsh. Math., (2023). https://doi.org/10.1007/s00605-023-01857-8
  * [16] S. Gupta and A. Vatwani, _Riesz type criteria for L-functions in the Selberg class_ , accepted, Canad. J. Math. arXiv:2211.02954.
  * [17] J. Hafner and J. Stopple, _A heat kernel associated to Ramanujan’s tau function,_ Ramanujan J., 4 (2000), 123–128.
  * [18] G. H. Hardy and J. E. Littlewood, _Contributions to the theory of the Riemann zeta-function and the theory of the distribution of primes_ , Acta Math., 41 (1916), 119–196.
  * [19] H. Iwaniec, E. Kowalski, Analytic number theory, American Mathematical Society Colloquium Publications, vol. 53, 2004.
  * [20] S. Kanemitsu and H. Tsukada, Contributions to the theory of zeta-functions. The modular relation supremacy. Series on Number Theory and its Applications, 10. World Sci., Singapore, 2015.
  * [21] A. Juyal, B. Maji, and S. Satyanarayana, _An exact formula for a Lambert series associated to a cusp form and the Möbious function_ , Ramanujan J., 57 (2022), 769–784.
  * [22] A. Juyal, B. Maji, and S. Satyanarayana, _An asymptotic expansion for a Lambert series associated to the symmetric square $L$-function,_ Int. J. Number Theory, 19 (2022), 553–567.
  * [23] W. Kohnen and N. P. Skoruppa, _A certain Dirichlet series attached to Siegel modular forms of degree two,_ Invent. Math., 95 (1989), no. 3, 541–558.
  * [24] W. Kohnen, _Estimates for Fourier coefficients of Siegel cusp forms of degree two. II._ Nagoya Math. J., 128 (1992), 171–176
  * [25] W. Kohnen, _Jacobi forms and Siegel modular forms: recent results and problems_ , Enseign. Math. 39 (1993), 121–136.
  * [26] W. Kohnen, _Estimates for Fourier coefficients of Siegel cusp forms of degree two_ , Compositio Math., 87 (1993), no. 2, 231–240.
  * [27] W. Kohnen, J. Sengupta, _Bounds for Fourier-Jacobi coefficients of Siegel cusp forms of degree two_ , L-functions and automorphic forms, 159–169, (eds J. Bruinier and W. Kohnen; Springer, Berlin, 2017) 159–169.
  * [28] W. Kohnen, _On the growth of the Petersson norms of Fourier-Jacobi coefficients of Siegel cusp forms_ , Bul. Lond. Math. Soc. 43 (2011), 2025-2052.
  * [29] B. Maji, P. Naskar, and S. Sathyanarayana, _A series associated to Rankin-Selberg L-function and modified K-Bessel function_ , submitted for publication, 2023.
  * [30] B. Maji, S. Sathyanarayana, and B. R. Shankar, _An asymptotic expansion for a twisted Lambert series associated to a cusp form and the Mobius function: level aspect_ , Results Math., 77, 123 (2022).
  * [31] B. Kumar and B. Paul, _Ramanujan-Petersson conjecture Fourier-Jacobi coefficient for Siegel cups forms_ , Bul. Lond. Math. Soc., 53 (2021), 274–284.
  * [32] H. Maass, Siegel’s modular forms and Dirichlet series, Lecture Notes in Mathematics, 216, Springer-Verlag, Berlin-New York, 1971.
  * [33] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, C. W. Clark, eds., NIST Handbook of Mathematical Functions, Cambridge University Press, Cambridge, 2010.
  * [34] S. Ramanujan, Notebooks of Ramanujan, Vol 2, Tata Institute of FundamentaI Research, Bombay, 1957.
  * [35] E. C. Titchmarsh, The Theory of the Riemann Zeta-function, Clarendon Press, Oxford, 1986.
  * [36] D. Zagier, _The Rankin-Selberg method for automorphic functions which are not of rapid decay,_ J. Fac. Sci. Univ. Tokyo IA Math. 28 (1981), 415–437.
  * [37] D. Zagier, Introduction to modular forms, From number theory to physics (Les Houches, 1989), Springer, 1992, 238–291.


# Scalability and Sample Efficiency Analysis of Graph Neural Networks for
Power System State Estimation
††thanks: This paper has received funding from the European Union’s Horizon
2020 research and innovation programme under Grant Agreement number 856967.

Ognjen Kundacina  The Institute for Artificial Intelligence
Research and Development of Serbia
Novi Sad, Serbia
<EMAIL_ADDRESS>Gorana Gojic  The Institute for Artificial
Intelligence
Research and Development of Serbia
Novi Sad, Serbia
<EMAIL_ADDRESS>Mirsad Cosovic  Faculty of Electrical Engineering
University of Sarajevo
Sarajevo, Bosnia and Herzegovina
<EMAIL_ADDRESS>Dragisa Miskovic  The Institute for Artificial
Intelligence
Research and Development of Serbia
Novi Sad, Serbia
<EMAIL_ADDRESS>Dejan Vukobratovic  Faculty of Technical Sciences
University of Novi Sad
Novi Sad, Serbia
<EMAIL_ADDRESS>
###### Abstract

Data-driven state estimation (SE) is becoming increasingly important in modern
power systems, as it allows for more efficient analysis of system behaviour
using real-time measurement data. This paper thoroughly evaluates a phasor
measurement unit-only state estimator based on graph neural networks (GNNs)
applied over factor graphs. To assess the sample efficiency of the GNN model,
we perform multiple training experiments on various training set sizes.
Additionally, to evaluate the scalability of the GNN model, we conduct
experiments on power systems of various sizes. Our results show that the GNN-
based state estimator exhibits high accuracy and efficient use of data.
Additionally, it demonstrated scalability in terms of both memory usage and
inference time, making it a promising solution for data-driven SE in modern
power systems.

###### Index Terms:

State Estimation, Graph Neural Networks, Machine Learning, Power Systems,
Real-Time Systems

##  I Introduction

Motivation and literature review: The state estimation (SE) algorithm is a key
component of the energy management system that provides an accurate and up-to-
date representation of the current state of the power system. Its purpose is
to estimate complex bus voltages using available measurements, power system
parameters, and topology information [1]. In this sense, the SE can be seen as
a problem of solving large, noisy, sparse, and generally nonlinear systems of
equations. The measurement data used by the SE algorithm usually come from two
sources: the supervisory control and data acquisition (SCADA) system and the
wide area monitoring system (WAMS) system. The SCADA system provides low-
resolution measurements that cannot capture system dynamics in real-time,
while the WAMS system provides high-resolution data from phasor measurement
units (PMUs) that enable real-time monitoring of the system. The SE problem
that considers measurement data from both WAMS and SCADA systems is formulated
in a nonlinear way and solved in a centralized manner using the Gauss-Newton
method [1]. On the other hand, the SE problem that considers only PMU data
provided by WAMS has a linear formulation, enabling faster, non-iterative
solutions.

In this work, we will focus on the SE considering only phasor measurements,
described with a system of linear equations [2], which is becoming viable with
the increasing deployment of PMUs. This formulation is usually solved using
linear weighted least-squares (WLS), which involve matrix factorizations and
can be numerically sensitive [3]. To address the numerical instability issues
that often arise when using traditional SE solvers, researchers have turned to
data-driven deep learning approaches [4, 5]. These approaches, when trained on
relevant datasets, are able to provide solutions even when traditional methods
fail. For example, in [4], a combination of feed-forward and recurrent neural
networks was used to predict network voltages using historical measurement
data. In the nonlinear SE formulation, the study [5] demonstrates the use of
deep neural networks as fast and quality initializers of the Gauss-Newton
method.

Both linear WLS and common deep learning SE methods at its best approach
quadratic computational complexity regarding the power system size. To fully
utilize high sampling rates of PMUs, there is a motivation to develop SE
algorithms with a linear computational complexity. One way of achieving this
could be using increasingly popular graph neural networks (GNNs) [6, 7]. GNNs
have several advantages when used in power systems, such as permutation
invariance, the ability to handle varying power system topologies, and
requiring fewer trainable parameters and less storage space compared to
conventional deep learning methods. One of the key benefits of GNNs is the
ability to perform distributed inference using only local neighbourhood
measurements, which can be efficiently implemented using the emerging 5G
network communication infrastructure and edge computing [8]. This allows for
real-time and low-latency decision-making even in large-scale networks, as the
computations are performed at the edge of the network, closer to the data
source, reducing the amount of data that needs to be transmitted over the
network. This feature is particularly useful for utilizing the high sampling
rates of PMUs, as it can reduce communication delays in PMU measurement
delivery that occur in centralized SE implementations.

GNNs are being applied in a variety of prediction tasks in the field of power
systems, including fault location [9], stability assessment [10], and load
forecasting [11]. GNNs have also been used for power flow problems, both in a
supervised [12] and an unsupervised [13] manner. A hybrid nonlinear SE
approach [14] combines a model and data-based approach using a GNN that
outputs voltages which are used a regularization term in the SE loss function.

Contributions: In our previous work [15], we proposed a data-driven linear
PMU-only state estimator based on GNNs applied over factor graphs. The model
demonstrated good approximation capabilities under normal operating conditions
and performed well in unobservable and underdetermined scenarios. This work
significantly extends our previous work in the following ways:

  * •

We conduct an empirical analysis to investigate how the same GNN architecture
could be used for power systems of various sizes. We assume that the local
properties of the graphs in these systems are similar, leading to local
neighbourhoods with similar structures which can be represented using the same
embedding space size and the same number of GNN layers.

  * •

To evaluate the sample efficiency of the GNN model, we run multiple training
experiments on different sizes of training sets. Additionally, we assess the
scalability of the model by training it on various power system sizes and
evaluating its accuracy, training convergence properties, inference time, and
memory requirements.

  * •

As a side contribution, the proposed GNN model is tested in scenarios with
high measurement variances, using which we simulate phasor misalignments due
to communication delays, and the results are compared with linear WLS
solutions of SE.

##  II Linear State Estimation with PMUs

The SE algorithm has a goal of estimating the values of the state variables
$\mathbf{x}$, so that they are consistent with measurements, as well as the
power system model defined by its topology and parameters. The power system’s
topology is represented by a graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$,
where $\mathcal{N}={1,\dots,n}$ is the set of buses and
$\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ is the set of branches.

PMUs measure complex bus voltages and complex branch currents, in the form of
magnitude and phase angle [16, Sec. 5.6]. PMUs placed at a bus measure the bus
voltage phasor and current phasors along all branches incident to the bus
[17]. The state variables are given as $\mathbf{x}$ in rectangular
coordinates, and therefore consist of real and imaginary components of bus
voltages. The PMU measurements are transformed from polar to rectangular
coordinate system, since then the SE problem can be formulated using a system
of linear equations [15]. The solution to this sparse and noisy system can be
found by solving the linear WLS problem:

$\left(\mathbf{H}^{T}\mathbf{\Sigma}^{-1}\mathbf{H}\right)\mathbf{x}=\mathbf{H}^{T}\mathbf{\Sigma}^{-1}\mathbf{z},$
(1)

where the Jacobian matrix $\mathbf{H}\in\mathbb{R}^{m\times 2n}$ is defined
according to the partial first-order derivatives of the measurement functions,
and $m$ is the total number of linear equations. The observation error
covariance matrix is $\mathbf{\Sigma}\in\mathbb{R}^{m\times m}$, while the
vector $\mathbf{z}\in\mathbb{R}^{m}$ contains measurement values in
rectangular coordinate system. The aim of the WLS-based SE is to minimize the
sum of residuals between the measurements and the corresponding values that
are calculated using the measurement functions [1].

This approach has the disadvantage of requiring a transformation of
measurement errors (magnitude and angle errors) from polar to rectangular
coordinates, making them correlated, resulting in a non-diagonal covariance
matrix $\mathbf{\Sigma}$ and increased computational effort. To simplify the
calculation, the non-diagonal elements of $\mathbf{\Sigma}$ are often ignored,
which can impact the accuracy of the SE [17]. We can use the classical theory
of propagation of uncertainty to compute variances in rectangular coordinates
from variances in polar coordinates [18]. The solution to (1) obtained by
ignoring the non-diagonal elements of the covariance matrix $\mathbf{\Sigma}$
to avoid its computationally demanding inversion is referred to as the
approximative WLS SE solution.

In the rest of the paper, we will explore whether using a GNN model trained
with measurement values, variances, and covariances labelled with the exact
solutions of (1) leads to greater accuracy compared to the approximative WLS
SE, which ignores covariances. The GNN model, once trained, scales linearly
with respect to the number of power system buses, allowing for lower
computation time compared to both the approximate and exact solvers of (1).

##  III Methods

In this section, we introduce spatial GNNs on a high-level and describe how
can they be applied to the linear SE problem.

###  III-A Spatial Graph Neural Networks

Spatial GNNs are a type of machine learning models that process graph-
structured data by iteratively applying message passing to local subsets of
the graph. The goal of GNNs is to transform the inputs from each node and its
connections into a higher-dimensional space, creating a $s$-dimensional vector
$\mathbf{h}\in\mathbb{R}^{s}$ for each node. GNNs contain $K$ layers, with
each layer representing a single iteration $k$ of the message passing process.
Each GNN layer includes trainable functions, which are implemented as neural
networks, such as a message function, an aggregation function, and an update
function, as shown in Fig. 1.

MessageMessage…Message$\operatorname{Aggregate}$$\operatorname{Update}$$\mathbf{h_{2}}^{k-1}$$\mathbf{h_{1}}^{k-1}$$\mathbf{h_{n_{j}}}^{k-1}$$\mathbf{m_{1,j}}^{k-1}$$\mathbf{m_{2,j}}^{k-1}$$\mathbf{m_{{n_{j}},j}}^{k-1}$$\mathbf{m_{j}}^{k-1}$$\mathbf{h_{j}}^{k}$$\mathbf{{h_{j}}}^{k-1}$
Figure 1: A GNN layer, which represents a single message passing iteration,
includes multiple trainable functions, depicted as yellow rectangles. The
number of first-order neighbours of the node $j$ is denoted as $n_{j}$.

The message function calculates the message
$\mathbf{m}_{i,j}\in\mathbb{R}^{u}$ between two node embeddings, the
aggregation function combines the incoming messages in a specific way,
resulting in an aggregated message $\mathbf{m_{j}}\in\mathbb{R}^{u}$, and the
update function calculates the update to each node’s embedding. The message
passing process is repeated a fixed number of times, with the final node
embeddings passed through additional neural network layers to generate
predictions. GNNs are trained by optimizing their parameters using a variant
of gradient descent, with the loss function being a measure of the distance
between the ground-truth values and the predictions.

###  III-B State Estimation using Graph Neural Networks

The proposed GNN model is designed to be applied over a graph with a SE factor
graph topology [19], which consists of factor and variable nodes with edges
between them. The variable nodes are used to create a $s$-dimensional
embedding for the real and imaginary parts of the bus voltages, which are used
to generate state variable predictions. The factor nodes serve as inputs for
measurement values, variances, and covariances. Factor nodes do not generate
predictions, but they participate in the GNN message passing process to send
input data to their neighbouring variable nodes. To improve the model’s
representation of a node’s neighbourhood structure, we use binary index
encoding as input features for variable nodes. This encoding allows the GNN to
better capture relationships between nodes and reduces the number of input
neurons and trainable parameters, as well as training and inference time,
compared to the one-hot encoding used in [15]. The GNN model can be applied to
various types and quantities of measurements on both power system buses and
branches, and the addition or removal of measurements can be simulated by
adding or removing factor nodes. In contrast, applying a GNN to the bus-branch
power system model would require assigning a single input vector to each bus,
which can cause problems such as having to fill elements with zeros when not
all measurements are available and making the output sensitive to the order of
measurements in the input vector.

Connecting the variable nodes in the $2$-hop neighbourhood of the factor graph
topology significantly improves the model’s prediction quality in unobservable
scenarios [15]. This is because the graph remains connected even when
simulating the removal of factor nodes (e.g., measurement loss), which allows
messages to be propagated in the entire $K$-hop neighbourhood of the variable
node. This allows for the physical connection between power system buses to be
preserved when a factor node corresponding to a branch current measurement is
removed.

The proposed GNN for a heterogeneous graph has two types of layers: one for
factor nodes and one for variable nodes. These layers, denoted as
$\operatorname{Layer^{f}}$ and $\operatorname{Layer^{v}}$, have their own sets
of trainable parameters, which allow them to learn their message, aggregation,
and update functions separately. Different sets of trainable parameters are
used for variable-to-variable and factor-to-variable node messages. Both GNN
layers use two-layer feed-forward neural networks as message functions, single
layer neural networks as update functions, and the attention mechanism [7] in
the aggregation function. Then, a two-layer neural network
$\operatorname{Pred}$ is applied to the final node embeddings $\mathbf{h}^{K}$
of variable nodes only, to create state variable predictions. The loss
function is the mean-squared error (MSE) between the predictions and the
ground-truth values, calculated using variable nodes only. All trainable
parameters are updated via gradient descent and backpropagation over a mini-
batch of graphs. The high-level computational graph of the GNN architecture
specialized for heterogeneous augmented factor graphs is depicted in Figure 2.

$\operatorname{Layer^{f}}$$\operatorname{Layer^{v}}$$\mathbf{h_{f}}^{K-1}$$\mathbf{h_{v2}}^{K-1}$$\mathbf{h_{v}}^{K}$$\operatorname{Layer^{v}}$$\operatorname{Pred}$Loss$...$$...$$...$$...$$output$$label$
Figure 2: Proposed GNN architecture for heterogeneous augmented factor graphs.
Variable nodes are represented by circles and factor nodes are represented by
squares. The high-level computational graph begins with the loss function for
a variable node, and the layers that aggregate into different types of nodes
have distinct trainable parameters.

The proposed model uses an inference process that requires measurements from
the $K$-hop neighbourhood of each node, allowing for computational and
geographical distribution. Additionally, since the node degree in the SE
factor graph is limited, the computational complexity for the inference
process is constant. As a result, the overall GNN-based SE has a linear
computational complexity, making it efficient and scalable for large networks.

##  IV Numerical Results

In this section, we conduct numerical experiments to investigate the
scalability and sample efficiency of the proposed GNN approach. By varying the
power system and training set sizes, we are able to assess the model’s memory
requirements, prediction speed, and accuracy and compare them to those of
traditional SE approaches.

We use the IEEE 30-bus system, the IEEE 118-bus system, the IEEE 300-bus
system, and the ACTIVSg 2000-bus system [20], with measurements placed so that
measurement redundancy is maximal. For the purpose of sample efficiency
analysis, we create training sets containing 10, 100, 1000, and 10000 samples
for each of the mentioned power systems. Furthermore, we use validation and
test sets comprising 100 samples. These datasets are generated by solving the
power flow problem using randomly generated bus power injections and adding
Gaussian noise to obtain the measurement values. All the data samples were
labelled using the traditional SE solver. An instance of the GNN model is
trained on each of these datasets.

In contrast to our previous work, we use higher variance values of $5\times
10^{-1}$ to examine the performance of the GNN algorithm under conditions
where input measurement phasors are unsynchronized due to communication delays
[21]. While this is usually simulated by using variance values that increase
over time, as an extreme scenario we fix the measurement variances to a high
value.

In all the experiments, the node embedding size is set to $64$, and the
learning rate is $4\times 10^{-4}$. The minibatch size is $32$, and the number
of GNN layers is $4$. We use the ReLU activation function and a gradient
clipping value of $5\times 10^{-1}$. The optimizer is Adam, and we use mean
batch normalization.

###  IV-A Properties of Power System Augmented Factor Graphs

For all four test power systems, we create augmented factor graphs using the
methodology described in Section III-B. Fig. 3 illustrates how the properties
of the augmented factor graphs, such as average node degree, average path
length, average clustering coefficient, along with the system’s maximal
measurement redundancy, vary across different test power systems.

301183002000$0$$2$$4$$6$$8$$10$$12$$14$Number of busesPower systemgraph
propertiesRedundancyAvg. DegreeAvg. Path LengthAvg. Cluster Coeff. Figure 3:
Properties of augmented factor graphs along with the system’s measurement
redundancy for different test power systems, labelled with their corresponding
number of buses.

The average path length is a property that characterizes the global graph
structure, and it tends to increase as the size of the system grows. However,
as a design property of high-voltage networks, the other graph properties such
as the average node degree, average clustering coefficient, as well as maximal
measurement redundancy do not exhibit a clear trend of change with respect to
the size of the power system. This suggests that the structures of local,
$K$-hop neighbourhoods within the graph are similar across different power
systems, and that they contain similar factor-to-variable node ratio.
Consequently, it is reasonable to use the same GNN architecture (most
importantly, the number of GNN layers and the node embedding size) for all
test power systems, regardless of their size. In this way, the proposed model
achieves scalability, as it applies the same set of operations to the local,
$K$-hop neighbourhoods of augmented factor graphs of varying sizes without
having to adapt to each individual case.

###  IV-B Training Convergence Analysis

First, we analyse the training process for the IEEE 30-bus system with four
different sizes of the training set. As mentioned in III-B, the training loss
is a measure of the error between the predictions and the ground-truth values
for data samples used in the training process. The validation loss, on the
other hand, is a measure of the error between the predictions and the ground-
truth values on a separate validation set. In this analysis, we used a
validation set of 100 samples.

$1$$30$$60$$90$$120$$150$00.050.1EpochValidation
loss$n=10$$n=10^{2}$$n=10^{3}$$n=10^{4}$ Figure 4: Validation losses for
trainings on four different training set sizes.

The training losses for all the training processes converged smoothly, so we
do not plot them for the sake of clarity. Figure 4 shows the validation losses
for 150 epochs of training on four different training sets. For smaller
training sets, the validation loss decreases initially but then begins to
increase, which is a sign of overfitting. In these cases, a common practice in
machine learning is to select the model with the lowest validation loss value.
As it will be shown in IV-C, the separate test set results for models created
using small training sets are still satisfactory. As the number of samples in
the training set increases, the training process becomes more stable. This is
because the model has more data to learn from and is therefore less prone to
overfitting.

Approx. SE (baseline)GNN SE

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(a) IEEE 30

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(b) IEEE 118

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(c) IEEE 300

$10$$10^{2}$$10^{3}$$10^{4}$$0$$0.02$$0.04$Number of training samplesTest set
MSE(d) ACTIVSg 2000

Figure 5: Test set results for various power systems and training set sizes.

Next, in Table I, we present the training results for the other power systems
and training sets of various sizes. The numbers in the table represent the
number of epochs after which either the validation loss stopped changing or
began to increase. Similarly to the experiments on the IEEE 30-bus system, the
trainings on smaller training sets exhibited overfitting, while others
converged smoothly. For the former, the number in the table indicates the
epoch at which the validation loss reached its minimum and stopped improving.
For the latter, the number in the table represents the epoch when there were
five consecutive validation loss changes less than $10^{-5}$.

TABLE I: Epoch until validation loss minimum for various power systems and
training set sizes.

Power system | IEEE 118 | IEEE 300 | ACTIVSg 2000
---|---|---|---
10 samples | $61$ | $400$ | $166$
100 samples | $38$ | $84$ | $200$
1000 samples | $24$ | $82$ | $49$
10000 samples | $12$ | $30$ | $15$

Increasing the size of the training set generally results in a lower number of
epochs until the validation loss reaches its minimum. However, the epochs
until the validation loss reaches its minimum vary significantly between the
different power systems. This could be due to differences in the complexity of
the systems or the quality of the data used for training.

###  IV-C Accuracy Assessment

Fig. 5 reports the mean squared errors (MSEs) between the predictions and the
ground-truth values on 100-sample sized test sets for all trained models and
the approximate WLS SE. These results indicate that even the GNN models
trained on small datasets outperform the approximate WLS SE, except for the
models trained on the IEEE 30-bus system with 10 and 100 samples. These
results suggest that the quality of the GNN model’s predictions and the
generalization capabilities improve as the amount of training data increases,
and the models with the best results (highlighted in bold) have significantly
smaller MSEs compared to the approximate WLS SE. While we use randomly
generated training sets in this analysis, using carefully selected training
samples based on historical load consumption data could potentially lead to
even better results with small datasets.

###  IV-D Inference Time and Memory Requirements

The plot in Fig. 6 shows the ratio of execution times between WLS SE and GNN
SE inference as a function of the number of buses in the system. These times
are measured on a test set of 100 samples. As expected, the difference in
computational complexity between GNN, with its linear complexity, and WLS,
with more than quadratic complexity, becomes apparent as the number of buses
increases. From the results, it can be observed that GNN significantly
outperforms WLS in terms of inference time on larger power systems.

301183002000121064Number of busesInference time ratio Figure 6: A ratio of the
execution times for WLS SE and GNN SE inference on a test set of 100 samples,
as a function of the power system size.

The number of trainable parameters in the GNN model remains relatively
constant, as the number of power system buses increases. The number of input
neurons for variable node binary index encoding does grow logarithmically with
the number of variable nodes. However, this increase is relatively small
compared to the total number of GNN parameters111In fact, all the GNN models
we train have the same number of trainable parameters: 49921, which equates to
0.19 MB of memory.. This indicates that the GNN approach is scalable and
efficient, as the model’s complexity does not significantly increase with the
size of the power system being analysed.

##  V Conclusions

In this study, we focused on thoroughly testing a GNN-based state estimation
algorithm in scenarios with large variances, and examining its scalability and
sample efficiency. The results showed that the proposed approach provides good
results for large power systems, with lower prediction errors compared to the
approximative SE. The GNN model used in this approach is also fast and
maintains constant memory usage, regardless of the size of the scheme.
Additionally, the GNN was found to be an effective approximation method for
WLS SE even with a relatively small number of training samples, particularly
for larger power systems, indicating its sample efficiency. Given these
characteristics, the approach is worthy of further consideration for real-
world applications.

## References

  * [1] A. Monticelli, “Electric power system state estimation,” _Proceedings of the IEEE_ , vol. 88, no. 2, pp. 262–282, 2000.
  * [2] M. Göl and A. Abur, “A fast decoupled state estimator for systems measured by pmus,” _IEEE Trans. Power Syst._ , vol. 30, no. 5, pp. 2766–2771, 2015.
  * [3] G. N. Korres and N. M. Manousakis, “State estimation and observability analysis for phasor measurement unit measured systems,” _IET Gener. Transm. Dis._ , vol. 6, no. 9, pp. 902–913, September 2012.
  * [4] L. Zhang, G. Wang, and G. B. Giannakis, “Real-time power system state estimation and forecasting via deep unrolled neural networks,” _IEEE Trans. Signal Process._ , vol. 67, no. 15, pp. 4069–4077, 2019.
  * [5] A. S. Zamzam, X. Fu, and N. D. Sidiropoulos, “Data-driven learning-based optimization for distribution system state estimation,” _IEEE Trans. Power Syst._ , vol. 34, no. 6, pp. 4796–4805, 2019.
  * [6] W. L. Hamilton, “Graph representation learning,” _Synthesis Lectures on Artificial Intelligence and Machine Learning_ , vol. 14, no. 3, pp. 1–159, 2020\.
  * [7] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio, “Graph Attention Networks,” in _Proc. ICLR_ , 2018.
  * [8] O. Kundacina, M. Forcan, M. Cosovic, D. Raca, M. Dzaferagic, D. Miskovic, M. Maksimovic, and D. Vukobratovic, “Near real-time distributed state estimation via ai/ml-empowered 5g networks,” in _Proc. SmartGridComm_. IEEE, 2022, pp. 284–289.
  * [9] K. Chen, J. Hu, Y. Zhang, Z. Yu, and J. He, “Fault location in power distribution systems via deep graph convolutional networks,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 1, pp. 119–131, 2020.
  * [10] R. Zhang, W. Yao, Z. Shi, L. Zeng, Y. Tang, and J. Wen, “A graph attention networks-based model to distinguish the transient rotor angle instability and short-term voltage instability in power systems,” _Int. J. Electr. Power Energy Syst._ , vol. 137, p. 107783, 2022.
  * [11] S. Arastehfar, M. Matinkia, and M. R. Jabbarpour, “Short-term residential load forecasting using graph convolutional recurrent neural networks,” _Eng. Appl. Artif. Intell._ , vol. 116, p. 105358, 2022.
  * [12] V. Bolz, J. Rueß, and A. Zell, “Power flow approximation based on graph convolutional networks,” in _Proc. ICMLA_. IEEE, 2019, pp. 1679–1686.
  * [13] T. B. Lopez-Garcia and J. A. Domínguez-Navarro, “Power flow analysis via typed graph neural networks,” _Eng. Appl. Artif. Intell._ , vol. 117, p. 105567, 2023.
  * [14] Q. Yang, A. Sadeghi, and G. Wang, “Data-driven priors for robust psse via gauss-newton unrolled neural networks,” _IEEE J. Emerg. Sel. Top. Circuits Syst._ , vol. PP, pp. 1–1, 01 2022.
  * [15] O. Kundacina, M. Cosovic, and D. Vukobratovic, “State estimation in electric power systems leveraging graph neural networks,” in _Proc. PMAPS_ , 2022, pp. 1–6.
  * [16] J. De La Ree, V. Centeno, J. S. Thorp, and A. G. Phadke, “Synchronized phasor measurement applications in power systems,” _IEEE Trans. Smart Grid_ , vol. 1, no. 1, pp. 20–27, 2010.
  * [17] A. Gomez-Exposito, A. Abur, P. Rousseaux, A. de la Villa Jaen, and C. Gomez-Quiles, “On the use of PMUs in power system state estimation,” in _Proc. PSCC_ , 2011.
  * [18] M. Zhou, V. Centeno, J. Thorp, and A. Phadke, “An alternative for including phasor measurements in state estimators,” _IEEE Trans. Power Syst._ , vol. 21, no. 4, pp. 1930–1937, 2006.
  * [19] M. Cosovic and D. Vukobratovic, “Distributed Gauss–Newton method for state estimation using belief propagation,” _IEEE Trans. Power Syst._ , vol. 34, no. 1, pp. 648–658, 2019.
  * [20] A. B. Birchfield, T. Xu, K. M. Gegner, K. S. Shetye, and T. J. Overbye, “Grid structural characteristics as validation criteria for synthetic networks,” _IEEE Trans. Power Syst._ , vol. 32, no. 4, pp. 3258–3265, 2017.
  * [21] J. Zhao, A. Gómez-Expósito, M. Netto, L. Mili, A. Abur, V. Terzija, I. Kamwa, B. Pal, A. K. Singh, J. Qi _et al._ , “Power system dynamic state estimation: Motivations, definitions, methodologies, and future work,” _IEEE Trans. Power Syst._ , vol. 34, no. 4, pp. 3188–3198, 2019.