# Upper bounds for moments of zeta sums

Peng Gao  School of Mathematical Sciences, Beihang University, Beijing 100191,
China<EMAIL_ADDRESS>
###### Abstract.

We establish upper bounds for moments of zeta sums using results on shifted
moments of the Riemann zeta function under the Riemann hypothesis.

Mathematics Subject Classification (2010): 11M06

Keywords: zeta sums, moments

##  1\. Introduction

Character sums have been extensively studied in the literature as they have
many important applications in number theory. In [Harper23], A. J. Harper
studied sizes of the sums given by

$\sum_{n\leq x}n^{it}\;\;\;\;\;\text{and}\;\;\;\;\;\sum_{n\leq x}\chi(n),$

where $t\in\mathbb{R}$ and $\chi(n)$ is a non-principal Dirichlet character
modulo a large prime $r$. Following the notation in [Harper23], we shall refer
the first sum above as a zeta sum.

Building on his work concerning moments of random multiplicative functions,
Harper [Harper23] showed that the low moments of zeta sums (and also character
sums) have “better than squareroot cancellation”. More precisely, he proved
that uniformly for $1\leq x\leq T$ and $0\leq k\leq 1$,

$\frac{1}{T}\int_{0}^{T}\Big{|}\sum_{n\leq
x}n^{it}\Big{|}^{2k}\mathrm{d}t\ll\bigg{(}\frac{x}{1+(1-k)\sqrt{\log\log(10L_{T})}}\bigg{)}^{k},$

where $L_{T}=\min\\{x,T/x\\}.$

In [Szab], B. Szabó obtained sharp upper bounds on shifted moments of
Dirichlet $L$-function at points on the critical line and then applied the
results to show under the generalized Riemann hypothesis (GRH) that for a
fixed real number $k>2$ and a large integer $q$, we have for $2\leq Y\leq
q^{1/2}$,

(1.1) $\displaystyle\sum_{\chi\in X_{q}^{*}}\bigg{|}\sum_{n\leq
Y}\chi(n)\bigg{|}^{2k}\ll_{k}\phi(q)Y^{k}(\log Y)^{(k-1)^{2}},$

where $X_{q}^{*}$ denotes the set of primitive Dirichlet characters modulo $q$
and $\phi$ denotes Euler’s totient function. A similar result is given in
[G&Zhao2024] for moments of quadratic Dirichlet character sums under GRH.

We note that the zeta sums behave very much like character sums. In fact,
other than periodicity, the function $n\mapsto n^{-it}$ for a fixed
$t\in\mathbb{R}$ is totally multiplicative and is unimodular. Thus, one
expects to establish results analogous to (1.1) for moments of zeta sums and
it is the aim of this paper to achieve this. For this, we define for real
numbers $m,T,Y>0$,

$\displaystyle S_{m}(T,Y):=\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\bigg{|}^{2m}dt.$

We are interested in bounding $S_{m}(T,Y)$ from the above. We first observe
that as pointed out in [Harper23] that using [Ivic, Lemma 1.2] that when $t$
is large and $x\geq t$, we have $\sum_{x<n\leq
2x}n^{it}=\frac{(2x)^{1+it}-x^{1+it}}{1+it}+O(1)\ll x/t$. Moreover, note that
by [Montgomery94, Chap. 7,(34)], we have $\sum_{n\leq x}n^{it}\ll t^{1/2}\log
t$ when $x\leq t$. As the term $x/t$ dominates $t^{1/2}\log t$ when $x\geq
t^{3/2}\log t$, we deduce that when $T$ is large enough and $Y\geq T^{3/2}\log
T$, we have for any real $m>0$,

$\displaystyle S_{m}(T,Y)\ll_{k}T^{1-2m}Y^{2m}.$

We may therefore focus on the case $Y<T^{3/2}\log T$. In fact, we shall assume
that $Y\leq(1-\varepsilon)T$ for any $\varepsilon>0$ throughout the paper as
this is often the most interesting case regarding character sums. For this
case, we establish the following result concerning the size of $S_{m}(T,Y)$
under the Riemann hypothesis (RH).

######  Theorem 1.1.

With the notation as above and assume the truth of RH. For any real number
$m>2$, large real numbers $T,Y$ such that $Y\leq(1-\varepsilon)T$ for any
$\varepsilon>0$, we have

(1.2) $\displaystyle S_{m}(T,Y)\ll TY^{m}(\log T)^{(m-1)^{2}}.$

We note that by Hölder’s inequality, we have for any real number $n>1$,

$\displaystyle S_{m}(T,Y)\ll T^{1-1/n}(S_{mn}(T,Y))^{1/n}.$

The above together with Theorem 1.1 then implies that $S_{m}(T,Y)\ll
TY^{m}(\log T)^{O(1)}$ for any $m>0$, upon choosing $n$ large enough. We
remark here that it is shown in [Harper23] that one has $S_{m}(T,Y)\ll
T^{m+1}$, so that our result above improves upon this when $Y$ is slightly
smaller than $T$.

Our proof of Theorem 1.1 follows the approaches in [Szab]. A key ingredient
used in the proof is a result of M. J. Curran [Curran] on shifted moments of
the Riemann zeta function $\zeta(s)$.

##  2\. Preliminaries

In this section, we include some results concerning shifted moments of the
Riemann zeta function. The first one is quoted from [Curran, Theorem 1.1].

######  Proposition 2.1.

With the notation as above and assume the truth of RH. Let $k\geq 1$ be a
fixed integer and $a_{1},\ldots,a_{k}$ be fixed non-negative real numbers. Let
$T$ be a large real number and let ${\bf b}=(b_{1},\ldots,b_{k})$ be a real
$k$-tuple with $|b_{j}|\leq(1-\varepsilon)T$ for a fixed $\varepsilon>0$. Then

$\displaystyle\begin{split}\int_{T}^{2T}\prod_{j=1}^{k}|\zeta(\tfrac{1}{2}+i(t+b_{k}))|^{a_{k}}dt\ll
T(\log T)^{(a_{1}^{2}+\cdots+a_{k}^{2})/4}\prod_{1\leq j<l\leq
l}|\zeta(1+i(b_{j}-b_{l})+1/\log T)|^{a_{j}a_{l}/2}.\end{split}$

Here the implied constant depends on $k$ and the $a_{j}$ but not on $T$ or the
$b_{j}$.

We remark here that [Curran, Theorem 1.1] is stated for $|b_{j}|\leq T/2$ but
an inspection of the proof indicates that it continues to hold for
$|b_{j}|\leq(1-\varepsilon)T$ with any $\varepsilon>0$. We also note that

$\Big{|}\zeta(1+1/\log
T+i\alpha)\Big{|}=\Big{|}\sum^{\infty}_{n=1}n^{-(1+1/\log
T+i\alpha)}\Big{|}\leq\Big{|}\sum^{\infty}_{n=1}n^{-(1+1/\log
T)}\Big{|}=\Big{|}\zeta(1+1/\log T)\Big{|}\ll\log T,$

where the last estimation above follows from [MVa1, Corollary 1.17]. Also by
[MVa1, Corollary 1.17], we see that for $\frac{1}{\log T}\leq|\alpha|\leq 10$,
we have

$|\zeta(1+1/\log T+i\alpha)|=\frac{1}{|1/\log
T+i\alpha|}+O(1)\ll\frac{1}{|\alpha|}.$

Moreover, by [MVa1, Corollary 13.16], we see that for $10\leq|\alpha|\leq
e^{T}$, we have under the RH that

$\log|\zeta(1+1/\log T+i\alpha)|\leq\log\log\log|\alpha|+O(1).$

Based on these observations, for $T$ be given as in Proposition 2.1, we now
introduce the function $g:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}$ defined by

(2.1) $\displaystyle\begin{split}g(x)=\begin{cases}\log T&\text{if
}x\leq\frac{1}{\log T}\text{ or }x\geq e^{T},\\\ \frac{1}{x}&\text{if
}\frac{1}{\log T}\leq x\leq 10,\\\ \log\log x&\text{if }10\leq x\leq
e^{T}.\end{cases}\end{split}$

The above discussions together with Proposition 2.1 allows us to derive the
following simplified version on shifted moments of the Riemann zeta function.

######  Corollary 2.2.

With the notation as above and assume the truth of RH. Let $k\geq 1$ be a
fixed integer and $a_{1},\ldots,a_{k}$ be fixed non-negative real numbers. Let
$T$ be a large real number and let ${\bf b}=(b_{1},\ldots,b_{k})$ be a real
$k$-tuple with $|b_{j}|\leq(1-\varepsilon)T$ for a fixed $\varepsilon>0$. Then

$\displaystyle\begin{split}\int_{T}^{2T}\prod_{j=1}^{k}|\zeta(\tfrac{1}{2}+i(t+b_{k}))|^{a_{k}}dt\ll
T(\log T)^{(a_{1}^{2}+\cdots+a_{k}^{2})/4}\prod_{1\leq j<l\leq
l}g(|b_{j}-b_{l}|)^{a_{j}a_{l}/2}.\end{split}$

Here the implied constant depends on $k,\varepsilon$ and the $a_{j}$ but not
on $T$ or the $b_{j}$.

We also note the following upper bounds on moments of the Riemann zeta
function, which can be obtained by modifying the proof of [Curran, Theorem
1.1].

######  Lemma 2.3.

With the notation as above and assume the truth of RH. Let $k\geq 1$ be a
fixed integer and $a_{1},\ldots,a_{k}$ be fixed non-negative real numbers. Let
$T$ be a large real number and let ${\bf b}=(b_{1},\ldots,b_{k})$ be a real
$k$-tuple with $|b_{j}|\leq(1-\varepsilon)T$ for any fixed $\varepsilon>0$.
Then for large real number $T$ and $\sigma\geq 1/2$,

$\displaystyle\int_{T}^{2T}\prod_{j=1}^{k}|\zeta(\sigma+i(t+b_{k}))|^{a_{k}}dt\ll$
$\displaystyle T(\log T)^{O(1)}.$

We end this section by including an estimation for an average of the moments
of the Riemann zeta function.

######  Proposition 2.4.

With the notation as above and assume the truth of RH. We have for any real
numbers $m\geq 2$, $10\leq E\leq(1-\varepsilon)T$ with $\varepsilon>0$ being
fixed,

(2.2)
$\displaystyle\begin{split}&\int^{2T}_{T}\bigg{(}\int_{0}^{E}|\zeta(1/2+i(\pm
s+t))|ds\bigg{)}^{2m}dt\\\ \ll&T\big{(}(\log T)^{(m-1)^{2}}E^{3}(\log\log
E)^{O(1)}+(\log T)^{m^{2}-3m+3}E^{2m}(\log\log E)^{O(1)}(\log\log
T)^{O(1)}\big{)}.\end{split}$

###### Proof.

Our proof follows closely that of [Szab, Proposition 3]. Without loss of
generality, we prove (2.2) only for the case where the sign $\pm$ in front of
$s$ is $+$ in what follows. We have by symmetry that for each fixed $t$ and
any fixed integer $k\geq 1$,

(2.3)
$\displaystyle\bigg{(}\int_{0}^{E}|\zeta(1/2+i(s+t))|ds\bigg{)}^{2m}\ll\int_{[0,E]^{k}}\prod_{a=1}^{k}|\zeta(1/2+i(t_{a}+t))|\cdot\bigg{(}\int_{\mathcal{D}}|\zeta(1/2+i(u+t))|du\bigg{)}^{2m-k}d\mathbf{t},$

where
$\mathcal{D}=\mathcal{D}(t_{1},\ldots,t_{k})=\\{u\in[0,E]:|t_{1}-u|\leq|t_{2}-u|\leq\ldots\leq|t_{k}-u|\\}$.

We let $\mathcal{B}_{1}=\big{[}-\frac{1}{\log T},\frac{1}{\log T}\big{]}$ and
$\mathcal{B}_{j}=\big{[}-\frac{e^{j-1}}{\log T},-\frac{e^{j-2}}{\log
T}\big{]}\cup\big{[}\frac{e^{j-2}}{\log T},\frac{e^{j-1}}{\log T}\big{]}$ for
$2\leq j<\lfloor\log\log T\rfloor+10:=K$. We further denote
$\mathcal{B}_{K}=[-E,E]\setminus\bigcup_{1\leq j<K}\mathcal{B}_{j}$.

Observe that for any $t_{1}\in[0,E]$, we have $\mathcal{D}\subset[0,E]\subset
t_{1}+[-E,E]\subset\bigcup_{1\leq j\leq K}t_{1}+\mathcal{B}_{j}$. Thus if we
denote $\mathcal{A}_{j}=\mathcal{B}_{j}\cap(-t_{1}+\mathcal{D})$, then
$(t_{1}+\mathcal{A}_{j})_{1\leq j\leq K}$ form a partition of $\mathcal{D}$.
We apply Hölder’s inequality twice to deduce that for $2m\geq k+1$,

(2.4)
$\displaystyle\begin{split}&\bigg{(}\int_{\mathcal{D}}|\zeta(1/2+i(u+t))|du\bigg{)}^{2m-k}\\\
\leq&\bigg{(}\sum_{1\leq j\leq K}\frac{1}{j}\cdot
j\int_{t_{1}+\mathcal{A}_{j}}|\zeta(1/2+i(u+t))|du\bigg{)}^{2m-k}\\\
\leq&\bigg{(}\sum_{1\leq j\leq
K}j^{2m-k}\bigg{(}\int_{t_{1}+\mathcal{A}_{j}}\big{|}\zeta(1/2+i(u+t))\big{|}du\bigg{)}^{2m-k}\bigg{)}\bigg{(}\sum_{1\leq
j\leq K}j^{-(2m-k)/(2m-k-1)}\bigg{)}^{2m-k-1}\\\ \ll&\sum_{1\leq j\leq
K}j^{2m-k}\bigg{(}\int_{t_{1}+\mathcal{A}_{j}}|\zeta(1/2+i(u+t))|du\bigg{)}^{2m-k}\\\
\leq&\sum_{1\leq j\leq
K}j^{2m-k}|\mathcal{B}_{j}|^{2m-k-1}\int_{t_{1}+\mathcal{A}_{j}}|\zeta(1/2+i(u+t))|^{2m-k}du.\end{split}$

We denote for $\mathbf{t}=(t_{1},\ldots,t_{k})$,

$\zeta(\mathbf{t},u)=\int^{2T}_{T}\prod_{a=1}^{k}|\zeta(1/2+i(t_{a}+t))|\cdot|\zeta(1/2+i(u+t))|^{2m-k}dt.$

We then deduce from (2.3) and (2.4) that

(2.5)
$\displaystyle\begin{split}\int^{2T}_{T}\bigg{(}\int_{0}^{E}|\zeta(1/2+i(s+t))|ds\bigg{)}^{2m}dt\ll&\sum_{1\leq
l_{0}\leq
K}l_{0}^{2m-k}|\mathcal{B}_{l_{0}}|^{2m-k-1}\int_{[0,E]^{k}}\int_{t_{1}+\mathcal{A}_{l_{0}}}\zeta(\mathbf{t},u)dud\mathbf{t}\\\
\ll&\sum_{1\leq l_{0},l_{1},\ldots l_{k-1}\leq
K}l_{0}^{2m-k}|\mathcal{B}_{l_{0}}|^{2m-k-1}\int_{\mathcal{C}_{l_{0},l_{1},\cdots,l_{k-1}}}\zeta(\mathbf{t},u)dud\mathbf{t},\end{split}$

where

$\mathcal{C}_{l_{0},l_{1},\cdots,l_{k-1}}=\\{(t_{1},\ldots,t_{k},u)\in[0,E]^{k+1}:u\in
t_{1}+\mathcal{A}_{l_{0}},\,|t_{i+1}-u|-|t_{i}-u|\in\mathcal{B}_{l_{i}},\
1\leq i\leq k-1\\}.$

We now distinguish two cases in the last summation of (2.5) according to the
size of $l_{0}$.

Case 1: $l_{0}<K$. First note that for any fixed $u$, $t_{1}$ is in a fixed
region of size $\ll\frac{e^{l_{0}}}{\log T}$. For fixed $u$ and $t_{1}$,
$t_{2}$ is in a fixed region of size $\ll E\frac{e^{l_{1}}}{\log T}$ as
$|t_{2}-u|\in|t_{1}-u|+\mathcal{B}_{l_{1}}$). Similar considerations then
imply that the volume of the region $\mathcal{C}_{l_{0},l_{1},\cdots,l_{k-1}}$
is $\ll E^{k}\frac{e^{l_{0}+l_{1}+\cdots+l_{k-1}}}{(\log T)^{k}}$. Also, by
the definition of $\mathcal{C}_{l_{0},l_{1},\cdots,l_{k-1}}$ we have
$\frac{e^{l_{0}}}{\log q}\ll|t_{1}-u|\ll E=(\log T)^{O(1)}$ so that
$g(|t_{1}-u|)\ll\frac{\log T}{e^{l_{0}}}\log\log E$, where $g$ is the function
defined in (2.1). We deduce from the definition of $\mathcal{A}_{j}$ that
$|t_{2}-u|\geq|t_{1}-u|$, so that
$E\gg|t_{2}-u|=|t_{1}-u|+(|t_{2}-u|-|t_{1}-u|)\gg\frac{e^{l_{0}}}{\log
T}+\frac{e^{l_{1}}}{\log T}$, which implies that $g(|t_{2}-u|)\ll\frac{\log
T}{e^{\max(l_{0},l_{1})}}\log\log E$. Similarly, we have
$g(|t_{i}-u|)\ll\frac{\log T}{e^{\max(l_{0},l_{1},\ldots,l_{i-1})}}\log\log E$
for any $1\leq i\leq k$. Moreover, we have
$\sum^{j-1}_{s=i}(|t_{s+1}-u|-|t_{s}-u|)\leq|t_{j}-t_{i}|$ for any $1\leq
i<j\leq k$, so that we have $g(|t_{j}-t_{i}|)\ll\frac{\log
T}{e^{\max(l_{i},\ldots,l_{j-1})}}\log\log E$. We then deduce from Theorem 2.1
that for $(t_{1},\ldots,t_{k},u)\in\mathcal{C}_{l_{0},l_{1},\cdots,l_{k-1}}$,

$\displaystyle\zeta(\mathbf{t},u)$ $\displaystyle\ll$ $\displaystyle T(\log
T)^{((2m-k)^{2}+k)/4}(\log\log E)^{O(1)}\bigg{(}\prod^{k-1}_{i=0}\frac{\log
T}{e^{\max(l_{0},l_{1},\ldots,l_{i})}}\bigg{)}^{(2m-k)/2}\bigg{(}\prod^{k-1}_{i=1}\prod^{k}_{j=i+1}\frac{\log
T}{e^{\max(l_{i},\ldots,l_{j-1})}}\bigg{)}^{1/2}$ $\displaystyle=$
$\displaystyle T(\log T)^{m^{2}}(\log\log
E)^{O(1)}\exp\Big{(}-\frac{2m-k}{2}\sum^{k-1}_{i=0}\max(l_{0},l_{1},\ldots,l_{i})-\frac{1}{2}\sum^{k-1}_{i=1}\sum^{k}_{j=i+1}\max(l_{i},\ldots,l_{j-1})\Big{)}.$

Here, we adopt the convention throughout the paper that any empty product is
defined to be $1$ and any empty sum is defined to be $0$. Observe that we have
$|\mathcal{B}_{l_{0}}|\ll\frac{e^{l_{0}}}{\log T}$, so that

(2.6) $\displaystyle\begin{split}&\sum_{\begin{subarray}{c}1\leq l_{0}<K\\\
1\leq l_{1},\ldots l_{k-1}\leq
K\end{subarray}}l_{0}^{2m-k}|\mathcal{B}_{l_{0}}|^{2m-k-1}\int_{\mathcal{C}_{l_{0},l_{1},\cdots,l_{k-1}}}\zeta(\mathbf{t},u)dud\mathbf{t}\\\
\ll&T(\log T)^{(m-1)^{2}}E^{k}(\log\log E)^{O(1)}\\\
&\cdot\sum_{\begin{subarray}{c}1\leq l_{0}<K\\\ 1\leq l_{1},\ldots l_{k-1}\leq
K\end{subarray}}l_{0}^{2m-k}\exp\Big{(}(2m-k-1)l_{0}+\sum^{k-1}_{i=0}l_{i}-\frac{2m-k}{2}\sum^{k-1}_{i=0}\max(l_{0},l_{1},\ldots,l_{i})-\frac{1}{2}\sum^{k-1}_{i=1}\sum^{k}_{j=i+1}\max(l_{i},\ldots,l_{j-1})\Big{)}\\\
=&T(\log T)^{(m-1)^{2}}E^{k}(\log\log E)^{O(1)}\\\
&\cdot\sum_{\begin{subarray}{c}1\leq l_{0}<K\\\ 1\leq l_{1},\ldots l_{k-1}\leq
K\end{subarray}}l_{0}^{2m-k}\exp\Big{(}\frac{2m-k}{2}l_{0}+\frac{1}{2}\sum^{k-1}_{i=1}l_{i}-\frac{2m-k}{2}\sum^{k-1}_{i=1}\max(l_{0},l_{1},\ldots,l_{i})-\frac{1}{2}\sum^{k-1}_{i=1}\sum^{k}_{j=i+2}\max(l_{i},\ldots,l_{j-1})\Big{)}.\end{split}$

We now set $k=4$ to see that in this case deduce from the above that

$\displaystyle\begin{split}&\frac{2m-k}{2}l_{0}+\frac{1}{2}\sum^{k-1}_{i=1}l_{i}-\frac{2m-k}{2}\sum^{k-1}_{i=1}\max(l_{0},l_{1},\ldots,l_{i})-\frac{1}{2}\sum^{k-1}_{i=1}\sum^{k}_{j=i+2}\max(l_{i},\ldots,l_{j-1})\\\
\ll&-(2m-4)\max(l_{0},\ldots,l_{k-1}).\end{split}$

We deduce from (2.6) and the above that

(2.7) $\displaystyle\begin{split}&\sum_{\begin{subarray}{c}1\leq l_{0}<K\\\
1\leq l_{1},\ldots l_{k-1}\leq
K\end{subarray}}l_{0}^{2m-k}|\mathcal{B}_{l_{0}}|^{2m-k-1}\int_{\mathcal{C}_{l_{0},l_{1},\cdots,l_{k-1}}}\zeta(\mathbf{t},u)dud\mathbf{t}\\\
\ll&T(\log T)^{(m-1)^{2}}E^{k}(\log\log
E)^{O(1)}\sum_{\begin{subarray}{c}1\leq l_{0}<K\\\ 1\leq l_{1},\ldots
l_{k-1}\leq
K\end{subarray}}l_{0}^{2m-k}\exp\Big{(}-(2m-4)\max(l_{0},\ldots,l_{k-1})\Big{)}\\\
\ll&T(\log T)^{(m-1)^{2}}E^{k}(\log\log E)^{O(1)},\end{split}$

where the last estimation above follows by noting that we have $m>2$.

Case 2 $l_{0}=K$. The volume of the region
$\mathcal{C}_{K,l_{1},\cdots,l_{k-1}}$ is $\ll
E^{k+1}\frac{e^{l_{1}+\cdots+l_{k-1}}}{(\log T)^{k-1}}$. For each $1\leq i\leq
k$, we have $g(|t_{i}-u|)\ll\log\log E$. Also, similar to Case 1, we have
$g(|t_{j}-t_{i}|)\ll\frac{\log T}{e^{\max(l_{i},\ldots,l_{j-1})}}$ for $1\leq
i\leq k$.

$\displaystyle\zeta(\mathbf{t},u)$ $\displaystyle\ll$ $\displaystyle T(\log
T)^{((2m-k)^{2}+k)/4}(\log\log
E)^{O(1)}\bigg{(}\prod^{k-1}_{i=1}\prod^{k}_{j=i+1}\frac{\log
T}{e^{\max(l_{i},\ldots,l_{j-1})}}\bigg{)}^{1/2}$ $\displaystyle=$
$\displaystyle T(\log T)^{((2m-k)^{2}+k^{2})/4}(\log\log
E)^{O(1)}\exp\Big{(}-\frac{1}{2}\sum^{k-1}_{i=1}\sum^{k}_{j=i+1}\max(l_{i},\ldots,l_{j-1})\Big{)}.$

As $|\mathcal{B}_{K}|\ll E$, we see that

(2.8) $\displaystyle\begin{split}&\sum_{\begin{subarray}{c}1\leq l_{1},\ldots
l_{k-1}\leq
K\end{subarray}}K^{2m-k}|\mathcal{B}_{K}|^{2m-k-1}\int_{\mathcal{C}_{K,l_{1},\cdots,l_{k-1}}}\zeta(\mathbf{t},u)dud\mathbf{t}\\\
\ll&T(\log T)^{((2m-k)^{2}+k^{2})/4-k+1}E^{2m}(\log\log E)^{O(1)}(\log\log
T)^{O(1)}\sum_{\begin{subarray}{c}1\leq l_{1},\ldots l_{k-1}\leq
K\end{subarray}}\exp\Big{(}\sum^{k-1}_{i=1}l_{i}-\frac{1}{2}\sum^{k-1}_{i=1}\sum^{k}_{j=i+1}\max(l_{i},\ldots,l_{j-1})\Big{)}.\end{split}$

We now set $k=4$ to see that in this case deduce from the above that

$\displaystyle\begin{split}&\sum^{k-1}_{i=1}l_{i}-\frac{1}{2}\sum^{k-1}_{i=1}\sum^{k}_{j=i+1}\max(l_{i},\ldots,l_{j-1})\ll
l_{1}/2.\end{split}$

We deduce from (2.8) and the above that

(2.9) $\displaystyle\begin{split}&\sum_{\begin{subarray}{c}1\leq l_{1},\ldots
l_{k-1}\leq
K\end{subarray}}K^{2m-k}|\mathcal{B}_{K}|^{2m-k-1}\int_{\mathcal{C}_{K,l_{1},\cdots,l_{k-1}}}\zeta(\mathbf{t},u)dud\mathbf{t}\\\
\ll&T(\log T)^{((2m-3)^{2}+3^{2})/4-2}E^{2m}(\log\log E)^{O(1)}(\log\log
T)^{O(1)}\sum_{\begin{subarray}{c}1\leq l_{1},\ldots l_{k-1}\leq
K\end{subarray}}\exp\Big{(}l_{1}/2\Big{)}\\\ \ll&T(\log
T)^{m^{2}-3m+3}E^{2m}(\log\log E)^{O(1)}(\log\log T)^{O(1)}.\end{split}$

We now deduce the estimation in (2.2) using (2.7) and (2.9). This completes
the proof of the proposition. ∎

##  3\. Proof of Theorem 1.1

###  3.1. Initial treatments

As we explained in the paragraph below Theorem 1.1, it suffices to establish
(1.2). We let $\Phi_{U}(tx)$ be a non-negative smooth function supported on
$(0,1)$, satisfying $\Phi_{U}(x)=1$ for $t\in(1/U,1-1/U)$ with $U$ a parameter
to be chosen later and such that $\Phi^{(j)}_{U}(x)\ll_{j}U^{j}$ for all
integers $j\geq 0$. We denote the Mellin transform of $\Phi_{U}$ by
$\widehat{\Phi}_{U}$ and we observe that repeated integration by parts gives
that, for any integer $i\geq 1$ and $\Re(s)\geq 1/2$,

(3.1) $\displaystyle\widehat{\Phi}_{U}(s)\ll U^{i-1}(1+|s|)^{-i}.$

We insert the function $\Phi_{U}(\frac{n}{Y})$ into the definition of
$S_{m}(T,Y)$ and apply the triangle inequality to obtain that

(3.2)
$\displaystyle\begin{split}S_{m}(T,Y)\leq&\int^{2T}_{T}\Big{|}\sum_{n}n^{-it}\Phi_{U}(\frac{n}{Y})\Big{|}^{2m}dt+\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{2m}dt.\end{split}$

We further apply the Mellin inversion to obtain that

$\displaystyle\begin{split}\int^{2T}_{T}\Big{|}\sum_{n}n^{-it}\Phi_{U}(\frac{n}{Y})\Big{|}^{2m}dt=&\int^{2T}_{T}\Big{|}\int\limits_{(2)}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt.\end{split}$

Observe that by [iwakow, Corollary 5.20] that under RH, we have for
$\Re(s)\geq 1/2$ and any $\varepsilon>0$,

(3.3) $\displaystyle\zeta(s)\ll|s|^{\varepsilon}.$

The bounds in (3.1) and (3.3) allow us to shift the line of integration in
(3.2) to $\Re(s)=1/2$ to obtain that

(3.4)
$\displaystyle\begin{split}\int^{2T}_{T}\Big{|}\sum_{n}n^{-it}\Phi_{U}(\frac{n}{Y})\Big{|}^{2m}dt=&\int^{2T}_{T}\Big{|}\int\limits_{(1/2)}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt.\end{split}$

We split the last integral above according to whether $|\Im(s)|\leq(\log
T)^{D}$ or not for some $D>0$ to be specified later, obtaining

$\displaystyle\begin{split}&\int^{2T}_{T}\Big{|}\sum_{n}n^{-it}\Phi_{U}(\frac{n}{Y})\Big{|}^{2m}dt\\\
\ll&\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|\leq(\log
T)^{D}\end{subarray}}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt+\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt.\end{split}$

We now set $U=(\log T)^{C}$ to deduce from (3.2), (3.4) and the above that

(3.5)
$\displaystyle\begin{split}S_{m}(T,Y)\ll&\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|\leq(\log
T)^{D}\end{subarray}}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt+\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt\\\
&+\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{2m}dt\\\
\ll&Y^{m}\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|s|\leq(\log
T)^{D}\end{subarray}}\Big{|}\zeta(1/2+i(s+t))\Big{|}\frac{1}{1+|s|}ds\Big{|}^{2m}dt+\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt\\\
&+\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{2m}dt.\end{split}$

It follows from the above that in order to establish Theorem 1.1, it remains
to prove the following results.

######  Lemma 3.2.

With the notation as above and assume the truth of RH. We have for $D\gg C$
large enough and any real number $m>2$,

(3.6) $\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt\ll
TY^{m}.$

######  Lemma 3.3.

With the notation as above and assume the truth of RH. We have for $D$ large
enough and any real number $m>2$,

(3.7) $\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|s|\leq(\log
T)^{D}\end{subarray}}\Big{|}\zeta(1/2+i(s+t))\Big{|}\frac{1}{1+|s|}ds\Big{|}^{2m}dt\ll
T(\log T)^{(m-1)^{2}}.$

######  Lemma 3.4.

With the notation as above and assume the truth of RH. We have for $C$ large
enough and any real number $m>2$,

(3.8) $\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{2m}dt\ll TY^{m}.$

###  3.5. Proof of Lemma 3.2

We apply (3.1) and Hölder’s inequality to deduce that, as $m\geq 1/2$,

(3.9)
$\displaystyle\begin{split}&\int^{2T}_{T}\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\zeta(s+it)Y^{s}\widehat{\Phi}_{U}(s)ds\Big{|}^{2m}dt\\\
\ll&\int^{2T}_{T}\Big{(}\Big{(}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\Big{|}\widehat{\Phi}_{U}(s)\Big{|}|\mathrm{d}s|\Big{)}^{2m-1}\cdot\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\Big{|}\zeta(s+it)Y^{s}\Big{|}^{2m}\Big{|}\widehat{\Phi}_{U}(s)\Big{|}|\mathrm{d}s|\Big{)}dt\\\
\ll&Y^{m}\Big{(}\int\limits_{\begin{subarray}{c}(1/2)\\\ |\Im(s)|>(\log
T)^{D}\end{subarray}}\Big{|}\widehat{\Phi}_{U}(s)\Big{|}|\mathrm{d}s|\Big{)}^{2m-1}\cdot\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(s+it)\Big{|}^{2m}\Big{|}\widehat{\Phi}_{U}(s)\Big{|}dt|\mathrm{d}s|.\end{split}$

Now we note that by (3.1),

(3.10) $\displaystyle\begin{split}\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\Big{|}\widehat{\Phi}_{U}(s)\Big{|}|\mathrm{d}s|\ll&\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\frac{U}{1+|s|^{2}}|\mathrm{d}s|\ll_{D}\frac{U}{(\log
T)^{D}}.\end{split}$

We also apply (3.3) to see that

(3.11) $\displaystyle\begin{split}&\int\limits_{\begin{subarray}{c}(1/2)\\\
|\Im(s)|>(\log
T)^{D}\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(s+it)\Big{|}^{2m}\cdot\Big{|}\widehat{\Phi}_{U}(s)\Big{|}dt|\mathrm{d}s|\\\
\ll&\int\limits_{\begin{subarray}{c}(1/2)\\\ |s|>(\log
T)^{D}\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(1/2+i(s+t))\Big{|}^{2m}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|\\\
\ll&\int\limits_{\begin{subarray}{c}(1/2)\\\ (\log T)^{D}<|s|\leq
5T\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(1/2+i(s+t))\Big{|}^{2m}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|+\int\limits_{\begin{subarray}{c}(1/2)\\\
|s|>5T\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(1/2+i(s+t))\Big{|}^{2m}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|\\\
\ll&\int\limits_{\begin{subarray}{c}(1/2)\\\ (\log T)^{D}<|s|\leq
5T\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(1/2+i(s+t))\Big{|}^{2m}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|+\int\limits_{\begin{subarray}{c}(1/2)\\\
|s|>5T\end{subarray}}\int^{2T}_{T}\Big{|}s+t\Big{|}^{\varepsilon}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|\\\
\ll&\int\limits_{\begin{subarray}{c}(1/2)\\\ (\log T)^{D}<|s|\leq
5T\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(1/2+i(s+t))\Big{|}^{2m}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|+\int\limits_{\begin{subarray}{c}(1/2)\\\
|s|>5T\end{subarray}}\int^{2T}_{T}\Big{|}s\Big{|}^{\varepsilon}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|\\\
\ll&\int\limits_{\begin{subarray}{c}(1/2)\\\ (\log T)^{D}<|s|\leq
5T\end{subarray}}\int^{2T}_{T}\Big{|}\zeta(1/2+i(s+t))\Big{|}^{2m}\cdot\frac{U}{1+|s|^{2}}dt|\mathrm{d}s|+UT^{\varepsilon}\\\
\ll&\int\limits_{\begin{subarray}{c}(1/2)\\\ (\log T)^{D}<|s|\leq
5T\end{subarray}}\int^{10T}_{-10T}\Big{|}\zeta(1/2+it)\Big{|}^{2m}\cdot\frac{U}{1+|s|^{2}}|\mathrm{d}t||\mathrm{d}s|+UT^{\varepsilon}\\\
\ll&TU(\log T)^{O(1)}(\log T)^{-D}+UT^{\varepsilon},\end{split}$

where the last estimation above follows from [Sound2009, Corollary B], which
asserts that

$\displaystyle\begin{split}\int^{10T}_{-10T}\Big{|}\zeta(1/2+it)\Big{|}^{2m}||\mathrm{d}t|\ll
T(\log T)^{O(1)}.\end{split}$

We now deduce the estimation given in (3.6) from (3.9)-(3.11), upon taking
$D\gg C$ large enough. This completes the proof of the lemma.

###  3.6. Proof of Lemma 3.3

We deduce from (3.5) by symmetry and Hölder’s inequality that,

$\displaystyle\begin{split}&\Big{|}\int\limits_{\begin{subarray}{c}(1/2)\\\
|s|\leq(\log
T)^{D}\end{subarray}}\Big{|}\zeta(1/2+i(s+t))|\frac{1}{1+|s|}ds\Big{|}^{2m}\\\
\ll&\Big{|}\int_{0}^{(\log
T)^{D}}\frac{|\zeta(1/2+i(s+t))|}{t+1}dt\Big{|}^{2m}\\\
\leq&\bigg{(}\sum_{n\leq D\log\log
T+1}n^{-2m/(2m-1)}\bigg{)}^{2m-1}\sum_{n\leq D\log\log
T+1}\bigg{(}n\int_{e^{n-1}-1}^{e^{n}-1}\frac{|\zeta(1/2+i(s+t))|}{s+1}ds\bigg{)}^{2m}\\\
\ll&\sum_{n\leq D\log\log
T+1}\frac{n^{2m}}{e^{2nm}}\bigg{(}\int_{e^{n-1}-1}^{e^{n}-1}|\zeta(1/2+i(s+t))|ds\bigg{)}^{2m}.\end{split}$

We apply Proposition 2.4 to see that for any integer $k\geq 1$ and any real
numbers $2m\geq k+1,\varepsilon>0$,

$\displaystyle\sum_{n\leq D\log\log
T+1}\frac{n^{2m}}{e^{2nm}}\int^{2T}_{T}\bigg{(}\int_{e^{n-1}-1}^{e^{n}-1}|\zeta(1/2+i(s+t))|ds\bigg{)}^{2m}dt$
$\displaystyle\ll$ $\displaystyle T\sum_{n\leq D\log\log
T+1}\frac{n^{2m}}{e^{2nm}}\Big{(}(\log T)^{(m-1)^{2}}e^{kn}(\log
2n)^{O(1)}+(\log T)^{m^{2}-3m+3}(\log 2n)^{O(1)}(\log\log
T)^{O(1)}e^{2mn}\Big{)}$ $\displaystyle\ll$ $\displaystyle T(\log
T)^{(m-1)^{2}}.$

We now deduce from the above that (3.7) holds. This completes the proof of the
lemma.

###  3.7. Proof of Lemma 3.4

We apply the Cauchy-Schwarz inequality to see that

(3.12) $\displaystyle\begin{split}&\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{2m}dt\\\
\leq&\bigg{(}\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{2}dt\bigg{)}^{1/2}\bigg{(}\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{4m-2}dt\bigg{)}^{1/2}.\end{split}$

We first note that it follows from [iwakow, Theorem 9.1] that for arbitrary
complex numbers $a_{n}$, we have for $T,Z\geq 2$ and any $\varepsilon>0$,

$\displaystyle\begin{split}&\int^{2T}_{T}\bigg{|}\sum_{n\leq
Z}a_{n}n^{-it}\bigg{|}^{2}dt\ll(T+Z)\sum_{\begin{subarray}{c}n\leq
Z\end{subarray}}|a_{n}|^{2}.\end{split}$

We apply the above to $Z=Y$ and keep in mind our assumption that $Y\leq T$ to
see that

(3.13) $\displaystyle\begin{split}&\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{2}dt\ll(T+Y)\sum_{\begin{subarray}{c}n\leq
Y\end{subarray}}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}^{2}\ll
T\sum_{\begin{subarray}{c}Y(1-1/U)\leq n\leq
Y\end{subarray}}1+T\sum_{\begin{subarray}{c}0\leq n\leq
Y/U\end{subarray}}1\ll\frac{TY}{U}.\end{split}$

We next note that

(3.14) $\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\big{(}1-\Phi_{U}(\frac{n}{Y})\big{)}\bigg{|}^{4m-2}dt\ll\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\bigg{|}^{4m-2}dt+\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\Phi_{U}(\frac{n}{Y})\bigg{|}^{4m-2}dt.$

By the remark made in the paragraph below Theorem 1.1, we see that

(3.15) $\int^{2T}_{T}\bigg{|}\sum_{n\leq
Y}n^{-it}\Phi_{U}(\frac{n}{Y})\bigg{|}^{4m-2}dt\ll TY^{2m-1}(\log T)^{O(1)}.$

To estimate the first expression on the right-hand side of (3.14), we apply
Perron’s formula as given in [MVa1, Corollary 5.3] to see that

(3.16) $\displaystyle\begin{split}\sum_{n\leq Y}n^{-it}=&\frac{1}{2\pi
i}\int_{1+1/\log Y-iY}^{1+1/\log
Y+iY}\zeta(s+it)\frac{Y^{s}}{s}ds+R_{1}+R_{2},\\\ =&\frac{1}{2\pi
i}\int_{1+1/\log Y-iY}^{1/2-iY}+\frac{1}{2\pi
i}\int_{1/2-iY}^{1/2+iY}+\frac{1}{2\pi i}\int_{1/2+iY}^{1+1/\log
Y+iY}\zeta(s+it)\frac{Y^{s}}{s}ds+R_{1}+R_{2},\end{split}$

where

(3.17)
$\displaystyle\begin{split}R_{1}=&O\Big{(}\sum_{\begin{subarray}{c}Y/2<n<2Y\\\
n\neq Y\end{subarray}}\min(1,\frac{1}{|n-Y|})\Big{)}=O(\log Y),\\\
R_{2}=&O\Big{(}\frac{4^{1+1/\log Y}+Y^{1+1/\log Y}}{Y}\zeta(1+1/\log
Y)\Big{)}=O(\log Y).\end{split}$

Here the last estimation above follows from [MVa1, Corollary 1.17]. We now
consider the moments of the horizontal integrals in (3.16). We may assume that
$Y\geq 10$, otherwise the lemma is trivial. By symmetry we only need to
consider only one of them. Note that we have $|Y^{s}/s|\ll 1$ in that range
and $m\geq 1$, which allows us to apply Hölder’s inequality to get

(3.18) $\displaystyle\begin{split}\int^{2T}_{T}\bigg{|}\int_{1/2+iY}^{1+1/\log
Y+iY}\zeta(s+it)\frac{Y^{s}}{s}ds\bigg{|}^{4m-2}dt\ll&\int^{2T}_{T}\bigg{(}\int_{1/2+iY}^{1+1/\log
Y+iY}|\zeta(s+it)||\mathrm{d}s|\bigg{)}^{4m-2}dt\\\
\ll&\int_{1/2+iY}^{1+1/\log
Y+iY}\int^{2T}_{T}|\zeta(s+it)|^{4m-2}dt|\mathrm{d}s|\\\ \ll&T(\log
T)^{O(1)},\end{split}$

where the last estimation above follows from Lemma 2.3, which implies that for
$1/2\leq\Re(s)\leq 1+1/\log Y$, we have under RH,

$\int^{2T}_{T}|\zeta(s+it)|^{4m-2}dt\ll T(\log T)^{O(1)}.$

We treat the moments of the vertical integral in (3.16) using Hölder’s
inequality (by noting that $4m-2>4$), Proposition 2.4 and the assumption
$Y\leq(1-\varepsilon)T$ to see that

(3.19)
$\displaystyle\begin{split}&\int^{2T}_{T}\bigg{|}\int_{1/2-iY}^{1/2+iY}\zeta(s+it)\frac{Y^{s}}{s}ds\bigg{|}^{4m-2}dt\\\
\ll&Y^{2m-1}\int^{2T}_{T}\bigg{(}\int_{0}^{Y}\frac{|\zeta(1/2+i(s+t))|}{s+1}ds\bigg{)}^{4m-2}dt\\\
\ll&Y^{2m-1}\sum_{n\leq\log
Y+2}\frac{n^{4m-2}}{e^{(4m-2)n}}\int^{2T}_{T}\bigg{(}\int_{e^{n-1}-1}^{e^{n}-1}|\zeta(1/2+i(s+t))|ds\bigg{)}^{4m-2}dt\\\
\ll&Y^{2m-1}T(\log T)^{O(1)}\Big{(}\sum_{n\leq\log
Y+2}\frac{n^{4m-2}}{e^{(4m-2)n}}e^{n}+\sum_{n\leq\log Y+2}n^{4m-2}\Big{)}\\\
\ll&Y^{2m-1}T(\log T)^{O(1)}.\end{split}$

We conclude from (3.16)-(3.19) that

(3.20) $\int^{2T}_{T}\bigg{|}\sum_{n\leq Y}n^{-it}\bigg{|}^{4m-2}dt\ll
Y^{2m-1}T(\log T)^{O(1)}.$

We then deduce from (3.12)-(3.15), (3.20) and recall that we have $U=(\log
T)^{C},Y\leq(1-\varepsilon)T$ to see that the estimation given in (3.8) is
valid. This completes the proof of the lemma.

Acknowledgments. The author is supported in part by NSFC grant 11871082. This
work grows out of discussions with Changhao Chen and Nankun Hong on large
inequalities for Dirichlet polynomials when the author visited Anhui
University in April 2024. The author is debt to them for the inspiration of
this paper and many helpful suggestions on the writing of the manuscript.

## References


# WIP: Federated Learning for Routing in Swarm Based Distributed Multi-Hop
Networks

Martha Cash, Joseph Murphy, Alexander Wyglinski  Department of Electrical and
Computer Engineering, Worcester Polytechnic Institute
Worcester, MA
{mcash, jrmurphy<EMAIL_ADDRESS>
###### Abstract

Unmanned Aerial Vehicles (UAVs) are a rapidly emerging technology offering
fast and cost-effective solutions for many areas, including public safety,
surveillance, and wireless networks. However, due to the highly dynamic
network topology of UAVs, traditional mesh networking protocols, such as the
Better Approach to Mobile Ad-hoc Networking (B.A.T.M.A.N.), are unsuitable. To
this end, we investigate modifying the B.A.T.M.A.N. routing protocol with a
machine learning (ML) model and propose implementing this solution using
federated learning (FL). This work aims to aid the routing protocol to learn
to predict future network topologies and preemptively make routing decisions
to minimize network congestion. We also present an FL testbed built on a
network emulator for future testing of the proposed ML aided B.A.T.M.A.N.
routing protocol.

###### Index Terms:

Federated Learning, Routing, UAV Networks, B.A.T.M.A.N., Machine Learning

##  I Introduction

Unmanned Aerial Vehicles (UAVs) are a rapidly developing technology that has
been used in numerous applications, including transportation, traffic control,
surveillance, search and rescue, and disaster management [1]. Although UAV
technology has many advantages, numerous challenges still need to be addressed
to implement networking protocols for UAV-based infrastructures [2]. For one,
UAV networks are highly dynamic. They do not have a consistent topology making
communication, control, and path planning protocols designed for less dynamic
mobile ad-hoc networks (MANETs) less effective [1]. These challenges motivate
a need for routing protocols catered to UAVs. Ideally, these protocols should
be simple, have low overhead, and not require extensive global topology
knowledge [3]. Further, because of the dynamic nature of UAVs, these protocols
should make decisions based on the expected network topology rather than just
the current state of the network. These requirements make AI-based routing
protocols for UAVs appealing.

AI-based routing protocols are not a new area of research. For example, in [4]
a supervised feed-forward neural network (FFNN) was proposed to learn network
traffic history to adaptively route packets and improve heterogeneous network
control. Using the Open Shortest Path First (OSPF) routing algorithm, the
model input was an array representing the number of packets that were
forwarded through each node in the network, while the model’s output was the
interface to forward the packet along. Simulation results demonstrated the
effectiveness of the proposed FFNN approach and outperformed the OSPF
baseline. Similarly, in [5], Boltzmann machines (RBM) were proposed where the
input was characterized as the traffic pattern observed at each router. Like
[4], this approach outperforms the baseline OSPF routing algorithm. Finally,
the research presented in [6] proposed a neural network (NN) trained at each
link in the network, and the output of the NN was the likelihood of successful
packet delivery if the packet was forwarded along that link. The authors
propose using buffer capacity, number of successful packet transfers, and node
popularity.

All of these approaches focus on applying common machine learning (ML)
techniques to the AI-based routing protocol problem in homogeneous networks,
which will not work UAV networks. However, an emerging ML technique, federated
learning (FL) [7], has yet to be explored as a solution for heterogeneous
networks. In this work, we propose an FL-based approach to the AI-based
routing protocol problem, specifically for UAV swarms. We narrow our focus to
the Better Approach to Mobile Ad-hoc Networking (B.A.T.M.A.N.) protocol [8],
propose modifying the algorithm using a NN model, and characterize the dataset
necessary for this problem. Finally, we present an FL emulation environment
built on the Extendable Mobile Ad-hoc Network Emulator (EMANE) [9] that will
be used for testing the proposed solution.

The remainder of this paper is organized as follows: Section II provides an
overview of the system model. Section III presents the proposed solution and
describes FL in more detail. Section IV shows the preliminary results of FL
setup in a network emulator and NN model, and Section V discusses conclusions
and future research directions.

##  II System Model

Figure 1: A two node UAV network. At each time step, the link cost between
$D_{B}$ and $D_{F}$ increases while the link cost between $D_{B}$ and $D_{F}$
is constant. The B.A.T.M.A.N. routing protocol will continue to select the
route with the lower link cost, branch $(b)$. The proposed solution, branch
$(a)$, would preemptively switch routes to manage network congestion and
select an alternate route despite being unfavorable at the current time step.

The initial system, shown in Fig. 1, is a simple two-node UAV network, which
will be extended to a multi-node heterogeneous UAV network in future work. The
network employs the B.A.T.M.A.N. routing protocol for communication between
UAVs. The aim of this study, illustrated branch $(a)$ in Fig. 1, is to
demonstrate the feasibility of an ML-aided B.A.T.M.A.N. protocol to improve
network congestion by predicting when to switch routes, even if switching to a
route with a higher cost is not immediately beneficial.

###  II-A B.A.T.M.A.N. Protocol

As a baseline, the network uses the B.A.T.M.A.N routing protocol [10].
B.A.T.M.A.N. was designed to address the challenges of routing in mobile ad-
hoc networks (MANETs), such as frequent topology changes and the lack of a
central authority to coordinate routing. Rather than maintain information
about the global network topology, B.A.T.M.A.N. only requires nodes to
maintain information about the best next hop to its immediate neighbors. The
network is flooded with originator messages (OGMs). OGMs routed through good
paths are received by nodes quicker than those transmitted on poor quality
routes, informing the nodes in the network which immediate neighbor has the
best route to transmit across. The routing tables are configured by selecting
the best next hop to the originator node [11].

However, there are a few drawbacks to the B.A.T.M.A.N. routing algorithm. For
one, if the network contains a substantial number of nodes, B.A.T.M.A.N. can
generate a large amount of overhead, as each node must re-broadcast the OGM to
its neighbors. This can lead to increased network congestion and reduced
overall efficiency. Additionally, in scenarios where the source or destination
of the packet is in motion, B.A.T.M.A.N. can suffer from higher delays, which
is undesirable if the network topology is highly dynamic. Finally, the
drawback we focus on in this work is that B.A.T.M.A.N. is a threshold-based
routing protocol, demonstrated in Fig. 1. As a result, the node will always
choose the next hop with the best route, even if conditions on the current
best route are degrading. Waiting to change routes until the threshold is met
can cause bottlenecks in the network [12].

##  III Proposed Solution

We propose two solutions. First, we integrate federated learning (FL) to a
wireless drone swarm network and propose an ML model and data set for
enhancing the B.A.T.M.A.N. routing protocol. Next, we demonstrate an FL
simulation environment built on the EMANE emulation environment, which will be
integrated with the proposed ML model in future work to investigate a larger
UAV network, and to introduce movement among the UAVs, since this will be a
feature compensated for by the ML model.

###  III-A Machine Learning Model & Dataset

We consider a supervised learning approach for this work. The objective of
supervised learning is to learn a mapping, or function approximation,
$\hat{\mathbb{F}}(\mathbf{x},\mathbf{y})$, between a set of samples, $x_{i}\in
X$, and their labels, $y_{i}\in Y$, where $X$ and $Y$ are the sample space and
label space, respectively. Ideally, $\hat{\mathbb{F}}(x,y)$ takes a set of new
samples, $\mathbf{x^{*}}$ and produces the correct label, $\mathbf{y^{*}}$.
The quality of the mapping is determined by the loss function,
$L(y^{*}_{i},\hat{y}_{i})$, where $y^{*}_{i}$ is the true label of the new
sample, and $\hat{y}_{i}$ is the output of
$\hat{\mathbb{F}}({x^{*}_{i}},\cdot)$ [13]. An accurate function approximation
is quantified by a low loss value.

Since the B.A.T.M.A.N. routing protocol does not maintain a history of route
conditions (i.e. link cost, throughput), we need to modify the B.A.T.M.A.N.
algorithm to include a memory element. The model should learn a history of the
prior link costs for each route, and the route the node selected. These
requirements make the long short-term memory (LSTM) model, a type of recurrent
neural network (RNN) that is designed to learn long-term dependencies in
sequential data, appropriate for this task [14].

The input to the LSTM model is a two dimensional array of the history of the
link cost at each neighboring route from $D_{B}$, (see Fig. 1). However, this
approach can be extended to $n$ dimensions for $n$ many neighbors in a more
complex network. The corresponding labels are a history of the selected route
for transmission. Instead of feeding the entire history of the network to the
model, we implement a windowing technique. For example, if the window size is
set to 4, then four prior time steps are fed into the model for training. We
can treat this as a classification problem and use the binary cross entropy
loss function since we are training the LSTM to select which of the two routes
from $D_{B}$ to transmit across.

###  III-B Federated Learning Approach

Traditional machine learning approaches are centralized, meaning a model is
trained on a central dataset that is collected and stored in a central
location. However, this approach assumes all devices on the network have the
same computational capabilities and network resources, which is often not the
case for UAV networks [2]. As a result, a FL-based distributed ML technique is
employed in this research where training of a global model is performed on
data distributed across many UAVs in various locations. For generality, we
assume there are $J$ UAVs in the network. Each UAV for $j\in J$ observes a
unique dataset $\mathbf{x_{j}}=\\{x_{j1},x_{j2},\ldots,x_{jN}\\}$. Since we
use a supervised learning approach, we assume a single input sample, $X_{jn}$
corresponds to a single output $y_{jn}\in\\{y_{j1},y_{j2},\ldots,y_{jN}\\}$.
The sets $\mathbf{x}_{j}$ and $\mathbf{y}_{j}$ are used to train the local ML
model at each UAV. Let $\mathbf{w}_{j}\in\mathbf{w}$ denote the corresponding
model parameters at the $j^{th}$ UAV. Then, the FL objective function can be
employed, which is defined as [15]:

$\operatorname*{arg\,min}_{\mathbf{w}\in\mathbb{R}^{d}}F(\mathbf{w})=\frac{1}{N}\sum_{j}^{J}\sum_{n=1}^{N_{j}}f(\mathbf{w}_{j},x_{jn},y_{jn}).$
(1)

We solve (1) via the following steps: First, all UAVs are initialized with
random parameters. Each UAV trains on its respective training sets, $x_{jn}$
and $y_{jn}$. After one epoch of training, the FL parameters at the $j^{th}$
UAV are sent to a central server. Once all parameters are received, the
central server aggregates the parameters according to the following
expression:

$w_{global}=\frac{1}{J}\sum_{j=1}^{J}\mathbf{w_{j}}$ (2)

The global parameters, $w_{global}$, are sent back to the UAVs, and the
training process repeats for $K$ epochs, or until $F(\mathbf{w})$ has
converged to the optimal parameters, $w^{*}$. Fig. 2 summarizes this approach.

###  III-C Federated Learning in EMANE

This works also aims to build an FL emulation environment, with the future
goal of integrating the emulator with the proposed ML model to test the
feasibility of the proposed solution. A system diagram of the emulator is
shown in Fig. 2. At the core of the emulator is EMANE, which allows for the
creation of Network Emulation Modules (NEMs) to model different radio
interface types. In turn, these can be incorporated into a real-time emulation
running in a distributed environment and allow the direct integration of
standard software, such as PyTorch, for handling ML tasks.

For the results in this work, we construct a simulator with three NEMs,
similar to the model setup in Fig. 1. One NEM is designated as the central
server. The remaining NEMs carry out the FL task. However, this could be
generalized to $M$ nodes, see Fig. 2.

Figure 2: Federated Learning Setup in EMANE

##  IV Preliminary Results

###  IV-A CIFAR10 Baseline

We use the CIFAR10 dataset [16] to test the performance of the network
emulator proposed in Section III-C and compare the results to a centralized ML
approach. To accurately model a distributed UAV environment, we randomly
sample the CIFAR10 dataset so that each NEM in the emulator is not training on
an identical dataset. We train the two edge nodes for 10 epochs. The results
are presented in Fig. 4, demonstrating that the FL emulator (’x’ curve)
approximated the traditional ML approach (’o’ curve). Because the nodes are
not training on an identical dataset, the FL training loss curve should
approximate the behavior of the centralized ML model training curve, but
should not identically match the curve. The results presented in Fig. 4
indicate the system setup proposed in Section III-C works and is suitable for
future integration with the ML model presented in Section III-A in addition to
expanding the number of NEMs included in the emulation.

###  IV-B Machine Learning Aided B.A.T.M.A.N.

We attempt to hand generate a simple testing and training set consisting of 50
time samples of link costs across the two routes branching from $D_{B}$ as
demonstrated in Fig. 3. We construct the LSTM to have two recurrent layers.
The input sequence length is set to four. We use a batch size of five and
train the model for ten epochs using the binary cross entropy loss function
and the ADAM optimizer with a learning rate of 0.01. Initial testing accuracy
results show the LSTM has an 100% classification accuracy. This is due to the
memoryless nature of the hand generated data set. Therefore, in this current
state, it is not possible to accurately estimate the practicality of the
proposed model.

Figure 3: Demonstration of Hand Generated Training Set Figure 4: Federated
Learning in EMANE Results: FL Model Approximates Centralized Model Performance

##  V Conclusion

This paper proposed a machine learning (ML) solution implemented via federated
learning (FL) to modify the B.A.T.M.A.N. routing protocol for unmanned aerial
vehicle networks. We presented a FL testbed built on the network emulator
EMANE and used the CIFAR10 dataset to compare the FL testbed to a traditional
centralized ML approach. The baseline results proved the testbed works as a
proof of concept for future work. We also propose modifying the B.A.T.M.A.N.
algorithm using a long short term memory model. However, current results do
not accurately reflect the viability of this modification. Therefore, future
work will need to generate a dataset from a simulation of the UAV swarm using
EMANE to assess this approach better.

## Acknowledgment

This research was sponsored by the DEVCOM Analysis Center and was accomplished
under Cooperative Agreement Number W911NF-22-2-0001. The views and conclusions
contained in this document are those of the authors and should not be
interpreted as representing the official policies, either expressed or
implied, of the Army Research Office or the U.S. Government. The U.S.
Government is authorized to reproduce and distribute reprints for Government
purposes, notwithstanding any copyright notation herein.

## References

  * [1] A. Rovira-Sugranes, A. Razi, F. Afghah, and J. Chakareski, “A review of ai-enabled routing protocols for uav networks: Trends, challenges, and future outlook,” _Ad Hoc Networks_ , vol. 130, p. 102790, 2022.
  * [2] L. Gupta, R. Jain, and G. Vaszkun, “Survey of important issues in uav communication networks,” _Commun. Surveys Tuts._ , vol. 18, no. 2, p. 1123–1152, apr 2016. [Online]. Available: https://doi.org/10.1109/COMST.2015.2495297
  * [3] M. Khaledi, A. Rovira-Sugranes, F. Afghah, and A. Razi, “On greedy routing in dynamic uav networks,” in _2018 IEEE International Conference on Sensing, Communication and Networking (SECON Workshops)_. IEEE, 2018, pp. 1–5.
  * [4] N. Kato, Z. M. Fadlullah, B. Mao, F. Tang, O. Akashi, T. Inoue, and K. Mizutani, “The deep learning vision for heterogeneous network traffic control: Proposal, challenges, and future perspective,” _IEEE Wireless Communications_ , vol. 24, no. 3, pp. 146–153, 2017.
  * [5] B. Mao, Z. M. Fadlullah, F. Tang, N. Kato, O. Akashi, T. Inoue, and K. Mizutani, “Routing or computing? the paradigm shift towards intelligent computer network packet transmission based on deep learning,” _IEEE Transactions on Computers_ , vol. 66, no. 11, pp. 1946–1960, 2017.
  * [6] D. K. Sharma, S. K. Dhurandher, I. Woungang, R. K. Srivastava, A. Mohananey, and J. J. P. C. Rodrigues, “A machine learning-based protocol for efficient routing in opportunistic networks,” _IEEE Systems Journal_ , vol. 12, no. 3, pp. 2207–2213, 2018.
  * [7] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-efficient learning of deep networks from decentralized data,” in _Artificial intelligence and statistics_. PMLR, 2017, pp. 1273–1282.
  * [8] A. Klein, L. Braun, and F. Oehlmann, “Performance study of the better approach to mobile adhoc networking (b.a.t.m.a.n.) protocol in the context of asymmetric links,” in _2012 IEEE International Symposium on a World of Wireless, Mobile and Multimedia Networks (WoWMoM)_ , 2012, pp. 1–7.
  * [9] U. N. R. Laboratory, “Extendable mobile ad-hoc network emulator (emane).” [Online]. Available: https://www.nrl.navy.mil/Our-Work/Areas-of-Research/Information-Technology/NCS/EMANE/
  * [10] M. Lindner, S. Eckelmann, S. Wunderlich, M. Hundebøll, A. Quartulli, and L. Lüssing, “The b.a.t.m.a.n. project.” [Online]. Available: http://www.open-mesh.org/
  * [11] B. Sliwa, S. Falten, and C. Wietfeld, “Performance evaluation and optimization of B.A.T.M.A.N. V routing for aerial and ground-based mobile ad-hoc networks,” _CoRR_ , vol. abs/1901.02298, 2019. [Online]. Available: http://arxiv.org/abs/1901.02298
  * [12] D. Seither, A. König, and M. Hollick, “Routing performance of wireless mesh networks: A practical evaluation of batman advanced,” in _2011 IEEE 36th Conference on Local Computer Networks_ , 2011, pp. 897–904.
  * [13] A. Valadarsky, M. Schapira, D. Shahaf, and A. Tamar, “A machine learning approach to routing,” _arXiv preprint arXiv:1708.03074_ , 2017.
  * [14] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [15] T. Zeng, O. Semiari, M. Mozaffari, M. Chen, W. Saad, and M. Bennis, “Federated learning in the sky: Joint power allocation and scheduling with uav swarms,” in _ICC 2020 - 2020 IEEE International Conference on Communications (ICC)_ , 2020, pp. 1–6.
  * [16] A. Krizhevsky, G. Hinton _et al._ , “Learning multiple layers of features from tiny images,” 2009.


# Universal ratios of critical amplitudes in the Potts model universality
class

B. Berche  P. Butera  W. Janke  L. Shchur  Laboratoire de Physique des
Matériaux, Université Henri Poincaré Nancy 1, BP 239, F-54506 Vandœuvre les
Nancy Cedex, France  Istituto Nazionale di Fisica Nucleare, Sezione di Milano-
Bicocca, Piazza delle Scienze 3, 20126, Milano, Italia  Institut für
Theoretische Physik, Universität Leipzig, Postfach 100 920, 04009 Leipzig,
Germany  Landau Institute for Theoretical Physics, 142432 Chernogolovka,
Russia<EMAIL_ADDRESS>
###### Abstract

Monte Carlo (MC) simulations and series expansions (SE) data for the energy,
specific heat, magnetization, and susceptibility of the three-state and four-
state Potts model and Baxter-Wu model on the square lattice are analyzed in
the vicinity of the critical point in order to estimate universal combinations
of critical amplitudes. We also form effective ratios of the observables close
to the critical point and analyze how they approach the universal critical-
amplitude ratios. In particular, using the duality relation, we show
analytically that for the Potts model with a number of states $q\leq 4$, the
effective ratio of the energy critical amplitudes always approaches unity
linearly with respect to the reduced temperature. This fact leads to the
prediction of relations among the amplitudes of correction-to-scaling terms of
the specific heat in the low- and high-temperature phases. It is a common
belief that the four-state Potts and the Baxter-Wu model belong to the same
universality class. At the same time, the critical behavior of the four-state
Potts model is modified by logarithmic corrections while that of the Baxter-Wu
model is not. Numerical analysis shows that critical amplitude ratios are very
close for both models and, therefore, gives support to the hypothesis that the
critical behavior of both systems is described by the same renormalization
group fixed point.

###### keywords:

Potts model; Baxter-Wu model; Critical exponents; Critical amplitudes;
Universality; Monte Carlo simulations; Series Expansions; Renormalization
Group

###### PACS:

0.50.+q, 75.10.-b

††journal: Computer Physics Communications

The fixed points of the renormalization group define the universal behavior of
a system through a set of critical exponents and universal combinations of
critical amplitudes [1]. The universality concept divides all systems at
criticality into a number of universality classes. It is instructive to know
the set of values of the critical exponents and of the universal combinations
of critical amplitudes for a given universality class.

The two-dimensional Potts model [2] is the simplest model which exhibits a
phase transition. It is solved exactly at the critical point for any number of
spin components $q$ and it is known that for $q\leq 4$ it has a continuous
phase transition while for $q>4$ the phase transition is of the first order.
The model has a great theoretical interest as new theories may be tested in
this model.

At the same time, these models may have some practical interest as they may be
realized in an adsorbate lattice placed onto a clean crystalline surface. The
full classification of such systems with continuous transitions is known
theoretically [3]. There are experiments in which some of them realize the
3-state and 4-state Potts models [4] and the critical exponents can be
experimentally estimated.

Critical exponents for the Potts model with $q\leq 4$ can be computed exactly
by different theoretical techniques [5, 6]. The values of the thermal critical
exponents and of the magnetic critical exponents follow from the
identification of the dimensions of the conformal algebra operators [6].

Nowadays, there is no doubt on the values of the leading critical exponents
whereas the values of the correction-to-scaling exponents are still under
discussion, as well as the values of the universal ratios of the critical
amplitudes. Our presentation is intended to give a short review of the
research on the subject.

The Potts model Hamiltonian [2] (see review [7] for details) can written as
$H=-\sum_{\langle ij\rangle}\delta_{s_{i}s_{j}}\;,$ where $s_{i}$ is a spin
variable taking integer values between $0$ and $q{-}1$, and the sum is
restricted to the nearest neighbor sites $\langle ij\rangle$ on the square
lattice.

Close to the critical temperature $T_{c}$ at which the continuous phase
transition occurs, the residual magnetization $M$ and the singular part of the
reduced susceptibility $\chi$ and of the specific heat $C$ of the system in
zero external field are characterized by the critical exponents $\beta$,
$\gamma$, and $\alpha$ and by the critical amplitudes $B$, $\Gamma_{\pm}$, and
$A_{\pm}$

$\displaystyle M(\tau)$ $\displaystyle\approx$ $\displaystyle
B(-\tau)^{\beta},\ \tau<0$ (1) $\displaystyle\chi_{\pm}(\tau)$
$\displaystyle\approx$ $\displaystyle\Gamma_{\pm}|\tau|^{-\gamma},$ (2)
$\displaystyle C_{\pm}(\tau)$ $\displaystyle\approx$
$\displaystyle\frac{A_{\pm}}{\alpha}|\tau|^{-\alpha},$ (3)

where $\tau$ is the reduced temperature $\tau=(T-T_{c})/T$ and the labels
$\pm$ refer to the high-temperature and low-temperature sides of the critical
temperature $T_{c}$. The critical amplitudes are not universal by themselves
but some combinations of them, f.e., $A_{+}/A_{-}$, $\Gamma_{+}/\Gamma_{-}$,
and $\Gamma_{+}A_{+}/B^{2}$, are universal [1] due to the scaling laws.

On the square lattice, in zero field, the model is self-dual. The duality
relation

$\left(e^{\beta}-1\right)\left(e^{\beta^{*}}-1\right)=q$ (4)

fixes the inverse critical temperature to $\beta_{c}{=}\ln(1{+}\sqrt{q})$. The
values $E(\beta)$ and $E(\beta^{*})$ of the internal energy at dual
temperatures are simply related through

$\left(1-e^{-\beta}\right)E(\beta)+\left(1-e^{-\beta^{*}}\right)E(\beta^{*})=2.$
(5)

Dual reduced temperatures $\tau$ and $\tau^{*}$ can be defined by
$\beta{=}\beta_{c}(1{-}\tau)$ and $\beta^{*}{=}\beta_{c}(1{+}\tau^{*})$. Close
to the critical point, $\tau$ and $\tau^{*}$ coincide through linear order,
since
$\tau^{*}{=}\tau{+}\frac{\ln(1{+}\sqrt{q})}{\sqrt{q}}\tau^{2}{+}O(\tau^{3})$.

The ratio of the free energy critical amplitudes $A_{+}/A_{-}$ is equal to
unity due to duality. Moreover, duality relations may be used to understand
the dependence on temperature of the effective amplitude functions which may
be constructed from the energy in the symmetric phase $E_{+}(\tau)$ and in the
ordered phase $E_{-}(\tau^{*})$

$\displaystyle A_{+}(\tau)$ $\displaystyle=$
$\displaystyle\alpha(1-\alpha)\beta_{c}(E_{+}(\tau)-E_{0})\tau^{\alpha-1},$
(6) $\displaystyle A_{-}(\tau^{*})$ $\displaystyle=$
$\displaystyle\alpha(1-\alpha)\beta_{c}(E_{0}-E_{-}(\tau^{*}))(\tau^{*})^{\alpha-1}$
(7)

as an effective amplitude ratio

$\frac{A_{+}(\tau)}{A_{-}(\tau^{*})}=\frac{(E_{+}(\tau)-E_{0})\tau^{\alpha-1}}{(E_{0}-E_{-}(\tau^{*}))(\tau^{*})^{\alpha-1}},$
(8)

where the constant $E_{0}$ is the value of the energy at the transition
temperature, $E_{0}{=}E(\beta_{c}){=}{-}1{-}1/\sqrt{q}$.

Evaluating expression (8) for small $\tau$ and denoting
$\alpha_{q}=-E_{0}\beta_{c}e^{-\beta_{c}}=\frac{\ln(1+\sqrt{q})}{\sqrt{q}}$,
we obtain

$\displaystyle\frac{A_{+}(\tau)}{A_{-}(\tau^{*})}$ $\displaystyle=$
$\displaystyle 1+(3-\alpha)\alpha_{q}\tau+O(\tau^{1+\alpha}).$ (9)

Note the linear dependence on $\tau$ of the effective amplitude ratio.

The 2-state Potts model is equivalent to the Ising model which was solved
exactly [8] (see Ref. [9] for details). The susceptibility behavior was
understood in the paper by Wu, McCoy, Tracy and Barouch [10]. It turns out
that there exist only integer corrections to scaling (for a recent and
detailed discussions we refer readers to Refs. [11, 12, 13]). Values of the
critical exponents and of some amplitude ratios are presented in Table 1.

Table 1: Exact values of critical exponents and ratios of critical amplitudes for the Ising model (2-state Potts model). $\nu$ | $\alpha$ | $\beta$ | $\gamma$ | $A_{+}/A_{-}$ | $\Gamma_{+}/\Gamma_{L}$ | $R_{C}^{+}=\Gamma_{+}A_{+}/B^{2}$
---|---|---|---|---|---|---
1 | 0 | $1/8$ | $7/4$ | 1 | 37.69365… | 0.318569…

The critical behavior of the susceptibility reads as

$\chi(\tau)=\Gamma_{\pm}|\tau|^{-\gamma}{\cal X}_{corr}(|\tau|^{\Delta})+{\cal
Y}_{bt}(\tau),$ (10)

where ${\cal X}_{corr}(|\tau|^{\Delta})$ is the correction-to-scaling function
and ${\cal Y}_{bt}(\tau)$ represents an analytic expression (“background
term”) which accounts for non-singular contributions to susceptibility.

Set of values of the thermal critical exponents $x_{\epsilon_{n}}$ and of the
magnetic critical exponents $x_{\sigma_{n}}$ are known analytically [5, 6]

$x_{\epsilon_{n}}=\frac{n^{2}+ny}{2-y},\;\;\;x_{\sigma_{n}}=\frac{\left(2n-1\right)^{2}-y^{2}}{4(2-y)}$
(11)

in terms of the variable $y$ linked to the number of states $q$ by
$\cos\frac{\pi y}{2}=\frac{1}{2}\sqrt{q}$.

For the 3-state Potts model there is a finite number of correction terms [6],
$x_{\epsilon_{2}}=14/5$, $x_{\epsilon_{2}}=6$, and $x_{\sigma_{2}}=4/3$. The
leading correction-to-scaling contribution is
$\Delta=-(2-x_{\epsilon_{2}})/(2-x_{\epsilon_{1}})=2/3$ and it was first
supported in numerical simulation [14].

Clear evidence for these leading correction to the scaling behavior may be
seen in Figure 1, where we plot the difference of the longitudinal and
transverse susceptibilities $(\chi_{L}-\chi_{T})$ (note the cancelation of
background terms (10) in the difference) multiplied by the leading behavior
factor $|\tau|^{\gamma}$, as a function of $|\tau|^{2/3}$. The ratio of the
susceptibilities $(\chi_{T}/\chi_{L})$ is shown in Figure 2 and may be used to
estimate the universal ratio of associated amplitudes $\Gamma_{T}/\Gamma_{L}$.

Analytical predictions for the amplitude ratios in Potts model for $q=2,3$,
and 4 were given in the papers [15, 16]. The values are shown in the Table 2
together with numerical estimations from Monte Carlo (MC) and series
expansions (SE) analyses [17, 18, 19]. The coincidence of the data is a good
indication for the validity of both two-kink approximation [15, 16] to the
exact scattering theory [20] for 3-state Potts model and of the analysis of MC
data and series expansions.

Table 2: Exact values of critical exponents and estimates of the ratios of critical amplitudes for the 3-state Potts model. $\nu$ | $\alpha$ | $\beta$ | $\gamma$ | $\Gamma_{+}/\Gamma_{L}$ | $\Gamma_{T}/\Gamma_{L}$ | $R_{C}^{+}{=}\Gamma_{+}A_{+}/B^{2}$ | Remark
---|---|---|---|---|---|---|---
5/6 | 1/3 | $1/9$ | $13/9$ | - | - | - | Exact result
|  |  |  | 13.848 | 0.327 | 0.1041 | [15, 16]
|  |  |  | 13.83(9) | 0.325(2) |  |  [19] \- SE
|  |  |  | 13.83(9) | 0.3272(7) | 0.1044(8) |  [18] \- SE
|  |  |  | 13.86(12) | 0.322(3) | 0.1049(29) |  [18] \- MC

The analysis of the 4-state Potts model is much more complicated because in
addition to the corrections to scaling there are confluent logarithmic
corrections [21, 22]. The result of the analysis of MC data and SE [23, 24,
19] is shown in the Table 3 together with the analytical estimates [15, 16].

Conformal field theory predicts the set of renormalization group (RG)
exponents $y_{\epsilon_{n}}{=}2{-}x_{\epsilon_{n}}{=}$3/2, 0, -5/2, … . The
leading correction-to-scaling exponent $y_{\epsilon_{2}}$ vanishes and gives
rise to a logarithmic behavior [21]. In our recent publication [24, 25] we
revised the renormalization group equations and included in our analysis the
known form of the logarithmic corrections and of the next-to-leading
corrections, taking into account the width of the temperature region window
examined. The set of magnetic exponents for the Potts model
$x_{\sigma_{n}}$=1/8, 9/8, 25/8, … translates into the magnetization exponent
$\beta{=}1/12$ and the leading correction-to-scaling exponent
$\Delta_{\sigma}{=}2/3$. Finally, the following behavior of the susceptibility
is assumed

$\chi_{+}(\tau)=\Gamma_{+}\tau^{-7/6}{\cal
G}^{3/4}(-\ln\tau)(1+a_{+}\tau^{2/3}+b_{+}\tau)+D_{+},$ (12)

where the function $\cal G$ contains a universal correction function $\cal E$
[22, 24, 25] and the leading nonuniversal correction function $\cal F$

$\displaystyle{\cal G}(-\ln|\tau|)$ $\displaystyle=$
$\displaystyle(-\ln|\tau|)\times{\cal E}(-\ln|\tau|)\times{\cal
F}(-\ln|\tau|),$ (13) $\displaystyle{\cal E}(-\ln|\tau|)$ $\displaystyle=$
$\displaystyle\left(1+\frac{3}{4}\frac{\ln(-\ln|\tau|)}{-\ln|\tau|}\right)$
(14) $\displaystyle\
\times\left(1-\frac{3}{4}\frac{\ln(-\ln|\tau|)}{-\ln|\tau|}\right)^{-1}\left(1+\frac{3}{4}\frac{1}{(-\ln|\tau|)}\right),$
$\displaystyle{\cal F}(-\ln|\tau|)$ $\displaystyle\simeq$
$\displaystyle\left(1+\frac{C_{1}}{-\ln|\tau|}+\frac{C_{2}\ln(-\ln|\tau|)}{(-\ln|\tau|)^{2}}\right)^{-1}.$
(15)

We fit our data to estimate the amplitude $\Gamma_{+}$, the coefficient of the
leading correction to scaling $a_{+}$ in Eq. (12), and coefficients $C_{1}$
and $C_{2}$ in Eq. (15).

It is obvious that the logarithmic corrections (the whole function $\cal
G(-\ln|\tau|)$) cancels in simple ratios, like $A_{+}/A_{-}$,
$\Gamma_{+}/\Gamma_{L}$, $\Gamma_{T}/\Gamma_{L}$, etc. This has been
demonstrated analytically for the effective ratio $A_{+}/A_{-}$ (see Eq. 9).
We note also that the RG analysis predicts [24, 25] powers of logarithmic
corrections to specific heat $\alpha^{\prime}{=}{-}1$, susceptibility
$\gamma^{\prime}{=}3/4$, and magnetization $\beta^{\prime}{=}{-}1/8$ such that
they cancel in all universal ratios. For example, the universal amplitude
ratio $R_{C}^{-}$ may be calculated as the limit of the ratio of functions

$R_{C}^{-}=\lim_{\tau\rightarrow
0}\tau\frac{(E_{-}(|\tau|)-E_{0})\chi_{-}(|\tau|)}{M(|\tau|)^{2}}\alpha(\alpha-1)\beta_{c}$
(16)

where $E_{0}=E(0)=\sqrt{2}$. One can check that in the ratio (16) not only
powers of $|\tau|$ cancel but also powers of $\cal E$. In the ratio, the
magnetization $M$ and the energy difference $E_{-}(|\tau|)-E_{0}$ have only
singular contributions and the only systematic deviation may come from the
background correction to susceptibility $\chi_{-}(|\tau|)$. It was shown in
[18] that the contribution from this background correction is negligible in
the critical region window, and the estimator (16) tends to the value
$0.0055(1)$ as $\tau\rightarrow 0$.

Figure 1: 3-state Potts model. Difference of susceptibilities
$(\chi_{L}-\chi_{T})|\tau|^{\gamma}$ as function of $|\tau|^{2/3}$. The almost
linear dependence supports the value 2/3 for the power of the leading
correction to scaling.

The Baxter-Wu [26] model is defined on a triangular lattice, with spins
$\sigma_{i}$ located at vertices. The three spins forming a triangular face
are coupled with a strength $J$, and the Hamiltonian reads

${\cal H}=-J\sum_{\rm faces}\sigma_{i}\sigma_{j}\sigma_{k},$ (17)

where the summation extends over all triangular faces of the lattice, both
pointing up and down. The ground state is four-fold degenerate and the
critical exponents are found to be the same as for the 4-state Potts model.

The exact behavior of the magnetization, energy, and specific heat are known
analytically [26, 27, 28]. An analysis of Monte Carlo data was performed by
two of us [29] and preliminary estimates shows that the values of the
susceptibility amplitude-ratio $\Gamma_{+}/\Gamma_{L}\approx 6.9$ and of the
ratio $R_{C}^{-}\approx 0.005$ are very near to those obtained from our
analysis of MC and SE data for the 4-state Potts model (see Table 3). We have
to note that logarithmic corrections to scaling are absent in the critical
behavior of Baxter-Wu model and this gives us more confidence in our analysis.

Delfino and Grinza [30] use the same analytical approach as in [15] to study
the Ashkin-Teller model which also belongs to the 4-state Potts model
universality class with some particular choice of coupling constants. This
leads to the estimatation $\Gamma_{+}/\Gamma_{L}\approx 4.02$ and
$\Gamma_{T}/\Gamma_{L}\approx 0.129$. This result is also very near to those
for 4-state Potts model (see second entry in Table 3.)

A possible explanation of the deviation of our results from the analytical
predictions may be explained as follows: the two-kink approximation is exact
for the 2-state Potts model, it gives good accuracy for the 3-state Potts
model, but it may be insufficient to produce accurate values for 4-state Potts
model. Further analyses have to be done to resolve the contradiction among
these results.

Table 3: Exact values of critical exponents and estimations of the ratios of critical amplitudes for the 4-state Potts model. $\nu$ | $\alpha$ | $\beta$ | $\gamma$ | $\Gamma_{+}/\Gamma_{L}$ | $\Gamma_{T}/\Gamma_{L}$ | $R_{C}^{-}{=}\Gamma_{-}A_{-}/B^{2}$ | Remark
---|---|---|---|---|---|---|---
2/3 | 2/3 | $1/12$ | $7/6$ | - | - | - | Exact result
|  |  |  | 4.013 | 0.129 | 0.00508 | [15, 16]
|  |  |  | 3.5(4) | 0.11(4) |  |  [19] \- SE
|  |  |  | 3.14(70) |  | 0.0068(9) |  [23] \- MC
|  |  |  | 6.93(6) | 0.1674(30) | 0.00512(13) |  [24, 25] \- MC
|  |  |  | 6.30(1) | 0.1511(24) | 0.00531(5) |  [24, 25] \- SE

BB and WJ acknowledge partial support within the Graduate School “Statistical
Physics of Complex Systems” of DFH-UFA under Contract No. CDFA-02-07.
Financial support within a common research program between the Landau
Institute and the Ecole Normale Supérieure de Paris, Paris Sud University is
also gratefully acknowledged.

## References

  * [1] V. Privman, P.C. Hohenberg, A. Aharony, in Phase Transitions and Critical Phenomena, Vol. 14, edited by C. Domb and J.L. Lebowitz (Academic, New York, 1991).
  * [2] R.B. Potts, Proc. Camb. Phil. Soc. 48 (1952) 106.
  * [3] E. Domany, M. Schick, J. Walker, and R.B. Griffiths, Phys. Rev. B 18 (1978) 2209; E. Domany and M. Schick, Phys. Rev. B 20 (1979) 3828; C. Rottman, Phys. Rev. B 24 (1981) 1482.
  * [4] A. Aharony, K. A. Muller and W. Berlinger, Phys. Rev. Lett. 38 (1977) 33;H. Pfnür and P. Piercy, Phys. Rev. B 41 (1990) 582; M. Sokolowski and H. Pfnür, Phys. Rev. Lett. 49 (1994) 7716; Y. Nakajima, et al, Phys. Rev. B 55(1997) 8129; C. Voges and H. Pfnür, Phys. Rev. B 57 (1998) 3345.
  * [5] M.P.M. den Nijs, J. Phys. A 12 (1979)1857; R.B. Pearson, Phys. Rev. B 22 (1980) 2579; B. Nienhuis, J. Phys. A 15 (1982) 199; B. Nienhuis, J. Stat. Phys. 34 (1984) 731; B. Nienhuis, in Phase Transitions and Critical Phenomena,Vol. 11, edited by C. Domb and J.L. Lebowitz(Academic Press, London, 1987).
  * [6] Vl.S. Dotsenko, Nucl. Phys. B 235 (1984) 54; Vl.S. Dotsenko and V.A. Fateev, Nucl. Phys. B 240 (1984) 312.
  * [7] F.Y. Wu, Rev. Mod. Phys. 54 (1982) 235.
  * [8] L. Onsager, Phys. Rev. ]65 (1944) 117.
  * [9] B. M. McCoy and T.T. Wu, The two-dimensional Ising model (Harvard Uni. Press, Camdridge, Massachusetts, 1973)
  * [10] E. Barouch, B.M. McCoy, and T.T. Wu, Phys. Rev. Lett. 31 (1973) 409; C.A. Tracy and B.M. McCoy, Phys. Rev. Lett. 31 (1973) 1500; T.T. Wu, B.M. McCoy, C.A. Tracy, and E. Barouch, Phys. Rev. B 13 (1976) 316.
  * [11] Orrick, W. P., Nickel, B. G., Guttmann, A. J., and Perk, J. H. H., Phys. Rev. Lett. 86 (2001) 4120; J. Stat. Phys. 102 (2001) 795.
  * [12] M. Caselle, M. Hasenbusch, A. Pelissetto, E. Vicari, J.Phys. A 35 (2002) 4861.
  * [13] G. Delfino, Phys.Lett. B450 (1999) 196
  * [14] G. von Gehlen, V. Rittenberg, and H. Ruegg, J. Phys. A 19 (1985) 107.
  * [15] G. Delfino and J.L. Cardy, Nucl. Phys. B 519 (1998) 551.
  * [16] G. Delfino, G.T. Barkema and J.L. Cardy, Nucl. Phys. B 565 (2000) 521.
  * [17] L.N. Shchur, P. Butera, and B. Berche, Nucl. Phys. B 620 (2002) 579.
  * [18] L.N. Shchur, L.N. Shchur, and P. Butera, Phys. Rev. B 77 (2008) 144410.
  * [19] I.G. Enting and A.J. Guttmann, Physica A 321 (2003) 90.
  * [20] L. Chim and A.B. Zamolodchikov, Int. J. Mod. Phys. A 7 (1992) 5317.
  * [21] M. Nauenberg and D.J. Scalapino, Phys. Rev. Lett. 44 (1980) 837; J. L. Cardy, N. Nauenberg and D.J. Scalapino, Phys. Rev. B 22 (1980) 2560.
  * [22] J. Salas and A. Sokal, J. Stat. Phys. 88 (1997) 567.
  * [23] M. Caselle, R. Tateo, and S. Vinci, Nucl. Phys. B 562 (1999) 549.
  * [24] L.N. Shchur, B. Berche and P. Butera, Europhys. Lett. bf 81 (2008) 30008.
  * [25] L.N. Shchur, B. Berche and P. Butera, unpublished.
  * [26] R. J. Baxter and F.Y. Wu, Phys. Rev. Lett. 31 (1973) 1294; Aust. J. Phys. 27 (1974) 357; R.J. Baxter, Aust. J. Phys. 27 (1974) 369.
  * [27] R.J. Baxter, Exactly Solved Models in Statistical Mechanics (New York, Academic Press, 1982).
  * [28] G.S. Joyce, Proc. R. Soc. Lond. A 343 (1975); ibid 345 (1975) 277.
  * [29] L.N. Shchur and W. Janke, unpublished.
  * [30] G. Delfino and P. Grinza, Nucl. Phys. B 682 (2004) 521.

Figure 2: 3-state Potts model. Ratio of transverse to longitudinal
susceptibilities.


# Performance Evaluation of Mixed-Precision Runge-Kutta Methods

Ben Burnett  University of Massachusetts Dartmouth
<EMAIL_ADDRESS>Sigal Gottlieb  University of Massachusetts Dartmouth
<EMAIL_ADDRESS>Zachary J. Grant  Oak Ridge National Laboratory
<EMAIL_ADDRESS>Alfa Heryudono  University of Massachusetts Dartmouth
<EMAIL_ADDRESS>
###### Abstract

Additive Runge-Kutta methods designed for preserving highly accurate solutions
in mixed-precision computation were proposed and analyzed in [8]. These
specially designed methods use reduced precision for the implicit computations
and full precision for the explicit computations. We develop a FORTRAN code to
solve a nonlinear system of ordinary differential equations using the mixed
precision additive Runge-Kutta (MP-ARK) methods on IBM POWER9 and Intel x86_64
chips. The convergence, accuracy, runtime, and energy consumption of these
methods is explored. We show that these MP-ARK methods efficiently produce
accurate solutions with significant reductions in runtime (and by extension
energy consumption).

###### Index Terms:

Mixed-precision, multiprecision, Runge-Kutta, numerical methods

##  I Introduction

The use of mixed precision to enhance performance of numerical algorithms is
gaining in popularity, as evidenced by the emergence of mixed precision
algorithms for common linear algebra operations[1] to aid in the evaluation of
large computational networks for deep learning. Mixed precision simulations
aim to combine the efficiency of low precision with the accuracy of high
precision computations. However, care is needed here to prevent the low
precision computations from degrading the overall accuracy of the method, or
the high precision computations from adversely impacting the efficiency.

In [2], the authors explored several mixed precision algorithms for solving
sparse linear systems and identified algorithms that benefit from a mixed
precision implementation and others that do not. However, while the
performance of mixed precision algorithms for linear algebra applications has
been studied [2] [3] [4] [5] [6] [7], mixed precision algorithms for solving
ordinary differential equations (ODES) have only recently been proposed. In
[8] Z. Grant proposed a numerical analysis framework for development of mixed-
precision Runge–Kutta methods by using an additive Runge–Kutta method approach
and framing the use of multiple precisions as a perturbation. Using this
rigorous approach, [8] developed novel methods that reduce the cost of the
computationally expensive implicit stages in the Runge–Kutta methods by
employing a low precision computation. Meanwhile, the structure of the MP-ARK
is designed to suppress the low precision errors either by introducing
inexpensive explicit high-precision correction terms, or by designing novel
methods that internally suppress the low precision perturbations.

The work in [8] focused on the development of a framework for analysis of MP-
ARK methods, and the development and convergence verification of such MP-ARK
methods, using a simplified ”chopping” routine to simulate low precision in
MATLAB. In this work, we revisit the novel methods in [8] and evaluate their
performance on multiple CPU architectures in single, double, and quadruple
precision. We first confirm that these methods converge as expected. Next, we
demonstrate that for a given level of accuracy, the mixed precision codes
provide reductions in runtime (and by extension energy consumption) ranging
from 2 times speed up to upwards of 15 times speed up.

##  II MP-ARK methods

For this study we looked at three of the methods outlined in [8]: the mixed
precision implicit midpoint method with and without correction steps, the
mixed-precision SDIRK method with and without correction steps, and a novel
MP-ARK method from [8] designed to suppress the low precision errors without
additional correction steps. In this section we describe these methods.

###  II-A Mixed Precision Methods Based on the Implicit Midpoint Rule

The implicit midpoint method was the main motivating method of the original
paper. It was used as the base example for mixed precision methods and a
demonstration of the perturbed Runge-Kutta analysis to them. The base method
is given by

$\displaystyle y^{(1)}$ $\displaystyle=$ $\displaystyle u^{n}+\frac{\Delta
t}{2}F\big{(}y^{(1)}\big{)}$ $\displaystyle u^{n+1}$ $\displaystyle=$
$\displaystyle u^{n}+\Delta tF\big{(}y^{(1)}\big{)}$

The mixed precision formulation of this method is then defined as

$\displaystyle y^{(1)}$ $\displaystyle=$ $\displaystyle u^{n}+\frac{\Delta
t}{2}F^{\epsilon}\big{(}y^{(1)}\big{)}$ (1a) $\displaystyle u^{n+1}$
$\displaystyle=$ $\displaystyle u^{n}+\Delta tF\big{(}y^{(1)}\big{)}$ (1b)

where $\epsilon$ represents the precision level for the low precision
evaluation.

The low precision evaluation of the implicit solver in the first stage of
Equation (1) requires a subroutine that maps the variables low precision,
performs a Newton iteration with tolerances that are consistent with the level
of the low precision, and returns the low precision value of the variable
$y^{(1)}$. This value is then cast to a high precision variable $y^{(1)}$ to
be used in the explicit second step of Equation (1).

Adding explicit high precision correction steps to the mixed precision
implicit midpoint rule (1) we define the corrected MP-IMR with $p-1$
correction terms

$\displaystyle y_{[0]}^{(1)}$ $\displaystyle=$ $\displaystyle
u^{n}+\frac{\Delta t}{2}F^{\epsilon}\big{(}y_{[0]}^{(1)}\big{)}$ (2a)
$\displaystyle y_{[k]}^{(1)}$ $\displaystyle=$ $\displaystyle
u^{n}+\frac{\Delta
t}{2}F\big{(}y_{[k-1]}^{(1)}\big{)}\;\text{for}\;k=1,...,p-1$ (2b)
$\displaystyle u^{n+1}$ $\displaystyle=$ $\displaystyle u^{n}+\Delta
tF\big{(}y_{[p-1]}^{(1)}\big{)}$ (2c)

The additional correction stages are performed at high precision but they are
cheap explicit computations, while the implicit solves, though performed in
low precision, are expected to be the computationally dominant. The rationale
behind the correction stage approach is that the gain in accuracy will be
well-worth the extra two cheap explicit stage. Of course, the purpose of our
numerical experiments will be to verify this assumption.

###  II-B Mixed Precision SDIRK method with and without corrections

The singly diagonally implicit Runge-Kutta (SDIRK) method is a higher order
than the implicit midpoint method and involves two implicit stages. The mixed
precision SDIRK method with $p-1$ correction terms is (where $k=1:p-1$)

$\displaystyle y_{[0]}^{(1)}$ $\displaystyle=$ $\displaystyle
u^{n}+\gamma\Delta tF^{\epsilon}\big{(}y_{[0]}^{(1)}\big{)}$ $\displaystyle
y_{[k]}^{(1)}$ $\displaystyle=$ $\displaystyle u^{n}+\gamma\Delta
tF\big{(}y_{[k-1]}^{(1)}\big{)}$ (3) $\displaystyle y_{[0]}^{(2)}$
$\displaystyle=$ $\displaystyle y^{n}+(1-2\gamma)\Delta
tF\big{(}y_{[p-1]}^{(1)}\big{)}+\gamma\Delta
tF^{\epsilon}\big{(}y_{[0]}^{(2)}\big{)}$ $\displaystyle y_{[k]}^{(2)}$
$\displaystyle=$ $\displaystyle y^{n}+(1-2\gamma)\Delta
tF\big{(}y_{[p-1]}^{(1)}\big{)}+\gamma\Delta tF\big{(}y_{[k-1]}^{(2)}\big{)}$
$\displaystyle y^{n+1}$ $\displaystyle=$ $\displaystyle u^{n}+\frac{\Delta
t}{2}F\big{(}y_{[p-1]}^{(1)}\big{)}+\frac{\Delta
t}{2}F\big{(}y_{[p-1]}^{(2)}\big{)}.$

###  II-C Mixed precision 4s3p MP-ARK method

The novel mixed precision additive Runge-Kutta method (MP-ARK) developed by Z.
Grant [8] is defined by

$\displaystyle y^{(1)}$ $\displaystyle=$ $\displaystyle u^{n}+\Delta
tA^{\epsilon}_{1,1}F^{\epsilon}\big{(}y^{(1)}\big{)}$ (4a) $\displaystyle
y^{(2)}$ $\displaystyle=$ $\displaystyle u^{n}+\Delta
tA_{2,1}F\big{(}y^{(1)}\big{)}$ (4b) $\displaystyle y^{(3)}$ $\displaystyle=$
$\displaystyle u^{n}+\Delta
t\big{[}A_{3,1}F\big{(}y^{(1)}\big{)}+A_{3,2}F\big{(}y^{(2)}\big{)}\big{]}$
$\displaystyle+$ $\displaystyle\Delta
t\big{[}A^{\epsilon}_{3,1}F^{\epsilon}\big{(}y^{(1)}\big{)}+A^{\epsilon}_{3,3}F^{\epsilon}\big{(}y^{(3)}\big{)}\big{]}$
$\displaystyle y^{(4)}$ $\displaystyle=$ $\displaystyle u^{n}+\Delta
t\big{[}A_{4,1}F\big{(}y_{(1)}\big{)}+A_{4,2}F\big{(}y^{(2)}\big{)}$
$\displaystyle+$ $\displaystyle A_{4,3}F\big{(}y^{(3)}\big{)}\big{]}$
$\displaystyle u^{n+1}$ $\displaystyle=$ $\displaystyle
u^{n}+\frac{1}{2}\Delta
t\big{[}\big{(}y^{(2)}\big{)}+\big{(}y^{(4)}\big{)}\big{]}$ (4e)

where $A$ and $A^{\epsilon}$ are

$\displaystyle A_{2,1}=0.211324865405187,$ $\displaystyle
A_{3,1}=0.709495523817170,$ $\displaystyle A_{3,2}=-0.865314250619423,$
$\displaystyle A_{4,1}=0.705123240545107,$ $\displaystyle
A_{4,2}=0.943370088535775,$ $\displaystyle A_{4,3}=-0.859818194486069,$
$\displaystyle A^{\epsilon}_{1,1}=0.788675134594813,$ $\displaystyle
A^{\epsilon}_{3,1}=0.051944240459852,$ $\displaystyle
A^{\epsilon}_{3,3}=0.788675134594813.$

##  III Implementation Details

The mixed-precision Runge-Kutta methods that we explored were implemented in
Fortran using the 2008 standard. The `iso_fortran_env` intrinsic module was
also used to import the `Real32`, `Real64`, and `Real128` derived types. The
implicit solvers were implemented using the Newton-Raphson method for
multivariate systems.

As in [8], we solve the Van der Pol system

$\displaystyle y_{1}^{\prime}$ $\displaystyle=$ $\displaystyle y_{2}$ (5a)
$\displaystyle y_{2}^{\prime}$ $\displaystyle=$ $\displaystyle
y_{2}(1-y_{1}^{2})-y_{1}$ (5b)

with initial conditions $y_{1}(0)=2.0$, $y_{2}(0)=0.0$, for time $t=[0,1]$. A
reference solution using an explicit fourth order Runge–Kutta at a small time
step was used for calculating the errors. The reference solution was computed
entirely in quadruple precision and all methods results regardless of
computational precision, were cast to quadruple precision in order to compute
the error between a method and this reference solution.

The tolerance of the implicit solvers were set using the Fortran intrinsic
machine epsilon $\epsilon$, where the tolerance value $(1+0.001)\epsilon$ was
chosen. The maximum number of iterations for the solver was set to $20$
regardless of the precision, however this fail-safe maximum was not attained.
These setting do have significant impact on the overall runtime of the method
and were determined to be fair for the consideration of mixed precision
methods, based on our tests showing that the implicit solvers on average only
used $4$ iterations to converge.

Higher precision computations (specifically quadruple precision but also
arbitrary precision) are generally evaluated in software on a majority of
platforms. Notable exceptions for this are the new IBM POWER9 that have
hardware level support for quadruple precision. In this work, we compare the
runtime performance on both software based high precision computations and
hardware based high precision computations. The experiments were performed on
a system with an Intel Xeon Gold 6126 CPU running CentOS 7 and an IBM POWER9
running CentOS 8. We will refer to these as x86 and POWER9, respectively. The
experiment performed looks at three levels of precision: single, double, and
quadruple. The Intel system is limited to hardware evaluation for single and
double precision and requires quadruple precision to be evaluated using
software. The POWER9 system is capable of hardware evaluation of all
precisions. This was done to show the usefulness of these methods under both
conditions and that the benefits are platform independent.

The codes were compiled on both architectures using gfortran version 11.1.0
using no optimization flags. On the POWER9 system, the GCC compilers were
compiled using the IBM AT14.0 libraries in order to enable hardware quadruple
precision support while the GCC compilers on the Intel system were compiled
using the default compiler settings. Timing on the method subroutine was done
using the intrinsic `cpu_time` Fortran subroutine.

The code for this is available on GitLab at
https://gitlab.com/bburnett6/mixed-precision-rk. This repository contains all
the necessary files to reproduce the study performed here. The only
dependencies are a working version of the GCC 11.1.0 compilers and an
environment with the Numpy and MatPlotLib python libraries for the experiment
driver and analysis files that were implemented in python.

##  IV Results

In this section we describe the errors, runtime, and energy consumed when
using the different approaches described in Section II to solve the non-stiff
Van der Pol system (5).

###  IV-A Mixed precision implicit midpoint method

####  IV-A1 Accuracy study

We begin with the mixed precision implicit midpoint rule (1). We reproduced
the convergence behavior reported in [8], as we show in Figure 1. We observe
that the single precision computation is not convergent for many values of
$\Delta t$: the error builds up due to roundoff for smaller $\Delta t$. The
mixed precision computations that use single precision (double-single or quad-
single) perform better, reducing the errors to below $10^{-7}$; however, they
still show evidence of degredation of accuracy due to the single precision
computations. The mixed precision quad-double computation performs as well as
the quad-quad computation for all values of $\Delta t$ tested. However, the
double precision computation performs just as well for all but the smallest
$\Delta t$ tested.

Figure 1: Errors from mixed precision implicit midpoint rule

Figure 2: Errors from mixed precision implicit midpoint rule with one
correction step Figure 3: Total energy used by the mixed precision implicit
midpoint rule with one correction scales linearly with the runtime.

The corrected mixed precision implicit midpoint rule (2) with only one
correction fixes the problems observed in the mixed precision computation of
(1). As we observe in Figure 2, all the mixed precision computations have the
same errors as their high-precision counterparts. An additional correction
(not shown) enhances the errors in the mixed precision simulation a bit
further.

####  IV-A2 Runtime & energy consumption

The energy consumption results for this simulation follow very closely the
runtime results. Figure 3 shows that for long enough runtimes, the energy
consumed scales linearly with the runtime for all precisions. In Figure 4 we
present the runtime of the mixed precision implicit midpoint rule with two
corrections, on the x86 and POWER9 chips, respectively. The runtime is not
impacted by the correction terms.

Figure 4: Mixed precision implicit midpoint method with two corrections: Time
to Completion. Left: x86, Right: POWER9.

Figure 5: Mixed precision implicit midpoint method with one correction step:
error vs runtime. Left: x86; right: POWER9

Remarkable in this figure is the significant reduction in runtime from the
quad precision to the mixed quad-double precision for the x86 chip, and from
the quad precision to the mixed quad-single precision for the POWER9 chip. The
cost of quad precision on the POWER9 chip is significantly less than on the
x86 chips, but the runtime savings is not impacted, as we see in the tables
below. The differences in the runtime behaviors is explained by the fact that
POWER9 chip has hardware-level support for quad precision.

####  IV-A3 Efficiency

To better illustrate the benefit of the mixed precision computation we plot
the method error against the runtime for the mixed precision implicit midpoint
method with one correction, in Figure 5.

Mixed precision IMR method with two corrections
---
Runtime (s) for x86
error | 64/64 | 64/32 | 128/128 | 128/64 | 128/32
$\approx 10^{-9}$ | 0.013 | 0.005 | 0.203 | 0.020 | 0.012
$\approx 10^{-15}$ | N/A | N/A | 108.4 | 13.91 | 7.165

On the x86 (Figure 5, left), we see that the mixed precision double-single
code is the most efficient. For an error level of $\approx 10^{-15}$, the
mixed precision quad-single code has a fifteen-fold reduction in runtime. The
mixed precision double-single code gives the same level of error at $\approx
10^{-9}$ as the double precision code with a runtime reduction of 2.6-fold.

On the POWER9 (Figure 5, right; also table below) the runtime is generally
better, but the runtime savings are similar. For the error level of $\approx
10^{-9}$, we see 3-fold savings from the mixed precision double-single over
the the double precision. For the error level of $\approx 10^{-15}$, we see
more than 14-fold savings from the mixed precision quad-single over the the
quad precision.

Mixed precision IMR method with two corrections
---
Runtime (s) for POWER9
error | 64/64 | 64/32 | 128/128 | 128/64 | 128/32
$\approx 10^{-9}$ | 0.021 | 0.007 | 0.0683 | 0.023 | 0.008
$\approx 10^{-15}$ | N/A | N/A | 34.64 | 14.00 | 2.395

Figure 6: Mixed precision SDIRK method errors vs runtime on x86. Top left: no
corrections; top right: three corrections.

###  IV-B Mixed precision SDIRK Method

The results for the mixed precision singly diagonally implicit Runge-Kutta
(MP-SDIRK) in Section II are similar to those of the implicit midpoint rule.
Figure 7 shows the convergence of this method with no corrections, one
correction, two corrections, and three corrections. The results show clearly
how with each successive correction the errors from the mixed precision
simulation are reduced, until the mixed precision results become as accurate
as their full-precision counterparts. Furthermore, we confirm that additional
corrections do not result in a detrimental increase in runtime.

Figure 7: Mixed precision SDIRK method errors. Top left: no corrections; Top
right: one correction; Bottom left: two corrections; Bottom right: three
corrections.

In Figure 6 we show the errors per runtime for the mixed precision SDIRK
method with no corrections (left) and with three corrections (right), on the
x86 chip. We observe that the method with three corrections is overall
significantly more efficient than the method with no corrections. For errors
above $\approx 10^{-13}$, the mixed precision double-single method performs
best. For errors below that level, the mixed precision quad-single method is
most efficient, as shown in the table below. The runtime savings on the POWER9
chip are similar.

Mixed precision SDIRK with three corrections
---
Runtime (s) for x86
error | 64/64 | 64/32 | 128/128 | 128/64 | 128/32
$\approx 10^{-13}$ | 0.018 | 0.009 | 0.253 | 0.045 | 0.031
$\approx 10^{-22}$ | N/A | N/A | 165.3 | 35.72 | 23.15
Runtime (s) for POWER9
error | 64/64 | 64/32 | 128/128 | 128/64 | 128/32
$\approx 10^{-13}$ | 0.031 | 0.016 | 0.092 | 0.034 | 0.019
$\approx 10^{-22}$ | N/A | N/A | 54.33 | 30.06 | 8.151

###  IV-C Novel Method 4s3pA

Finally, we present the results from the novel four stage third order method
presented in [8] as 4s3pA. The errors per $\Delta t$ are presented in Figure
8. We observe that the mixed precision double-single has essentially the same
errors as the fully double precision computation, at almost 3-fold savings in
runtime. The quad-double has the same errors as the fully quad precision
computation at runtime savings factors of 2.6 to 5.75 depending on the
platform.

Figure 8: Novel mixed-precision method 4s3pA Method Error

Figure 9: Novel mixed-precision method 4s3pA Error vs runtime

Mixed precision novel 4s3pA method
---
Runtime (s) for x86
error | 64/64 | 64/32 | 128/128 | 128/64 | 128/32
$\approx 10^{-13}$ | 0.026 | 0.009 | 0.4300 | 0.052 | 0.035
$\approx 10^{-22}$ | N/A | N/A | 231.4 | 40.71 | N/A
Runtime (s) for POWER9
error | 64/64 | 64/32 | 128/128 | 128/64 | 128/32
$\approx 10^{-13}$ | 0.047 | 0.017 | 0.140 | 0.049 | 0.018
$\approx 10^{-22}$ | N/A | N/A | 74.68 | 29.42 | N/A

##  V Conclusions

In this paper we have presented a numerical study of the new mixed precision
Runge–Kutta methods on two different platforms. The results presented confirm
that the mixed precision Runge-Kutta methods on a system of nonlinear ODEs
provide significant runtime (and energy) savings per level of error. This is
of particular interest when a high level of accuracy is required, so that
expensive high precision must be used. The nonlinear ODE system explored here
is a small and relatively simple problem; we expect that these savings would
be extended (and indeed improved) on more complicated problems.

## Acknowledgment

Z. Grant acknowledges the support of the Department of Computational and
Applied Mathematics, Oak Ridge National Laboratory, Oak Ridge TN 37830. This
material is based upon work supported by the U.S. Department of Energy, Office
of Science, Office of Advanced Scientific Computing Research, as part of their
Applied Mathematics Research Program. The work was performed at the Oak Ridge
National Laboratory, which is managed by UT-Battelle, LLC under Contract No.
De- AC05-00OR22725. The United States Government retains and the publisher, by
accepting the article for publication, acknowledges that the United States
Government retains a non-exclusive, paid-up, irrevocable, world-wide license
to publish or reproduce the published form of this manuscript, or allow others
to do so, for the United States Government purposes. The Department of Energy
will provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan http://energy.gov/downloads/doe-
public-access-plan.

Gottlieb’s work was partially supported by AFOSR grant FA9550-18-1-0383.
Burnett, Gottlieb, and Heryudono’s work were partially supported by ONR UMass
Dartmouth Marine and UnderSea Technology (MUST) grant N00014-20-1-2849 under
the project S31320000049160.

All computations for this project were performed on CARNIE, the UMass-
Dartmouth HPC cluster purchased by ONR DURIP grant N00014-18-1-2255.

We are grateful to the IBM Linux team that extended a warm and helping hand to
figure out how our POWER9 system was running, and what will be needed to get
hardware quadruple precision support.

## References

  * [1] Suyog Gupta, Ankur Agrawal, Kailash Gopalakrishnan, and Pritish Narayanan. Deep Learning with Limited Numerical Precision. In Proceedings of the 32Nd International Conference on International Conference on Machine Learning - Volume 37, ICML’15, pages 1737–1746. JMLR.org, 2015.
  * [2] Mawussi Zounon, Nicholas J. Higham, Craig Lucas, and Françoise Tisseur. 2020. Performance Evaluation of Mixed Precision Algorithms for Solving Sparse Linear Systems. 1, 1 (September 2020), 18 pages. https://doi.org/10.1145/nnnnnnn.nnnnnnn.
  * [3] Azzam Haidar, Stanimire Tomov, Jack Dongarra, and Nicholas J. Higham. Harnessing GPU Tensor Cores for Fast FP16 Arithmetic to Speed Up Mixed-precision Iterative Refinement Solvers. In Proceedings of the International Conference for High Performance Computing, Networking, Storage, and Analysis, SC ’18, pages 47:1–47:11, Piscataway, NJ, USA, 2018. IEEE Press.
  * [4] Christian Richter, Sebastian Schops, and Markus Clemens. Gpu-accelerated mixed precision algebraic multigrid preconditioners for discrete elliptic field problems. IEEE Transactions on Magnetics, 50(2), 2014.
  * [5] Ahmad Abdelfattah, Stanimire Tomov, and Jack J. Dongarra. Fast Batched Matrix Multiplication for Small Sizes Using Half-Precision Arithmetic on GPUs. In 2019 IEEE International Parallel and Distributed Processing Symposium, IPDPS 2019, Rio de Janeiro, Brazil, May 20-24, 2019, pages 111–122. IEEE, 2019.
  * [6] Nicholas J. Higham. Error analysis for standard and GMRES-based iterative refinement in two and three-precisions. MIMS EPrint 2019.19, Manchester Institute for Mathematical Sciences, The University of Manchester, November 2019.
  * [7] Petschow, M. and Quintana-Ort, E. S. and Bientinesi, P. Improved Accuracy and Parallelism for MRRR-Based Eigensolvers—A Mixed Precision Approach. SIAM Journal on Scientific Computing, 36(2):C240–C263, 2014. https://doi.org/10.1137/130911561.
  * [8] Zachary J. Grant, “Perturbed Runge-Kutta Methods for Mixed Precision Applications,” https://arxiv.org/pdf/2012.13055.pdf, 24 December, 2020.