#  On off-critical zeros of lattice energies in the neighborhood of the
Riemann zeta function

Laurent Bétermin  Institut Camille Jordan, Université Claude Bernard Lyon 1,
69622 Villeurbanne, France  Ladislav Šamaj  Igor Travěnec  Institute of
Physics, Slovak Academy of Sciences, Dúbravská cesta 9, 84511 Bratislava,
Slovakia

###### Abstract

The Riemann zeta function $\zeta(s):=\sum_{n=1}^{\infty}1/n^{s}$ can be
interpreted as the energy per point of the lattice $\mathbb{Z}$, interacting
pairwisely via the Riesz potential $1/r^{s}$. Given a parameter
$\Delta\in(0,1]$, this physical model is generalized by considering the energy
per point $E(s,\Delta)$ of a periodic one-dimensional lattice alternating the
distances between the nearest-neighbour particles as $2/(1+\Delta)$ and
$2\Delta/(1+\Delta)$, keeping the lattice density equal to one independently
of $\Delta$. This energy trivially satisfies $E(s,1)=\zeta(s)$ at $\Delta=1$,
it can be easily expressed as a combination of the Riemann and Hurwitz zeta
functions, and extended analytically to the punctured $s$-plane
$\mathbb{C}\setminus\\{1\\}$. In this paper, we perform numerical
investigations of the zeros of the energy $\\{\rho=\rho_{x}+{\rm
i}\rho_{y}\\}$, which are defined by $E(\rho,\Delta)=0$. The numerical results
reveal that in the Riemann limit $\Delta\to 1^{-}$ theses zeros include the
anticipated critical zeros of the Riemann zeta function with
$\Re(\rho_{x})=\frac{1}{2}$ as well as an unexpected – comparing to the
Riemann Hypothesis – infinite series of off-critical zeros. The analytic
treatment of these off-critical zeros shows that their imaginary components
are equidistant and their real components diverge logarithmically to $-\infty$
as $\Delta\to 1^{-}$, i.e., they become invisible at the Riemann’s $\Delta=1$.

###### keywords:

Riemann zeta function; Hurwitz zeta function; critical and off-critical zeros;
Riemann hypothesis

###### MSC:

[2010]11E45

##  1 Introduction and main results

Let two points at distance $r$ interact via the Riesz potential $1/r^{s}$ with
real $s$ [4]. If the points are located on the lattice $\mathbb{Z}$ and
interact pairwisely by the Riesz potential where $s>1$, the energy per point
is given by the Riemann zeta function [20]

$\zeta(s):=\frac{1}{2}\sum_{n\in\mathbb{Z}^{*}}\frac{1}{|n|^{s}}=\sum_{n=1}^{\infty}\frac{1}{n^{s}}\qquad
s>1,$ (1)

where the prefactor $\frac{1}{2}$ is due to the fact that each interaction
energy is shared by a pair of points. The function $\zeta$ can be analytically
continued to the whole complex $s$-plane, with a simple pole at $s=1$. The
Riemann zeta function plays a fundamental role in the algebraic and analytic
number theories [12, 13, 21, 14, 16, 24, 22], see monographs [7, 17, 25]. The
so-called Riemann Hypothesis about the location of its nontrivial zeros
exclusively on the critical line $\Re(s)=\frac{1}{2}$ (the symbol $\Re$ means
the real part) is one of the Hilbert and Clay Millennium Prize problems [18].
Throughout the present paper we assume that the Riemann Hypothesis holds. The
Riemann zeta function and its Epstein’s [9, 10, 5, 26], Hurwitz’s [15, 23,
19], Barnes’s [2], etc. generalisations have numerous applications both in
mathematics (prime numbers, applied statistics [3]) and in physics [8].

Let the Riemann zeta function be a member of a family of functions which
exhibit nontrivial zeros off the critical line. Possible mechanisms of the
disappearance of these off-critical zeros at the Riemann’s point might of
general interest. In this paper, we propose a natural extension of the Riemann
zeta function as the energy of a unit density lattice $L_{\Delta}$ with
alternating distances between the nearest neighbours, say $2/(1+\Delta)$ and
$2\Delta/(1+\Delta)$; due to the $\Delta\to 1/\Delta$ symmetry of the problem,
it is sufficient to restrict oneself to $\Delta$ from the interval $(0,1]$.

In analogy with the original model with constant unit spacing, each point
interacts pairwisely with the other points via the Riesz interaction
$1/r^{s}$, $s>1$ and the lattice energy per point is therefore given (see
Proposition 2.1) by

$E(s,\Delta)=\frac{1}{2^{s}}\zeta(s)+\frac{1}{2^{s+1}}\left[\zeta\left(s,\frac{1}{1+\Delta}\right)+\zeta\left(s,\frac{\Delta}{1+\Delta}\right)\right],\
s>1,\Delta\in(0,1],$

where

$\zeta(s,a)=\sum_{n=0}^{\infty}\frac{1}{(n+a)^{s}},\qquad s>1,$ (2)

is the Hurwitz zeta function with the real (positive) parameter $a$. Remark
that this lattice energy, as a combination of Riemann and Hurwitz zeta
functions, has an analytic continuation on $\mathbb{C}\backslash\\{1\\}$ (see
Proposition 2.3). For given $\Delta\in(0,1]$, the set of zeros of the lattice
energy is defined as $Z_{\Delta}:=\\{\rho=\rho_{x}+{\rm
i}\rho_{y}\in\mathbb{C},(\rho_{x},\rho_{y})\in\mathbb{R}^{2}:E(\rho,\Delta)=0\\}$,
noticing that $Z_{1}$ is the set of zeros of the Riemann zeta function
composed by trivial zeros (i.e., $\rho\in-2\mathbb{N}$) and critical zeros
(i.e., $\Re(\rho)=\frac{1}{2}$) assuming that Riemann Hypothesis holds.
Furthermore, for specific values of the parameter
$\Delta\in\\{1/5,1/3,1/2\\}$, the energy can be factorized as
$E(s,\Delta)=f_{\Delta}(s)\zeta(s)$ where $f_{\Delta}$ is a sum of $p^{s}$
with integers $p$. This automatically gives us critical zeros (i.e. solutions
of $\zeta(\rho)=0$ assuming the Riemann Hypothesis) and possible off-critical
zeros (i.e. solutions of $f_{\Delta}(\rho)=0$), as shown in Proposition 2.4.

The goal of this paper is to study, both numerically and analytically, the set
of zeros $Z_{\Delta}$ when $\Delta$ is in a neighborhood of $1$, i.e. when
$E(\cdot,\Delta)$ is in the neighborhood of the Riemann zeta function.
Numerical and analytic analysis shows that approaching $\Delta\to 1^{-}$ the
zeros of $E(\rho,\Delta)$ involve the anticipated critical zeros of the
Riemann zeta function with $\Re(\rho_{x})=\frac{1}{2}$ as well as an infinite
series of unexpected off-critical zeros with the following asymptotics for
their components, as $\Delta\to 1^{-}$ (see Theorem 4.7):

$\displaystyle\rho_{x}(\Delta)=\frac{2}{\ln
2}\ln(1-\Delta)+O(1-\Delta)\to-\infty,$
$\displaystyle\rho_{y}(\Delta)=\frac{(2k+1)\pi}{\ln
2}+O\left((1-\Delta)^{\frac{2\ln 3}{\ln 2}}\right),\quad k\in\mathbb{Z}.$

This means that, asymptotically, there is a infinite sequence of equidistant
zero components along the $\rho_{y}$ axis. Furthermore, the divergence of
$\rho_{x}$ to $-\infty$ as $\Delta\to 1^{-}$ is an example of the
disappearance of off-critical zeros when approaching the Riemann’s point.
Moreover, the behavior of these zero components with respect to
$\Delta\in(0,1]$ is numerically studied (see Figures 1 and 2).

Plan of the paper. The generalized 1D model for Riesz points with alternating
lattice spacings is presented in section 2. The energy per particle
$E(s,\Delta)$ is expressed as a combination of Hurwitz zeta functions in
section 2.1. The properties of the Hurwitz zeta function are discussed in
section 2.2. Special values of the parameters $\Delta$ when the energy
$E(s,\Delta)$ factorizes itself onto the product of the Riemann zeta function
and some simple function are given in section 2.3. Numerical results for zeros
at any $0<\Delta<1$, together with tests at the special values of
$\Delta=1/5,1/3,1/2$ are presented in section 3. The spectrum of critical and
off-critical zeros in the Riemann’s limit $\Delta\to 1^{-}$ is discussed in
section 4.

##  2 The generalized one-dimensional model

###  2.1 Definition of the model

Given $\Delta\in(0,1]$, we consider the infinite set of points
$L_{\Delta}\subset\mathbb{R}$ given by

$L_{\Delta}:=2\mathbb{Z}\cup\left(2\mathbb{Z}+\frac{2\Delta}{1+\Delta}\right),$

which is the unit density periodic configuration with alternating distances
$2/(1+\Delta)$ and $2\Delta/(1+\Delta)$, since
$2-2/(1+\Delta)=2\Delta/(1+\Delta)$. Assuming that each pair of points in
$L_{\Delta}$ interacts via the Riesz potential $1/r^{s}$, $s>1$, the total
energy per point of this system is therefore

$E(s,\Delta):=\frac{1}{4}\sum_{k\in\left\\{0,\frac{2\Delta}{1+\Delta}\right\\}}\sum_{p\in
L_{\Delta}\atop p\neq k}\frac{1}{|p-k|^{s}}.$

The following proposition shows how to write this energy in terms of Riemann
and Hurwitz zeta functions.

######  Proposition 2.1.

For any $s>1$ and any $\Delta\in(0,1]$, we have

$E(s,\Delta)=\frac{1}{2^{s}}\zeta(s)+\frac{1}{2^{s+1}}\left[\zeta\left(s,\frac{1}{1+\Delta}\right)+\zeta\left(s,\frac{\Delta}{1+\Delta}\right)\right].$
(3)

###### Proof.

We simply compute the above double sum as follows:

$\displaystyle E(s,\Delta):$
$\displaystyle=\frac{1}{4}\sum_{k\in\left\\{0,\frac{2\Delta}{1+\Delta}\right\\}}\sum_{p\in
L_{\Delta}\atop p\neq k}\frac{1}{|p-k|^{s}}$
$\displaystyle=\frac{1}{4}\sum_{k\in\left\\{0,\frac{2\Delta}{1+\Delta}\right\\}}\left(\sum_{p\in
2\mathbb{Z}\atop p\neq k}\frac{1}{|p-k|^{s}}+\sum_{p\in
2\mathbb{Z}+\frac{2\Delta}{1+\Delta}\atop p\neq k}\frac{1}{|p-k|^{s}}\right)$
$\displaystyle=\frac{1}{4}\sum_{n\in\mathbb{Z}\atop n\neq
0}\frac{1}{|2n|^{s}}+\frac{1}{4}\sum_{n\in\mathbb{Z}\atop
2n\neq\frac{2\Delta}{1+\Delta}}\frac{1}{\left|2n-\frac{2\Delta}{1+\Delta}\right|^{s}}$
$\displaystyle+\frac{1}{4}\sum_{n\in\mathbb{Z}\atop
2n\neq-\frac{2\Delta}{1+\Delta}}\frac{1}{\left|2n+\frac{2\Delta}{1+\Delta}\right|^{s}}+\frac{1}{4}\sum_{n\in\mathbb{Z}\atop
n\neq 0}\frac{1}{|2n|^{s}}$
$\displaystyle=\frac{1}{2^{s+1}}\sum_{n\in\mathbb{Z}^{*}}\frac{1}{|n|^{s}}+\frac{1}{2^{s+2}}\sum_{n\in\mathbb{Z}}\frac{1}{\left|n+\frac{\Delta}{1+\Delta}\right|^{s}}+\frac{1}{2^{s+2}}\sum_{n\in\mathbb{Z}}\frac{1}{\left|n-\frac{\Delta}{1+\Delta}\right|^{s}}.$

We now split the two last sums in order to get two Hurwitz zeta functions and
two rests that we write again in terms of the same Hurwitz zeta functions:

$\displaystyle E(s,\Delta):$
$\displaystyle=\frac{1}{2^{s}}\zeta(s)+\frac{1}{2^{s+2}}\left[\zeta\left(s,\frac{1}{1+\Delta}\right)+\zeta\left(s,\frac{\Delta}{1+\Delta}\right)\right]$
$\displaystyle\qquad+\frac{1}{2^{s+2}}\left(\sum_{n=1}^{+\infty}\frac{1}{\left|-n+\frac{\Delta}{1+\Delta}\right|^{s}}+\sum_{n=1}^{+\infty}\frac{1}{\left|-n+\frac{1}{1+\Delta}\right|^{s}}\right)$
$\displaystyle=\frac{1}{2^{s}}\zeta(s)+\frac{1}{2^{s+2}}\left[\zeta\left(s,\frac{1}{1+\Delta}\right)+\zeta\left(s,\frac{\Delta}{1+\Delta}\right)\right]$
$\displaystyle\qquad+\frac{1}{2^{s+2}}\left(\sum_{n=1}^{+\infty}\frac{1}{\left(n-\frac{\Delta}{1+\Delta}\right)^{s}}+\sum_{n=1}^{+\infty}\frac{1}{\left(n-\frac{1}{1+\Delta}\right)^{s}}\right).$

Since we have, by the change of variables $n=k+1$,

$\sum_{n=1}^{+\infty}\frac{1}{\left(n-\frac{\Delta}{1+\Delta}\right)^{s}}=\sum_{k=0}^{+\infty}\frac{1}{\left(k+\frac{1}{1+\Delta}\right)^{s}}\
\textnormal{and}\
\sum_{n=1}^{+\infty}\frac{1}{\left(n-\frac{1}{1+\Delta}\right)^{s}}=\sum_{k=0}^{+\infty}\frac{1}{\left(k+\frac{\Delta}{1+\Delta}\right)^{s}},$

we obtain

$\displaystyle E(s,\Delta)$
$\displaystyle=\frac{1}{2^{s}}\zeta(s)+\frac{2}{2^{s+2}}\left[\zeta\left(s,\frac{1}{1+\Delta}\right)+\zeta\left(s,\frac{\Delta}{1+\Delta}\right)\right]$
$\displaystyle=\frac{1}{2^{s}}\zeta(s)+\frac{1}{2^{s+1}}\left[\zeta\left(s,\frac{1}{1+\Delta}\right)+\zeta\left(s,\frac{\Delta}{1+\Delta}\right)\right]$

and the proof is complete. ∎

Notice that the energy satisfies the required symmetry relation
$E(s,\Delta)=E(s,1/\Delta)$.

######  Remark 2.2 (Crystallization result).

It is known (see e.g. [28]), by a convexity argument (or by the so-called one-
dimensional “Universal Optimality” of $\mathbb{Z}$, see [6]) that, for all
$s>0$,

$\min_{\Delta\in(0,1]}E(s,\Delta)=E(s,1)=\zeta(s),$

with equality if and only if $\Delta=1$. From our results (see Theorem 4.7),
the Riemann zeta function is therefore at the same time the minimal value of
our energy and the only one in its $\Delta$-neighborhood for which the non-
trivial zeros are strictly located on the critical line $\Re(s)=\frac{1}{2}$.
It might be interesting to investigate other lattice energies to understand
how universal this phenomenon is.

###  2.2 The Hurwitz zeta function and the analytic continuation of
$E(\cdot,\Delta)$

The Hurwitz zeta function (2) is a generalization of the Riemann zeta-function
(1) via a shift $a>0$. In particular,

$\zeta(s,1)=\zeta(s).$ (4)

In the symbolic computer language Mathematica, the Riemann and Hurwitz zeta
functions are tabulated under the symbols Zeta[$s$] and Zeta[$s,a$],
respectively.

The Hurwitz zeta function satisfies two easily verifiable important
equalities,

$\forall x\in\left[0,\frac{1}{2}\right],\forall
s>1,\quad\zeta(s,x)+\zeta(s,1/2+x)=2^{s}\zeta(s,2x),$ (5)

and the multiplication theorem

$k\in\mathbb{N},\forall s>1,\quad
k^{s}\zeta(s)=\sum_{n=1}^{k}\zeta\left(s,\frac{n}{k}\right).$ (6)

The second relation easily implies that

$\zeta(s,1/3)+\zeta(s,2/3)=(3^{s}-1)\zeta(s)$ (7)

as well as

$\zeta(s,1/2)=(2^{s}-1)\zeta(s).$ (8)

From (5) with $x=1/4$ and (8), we therefore obtain

$\zeta(s,1/4)+\zeta(s,3/4)=(4^{s}-2^{s})\zeta(s).$ (9)

From (3) and (8), it is also straightforward to check that $E(s,1)=\zeta(s)$.

Furthermore, it is clear that (3) holds for all $s\in\mathbb{C}$ such that
$\Re(s)>1$ and we get the following result by directly applying the classical
one (see e.g. [11, 15, 26]) on the analytic continuation of the Riemann and
Hurwitz zeta functions.

######  Proposition 2.3.

For all $\Delta\in(0,1]$, the function $s\mapsto E(s,\Delta)$ has an analytic
continuation on $\mathbb{C}\backslash\\{1\\}$. Furthermore, we have, for all
$s$ such that $\Re(s)<1$ and all $\Delta\in(0,1]$,

$\displaystyle\pi^{-\frac{s}{2}}\Gamma\left(\frac{s}{2}\right)E(s,\Delta)$
$\displaystyle=\frac{1}{2^{s+1}}\left\\{\int_{0}^{\infty}\left[\vartheta\left(\frac{1}{1+\Delta},it\right)-1\right]t^{\frac{-1-s}{2}}{\rm
d}t+f(s)\right\\},$

where $\displaystyle\vartheta(z,it)=\sum_{n\in\mathbb{Z}}e^{-\pi
n^{2}t}e^{2i\pi nz}$ is the Jacobi theta function defined for $t>0$ and
$z\in\mathbb{C}$, and

$\displaystyle f(s)$
$\displaystyle=\int_{0}^{\infty}\left[\vartheta(0,it)-1-\frac{1}{\sqrt{t}}\right]t^{\frac{s}{2}-1}{\rm
d}t\qquad\mbox{for $0<\Re(s)<1$,}$
$\displaystyle=\int_{0}^{\infty}\left[\vartheta(0,it)-\frac{1}{\sqrt{t}}\right]t^{\frac{s}{2}-1}{\rm
d}t\qquad\mbox{for $\Re(s)<0$.}$

###### Proof.

Recall that, for $\Delta\in(0,1]$ and $\Re(s)>1$, we have

$E(s,\Delta)=\frac{1}{2^{s}}\zeta(s)+\frac{1}{2^{s+1}}\left[\zeta\left(s,\frac{1}{1+\Delta}\right)+\zeta\left(s,\frac{\Delta}{1+\Delta}\right)\right].$

It has been shown in [11, 15, 26] that $s\mapsto\zeta(s)$ and
$s\mapsto\zeta(s,a)$, $a>0$, admit an analytic continuation to
$\mathbb{C}\backslash\\{1\\}$, which implies the same for $s\mapsto
E(s,\Delta)$. Furthermore, writing $z=\frac{1}{1+\Delta}$, we have

$E(s,\Delta)=\frac{1}{2^{s}}\zeta(s)+\frac{1}{2^{s+1}}\left[\zeta\left(s,z\right)+\zeta\left(s,1-z\right)\right].$

Considering the analytic continuation of $\zeta(1-\alpha,a)$, $a>0$, the
following formula is well-known [11, 15] for all $\alpha$ such that
$\Re(\alpha)>0$:

$\pi^{-\frac{1-\alpha}{2}}\Gamma\left(\frac{1-\alpha}{2}\right)\left[\zeta\left(1-\alpha,z\right)+\zeta\left(1-\alpha,1-z\right)\right]=\int_{0}^{\infty}[\vartheta(z,it)-1]t^{\frac{\alpha}{2}}\frac{{\rm
d}t}{t}$ (10)

holding for $z\notin\mathbb{Z}$ (see the remark in Theorem 12.6 on page 257 of
[1]). Replacing $\alpha$ by $s=1-\alpha$ (so that $\Re(s)<1$ for
$\Re(\alpha)>0$), we write

$\displaystyle\pi^{-\frac{s}{2}}\Gamma\left(\frac{s}{2}\right)E(s,\Delta)$
$\displaystyle=\frac{\pi^{-\frac{s}{2}}}{2^{s}}\Gamma\left(\frac{s}{2}\right)\zeta(s)+\frac{\pi^{-\frac{s}{2}}}{2^{s+1}}\Gamma\left(\frac{s}{2}\right)\left[\zeta\left(s,z\right)+\zeta\left(s,1-z\right)\right]$
$\displaystyle=\frac{2\pi^{-\frac{s}{2}}}{2^{s+1}}\Gamma\left(\frac{s}{2}\right)\zeta(s)+\frac{\pi^{-\frac{s}{2}}}{2^{s+1}}\Gamma\left(\frac{s}{2}\right)\left[\zeta\left(s,z\right)+\zeta\left(s,1-z\right)\right].$

Next we replace
$\pi^{-\frac{s}{2}}\Gamma\left(\frac{s}{2}\right)\left[\zeta\left(s,z\right)+\zeta\left(s,1-z\right)\right]$
by the integral given in (10) and
$2\pi^{-\frac{s}{2}}\Gamma\left(\frac{s}{2}\right)\zeta(s)$ by the $d=1$
integral in Eq. (21) of [26] for $0<\Re(s)<1$ and by the $d=1$ integral in Eq.
(22) of [26] for $\Re(s)<0$ to complete the proof. ∎

Therefore, we can consider the zeros of $s\mapsto E(s,\Delta)$ in
$\mathbb{C}\backslash\\{1\\}$ defined as

$Z_{\Delta}:=\\{\rho=\rho_{x}+{\rm
i}\rho_{y}\in\mathbb{C},(\rho_{x},\rho_{y})\in\mathbb{R}^{2}:E(\rho,\Delta)=0\\},$

noticing that $Z_{1}$ is the set of zeros of the Riemann zeta function. We
recall that, according to the Riemann Hypothesis,

$Z_{1}=-2\mathbb{N}\cup Z^{C},\quad Z^{C}\subset\\{\Re(z)=1/2\\},$

where $Z^{C}$ is called the set of critical zeros of $\zeta$ and
$-2\mathbb{N}$ are the trivial zeros of $\zeta$. We are going to see in the
next sections that $Z_{\Delta}$ can have other nontrivial off-critical zeros.

###  2.3 Factorization and zeros of the energy for special values of $\Delta$

There exist special values of $\Delta$ for which the energy $E(s,\Delta)$
factorizes itself into a product of the Riemann zeta function $\zeta(s)$ and
some simple functions of $s$, by using the previously presented relations (5)
and (7). For these cases, both critical and off-critical zeros can be found
easily. The most obvious choice of $\Delta$ is $\Delta=1$ for which we have

$E(s,1)=\zeta(s).$ (11)

In the cases $\Delta\in\left\\{\frac{1}{2},\frac{1}{3}\right\\}$ we have the
following result giving the zeros of $E(s,\Delta)$ as well as the
factorization of the energy.

######  Proposition 2.4.

For all $s\in\mathbb{C}\backslash\\{1\\}$, we have

$E(s,1/2)=\frac{1}{2^{s+1}}(1+3^{s})\zeta(s)\quad\textnormal{and}\quad
E(s,1/3)=\frac{1}{2^{s+1}}(2-2^{s}+4^{s})\zeta(s).$ (12)

Furthermore, the zeros of $E(s,1/2)$ and $E(s,1/3)$ are

$\displaystyle Z_{\frac{1}{2}}=Z_{1}\cup\left\\{\frac{(2k+1){\rm i}\pi}{\ln
3}\right\\}_{k\in\mathbb{Z}},$ $\displaystyle
Z_{\frac{1}{3}}=Z_{1}\cup\left\\{\frac{1}{\ln{2}}\left[\ln\left(\frac{1\pm{\rm
i}\sqrt{7}}{2}\right)+2\pi{\rm i}k\right]\right\\}_{k\in\mathbb{Z}}.$

###### Proof.

For $\displaystyle\Delta=\frac{1}{2}$, one has
$\displaystyle\frac{1}{1+\Delta}=\frac{2}{3}$ and
$\displaystyle\frac{\Delta}{1+\Delta}=\frac{1}{3}$ and therefore, using (7),
it holds that

$E(s,1/2)=\frac{1}{2^{s+1}}(1+3^{s})\zeta(s).$ (13)

where the function $1+3^{s}$ yields an infinite sequence of (purely imaginary)
off-critical zeros

$\rho_{k}=\frac{(2k+1){\rm i}\pi}{\ln{3}},\quad k\in\mathbb{Z}.$ (14)

Furthermore, for $\displaystyle\Delta=\frac{1}{3}$, we have
$\displaystyle\frac{1}{1+\Delta}=\frac{3}{4}$ and
$\displaystyle\frac{\Delta}{1+\Delta}=\frac{1}{4}$ and therefore, applying
(9), it holds that

$E(s,1/3)=\frac{1}{2^{s+1}}(2-2^{s}+4^{s})\zeta(s).$ (15)

The function $2-2^{s}+4^{s}$ yields the following zeros

$\rho_{k}=\frac{1}{\ln{2}}\left[\ln\left(\frac{1\pm{\rm
i}\sqrt{7}}{2}\right)+2\pi{\rm i}k\right],\quad k\in\mathbb{Z}.$ (16)

∎

######  Remark 2.5.

In the $\Delta=1/3$ case, since all these zeros have the real part equal to
$\frac{1}{2}$, the energy for $\Delta=1/3$ exhibits only critical zeros.

The last factorization we are considering in our paper corresponds to
$\Delta=1/5$:

$E(s,1/5)=\frac{1}{2^{s+1}}(3-2^{s}-3^{s}+6^{s})\zeta(s).$ (17)

The function $3-2^{s}-3^{s}+6^{s}$ exhibits only off-critical zeros which can
be found only numerically, e.g., $s\approx 0.635084\pm 1.07885{\rm i}$.

##  3 Numerical results

The starting point of our numerical determination of zeros of the energy
$E(s,\Delta)$ was the case $\Delta=1/2$, with the factorization form (12),
whose spectrum of zeros involves both the critical zeros of the Riemann zeta
function as well as an infinite set of off-critical zeros (14). It was checked
that the accuracy of determination of complex zeros by using the symbolic
language Mathematica is 34-35 decimal digits for both real and imaginary
components.

Then we proceeded to the left and right from this point $\Delta=1/2$ by
changing successively $\Delta$ by a small amount to avoid an uncontrolled skip
between neighbouring branches of zeros. Our numerical experience indicates
that changing $\Delta$ by 0.01 is certainly safe from this point of view. We
observe the following:

  * 1.

The zeros of $E(s,\Delta)$ form continuous non-crossing curves in the complete
$(\rho_{x},\rho_{y},\Delta)$ space.

  * 2.

Nevertheless, the curves may intersect in the reduced spaces
$(\rho_{y},\Delta)$ and $(\rho_{x},\Delta)$, see Figures 1 and 2,
respectively.

Notice that the values of $\Delta\in\\{1/5,1/3,1/2\\}$ serve as test points of
our numerical calculations.

Figure 1: Imaginary parts $\rho_{y}$ of the energy zeros $E(\rho,\Delta)=0$.
The special cases $\Delta\in\\{1/5,1/3,1/2,1\\}$ are visualized by vertical
dashed lines. The points which correspond to critical zeros with
$\rho_{x}=\frac{1}{2}$ are denoted by red colour.

The dependence of the imaginary component of zeros $\rho_{y}$, in the range of
its values $[0,25]$, on the parameter $\Delta\in(0,1]$ is pictured in Figure
1. The special cases $\Delta\in\\{1/5,1/3,1/2\\}$, when the energy factorizes
itself onto a product of the Riemann zeta function and a simple function, are
visualized by vertical dashed lines. These cases yield us precise values of
the corresponding zeros and also make us sure not to miss any curve of zeros.
The critical zeros with $\rho_{x}=\frac{1}{2}$ are denoted by red colour, all
other zeros are off-critical; the off-critical zeros for the special values of
$\Delta\in\\{1/5,1/2\\}$ are denoted by blue colour. We observe two types of
zero curves:

  * 1.

The first three “standard” zero curves $\rho_{y}(\Delta)$, which end up at the
Riemann critical zeros at $\Delta=1$, are represented by full symbols.

  * 2.

The first three “non-standard” curves $\rho_{y}(\Delta)$ are represented by
open symbols. Since in the limit $\Delta\to 1^{-}$ these curves tend to off-
critical zeros with the divergent real component $\rho_{x}\to-\infty$, the
curves end up with crosses indicating the absence of off-critical zeros at the
Riemann’s $\Delta=1$.

Figure 2: The dependence of the real part $\rho_{x}$ of zeros
$E(\rho,\Delta)=0$, corresponding to the three “nonstandard” curves denoted by
open symbols in Figure 1, on the parameter $\Delta\in[0,1]$. The critical
zeros with $\rho_{x}=\frac{1}{2}$ are denoted by red colour.

The dependence of the real part $\rho_{x}$ of the first three “non-standard”
energy zeros on $\Delta\in(0,1]$ is presented in Figure 2 by open symbols, in
close analogy with Figure 1. As before, the critical zeros with
$\rho_{x}=\frac{1}{2}$ are denoted by red colour. We observe the following:

  * 1.

For $\Delta\gtrapprox 0.75$, the three curves coincide on the considered scale
and go to $-\infty$ as $\Delta\to 1^{-}$.

##  4 Analytic results in the limit $\Delta\to 1^{-}$

Approaching $\Delta\to 1^{-}$, with regard to Eq. (11) one anticipates the
presence of critical zeros of the Riemann zeta function for $E(s,\Delta)$.
Surprisingly, as was already indicated, there are also additional curves of
off-critical zeros.

To derive coordinates of these off-critical zeros, we set
$\Delta=1-\varepsilon$ in (3) and expand the energy in Taylor series in the
small positive $\varepsilon\to 0^{+}$ up to the order $\varepsilon^{5}$.

######  Proposition 4.6.

Let $s\in\mathbb{C}\backslash\\{1\\}$, then, as $\varepsilon\to 0^{+}$,

$\displaystyle E(s,1-\varepsilon)$
$\displaystyle=\zeta(s)+\frac{2^{2+s}-1}{2^{5+s}}s(1+s)\zeta(2+s)\left(\varepsilon^{2}+\varepsilon^{3}+\frac{3}{4}\varepsilon^{4}+\frac{1}{2}\varepsilon^{5}\right)$
$\displaystyle+\frac{1}{3}\frac{2^{4+s}-1}{2^{11+s}}s(1+s)(2+s)(3+s)\zeta(4+s)\left(\varepsilon^{4}+2\varepsilon^{5}\right)+O(\varepsilon^{6}).$

###### Proof.

It directly follows from the Taylor expansion of the Hurwitz zeta function
(see e.g. [27]): for $|a|<1$,

$\zeta(s,a)=\frac{1}{a^{s}}+\sum_{n=0}^{\infty}(-a)^{n}\binom{s+n-1}{n}\zeta(s+n),$
(18)

where the binomial coefficient for a complex $s$ has to be understood as

$\binom{s+n-1}{n}=\frac{s(s+1)(s+2)\cdots(s+n-1)}{n!},$ (19)

as well as the analytic continuation of $s\mapsto E(s,\Delta)$. ∎

It is clear that in the limit $\varepsilon\to 0^{+}$ the zeros of
$E(s,1-\varepsilon)$ coincide trivially with the critical ones of the Riemann
zeta function $\zeta(s)$. Let us compute the $\varepsilon\to 0^{+}$
asymptotics of the other nontrivial zeros.

######  Theorem 4.7.

The nontrivial off-critical zeros of $E(s,1-\varepsilon)$ are given, as
$\varepsilon\to 0^{+}$, by $\\{\rho(k)=\rho_{x}(k)+{\rm
i}\rho_{y}(k)\\}_{k\in\mathbb{Z}}$ where

$\displaystyle\rho_{x}(k)=\frac{2}{\ln{2}}\ln\varepsilon+\left(-3+\frac{2}{\ln
2}\ln\pi\right)+\frac{1}{\ln 2}\varepsilon$ $\displaystyle+\frac{1}{\ln
2}\left(\frac{1}{4}+\frac{7\pi^{2}}{24}\right)\varepsilon^{2}+\frac{1}{\ln
2}\left(\frac{1}{12}+\frac{7\pi^{2}}{24}\right)\varepsilon^{3}$
$\displaystyle+\frac{8}{3\ln 2}\left(\frac{\pi^{2}}{8}\right)^{\frac{\ln
3}{\ln 2}}\cos\left[\frac{\ln
3}{\ln{2}}(2k+1)\pi\right]\varepsilon^{2\frac{\ln 3}{\ln
2}}+o\left(\varepsilon^{2\frac{\ln 3}{\ln 2}}\right),$ (20)
$\displaystyle\rho_{y}(k)=\frac{1}{\ln{2}}(2k+1)\pi+\frac{8}{3\ln
2}\left(\frac{\pi^{2}}{8}\right)^{\frac{\ln 3}{\ln 2}}\sin\left[\frac{\ln
3}{\ln{2}}(2k+1)\pi\right]\varepsilon^{2\frac{\ln 3}{\ln 2}}$ (21)
$\displaystyle+o\left(\varepsilon^{2\frac{\ln 3}{\ln 2}}\right).$ (22)

In particular,

  1. 1.

Vanishing of off-critical zeros: we have $\displaystyle\lim_{\varepsilon\to
0^{+}}\rho_{x}(k)=-\infty$;

  2. 2.

Asymptotic crystallization of their imaginary parts: at first order, as
$\varepsilon\to 0^{+}$, the imaginary parts of off-critical zeros are
equidistributed on the lattice $\displaystyle\frac{(2\mathbb{Z}+1)\pi}{\ln
2}$.

###### Proof.

The other nontrivial zeros $\\{\rho\\}$, besides the one of the Riemann zeta
function, correspond to solutions of the equation

$\displaystyle 2^{\rho}$ $\displaystyle=$
$\displaystyle\frac{1-2^{2+\rho}}{2^{5}}\frac{\rho(1+\rho)\zeta(\rho+2)}{\zeta(\rho)}\left(\varepsilon^{2}+\varepsilon^{3}+\frac{3}{4}\varepsilon^{4}+\frac{1}{2}\varepsilon^{5}\right)$
(23)
$\displaystyle+\frac{1}{3}\frac{1-2^{4+\rho}}{2^{11}}\frac{\rho(1+\rho)(2+\rho)(3+\rho)(4+\rho)\zeta(4+\rho)}{\zeta(\rho)}\left(\varepsilon^{4}+2\varepsilon^{5}\right)$
$\displaystyle+O(\varepsilon^{6}).$

As will be showed later, the component $\rho_{x}$ of $\rho=\rho_{x}+{\rm
i}\rho_{y}$ goes to $-\infty$ as $\varepsilon\to 0^{+}$. To simplify our
computations, one applies the well known duality transformation

$\pi^{-s/2}\Gamma\left(\frac{s}{2}\right)\zeta(s)=\pi^{(s-1)/2}\Gamma\left(\frac{1-s}{2}\right)\zeta(1-s)$
(24)

to each Riemann zeta function in (23). Using then the formula
$\Gamma(x+1)=x\Gamma(x)$, one ends up with the result

$\displaystyle 2^{\rho}$ $\displaystyle=$
$\displaystyle-\frac{1-2^{2+\rho}}{2^{3}}\pi^{2}\frac{\zeta(-1-\rho)}{\zeta(1-\rho)}\left(\varepsilon^{2}+\varepsilon^{3}+\frac{3}{4}\varepsilon^{4}+\frac{1}{2}\varepsilon^{5}\right)$
(25)
$\displaystyle+\frac{1}{3}\frac{1-2^{4+\rho}}{2^{7}}\pi^{4}\frac{\zeta(-3-\rho)}{\zeta(1-\rho)}\left(\varepsilon^{4}+2\varepsilon^{5}\right)+O(\varepsilon^{6}).$

In the limit $\varepsilon\to 0^{+}$, the r.h.s. of this equation vanishes and,
consequently, the component $\rho_{x}$ of $\rho=\rho_{x}+{\rm i}\rho_{y}$ must
go to $-\infty$ as indicated before. In the limit $\rho_{x}\to-\infty$, the
ratios of Riemann zeta functions in (25) can be expanded as follows

$\frac{\zeta(-1-\rho)}{\zeta(1-\rho)}=\frac{1+2^{\rho+1}+3^{\rho+1}+\sum_{k\geq
4}k^{\rho+1}}{1+2^{\rho-1}+3^{\rho-1}+\sum_{k\geq
4}k^{\rho-1}}=1+\frac{3}{2}2^{\rho}+\frac{8}{3}3^{\rho}+O(4^{\rho})$ (26)

and

$\frac{\zeta(-3-\rho)}{\zeta(1-\rho)}=\frac{1+2^{\rho+3}+3^{\rho+3}+\sum_{k\geq
4}k^{\rho+3}}{1+2^{\rho-1}+3^{\rho-1}+\sum_{k\geq
4}k^{\rho-1}}=1+\frac{15}{2}2^{\rho}+O(3^{\rho}).$ (27)

In the leading order of the smallness parameter $\varepsilon$, it holds that

$2^{\rho_{x}+{\rm
i}\rho_{x}}=-\frac{\pi^{2}}{8}\varepsilon^{2}+o(\varepsilon^{2}),$ (28)

Since the right-hand side of this equation is real and negative, the leading
order of the $\rho_{y}$-component is given by $2^{{\rm i}\rho_{y}}=-1+o(1)$,
or, equivalently,

$\rho_{y}(k)=\frac{1}{\ln{2}}(2k+1)\pi+o(1),\qquad k\in\mathbb{Z}.$ (29)

This means that in the limit $\varepsilon\to 0^{+}$ there exists an infinite
sequence of equidistant zero components along the $\rho_{y}$ axis. As follows
from (28), the $x$-component of these zeros diverges logarithmically as
$\varepsilon\to 0^{+}$:

$\rho_{x}(k)=\frac{2}{\ln{2}}\ln\varepsilon+\left(-3+\frac{2}{\ln
2}\ln\pi\right)+O(\varepsilon).$ (30)

Note that the leading terms are the same for any value of $k$. This behavior
can be seen in Figure 2.

Higher orders of the expansion of $\rho_{y}(k)$ and $\rho_{x}(k)$ in
$\varepsilon$ can be obtained by inserting the leading order expressions (29)
and (30) directly into the basic relation (25). Performing the expansion
procedure in $\varepsilon$ it is important to realize that

$3^{s}=\left(\frac{\pi^{2}}{8}\varepsilon^{2}\right)^{\frac{\ln 3}{\ln
2}}\exp\left[{\rm i}\frac{\ln 3}{\ln
2}(2k+1)\pi\right]+o\left(\varepsilon^{2\frac{\ln 3}{\ln 2}}\right)$ (31)

is of order $2\ln 3/\ln 2\approx 3.17>3$. After simple algebra one obtains the
desired asymptotics for $\rho_{x}(k)$ and $\rho_{y}(k)$. ∎

Figure 3: Imaginary parts of the first three $(k=0,1,2)$ off-critical zeros
$E(\rho,\Delta)=0$. The deviation $\delta\rho_{y}(k)$ is defined by (32). The
numerical results are depicted by open circles $(k=0)$, squares $(k=1)$ and
triangles $(k=2)$. The plots of $\delta\rho_{y}(k)$ in the region of the small
anisotropy parameter $\varepsilon=1-\Delta$, anticipated to behave according
to the asymptotic formula (33), are represented by the corresponding dashed
curves.

Comparison with our numerics. To check numerically our expansion in
$\varepsilon$ for the imaginary parts of the first three $(k=0,1,2)$ off-
critical zeros, let us define the deviations from their $\varepsilon=0$ values
as follows

$\delta\rho_{y}(k):=\rho_{y}(k)-\frac{1}{\ln{2}}(2k+1)\pi.$ (32)

We know from (22) that the deviations are expected to behave in the region of
the small anisotropy parameter $\varepsilon\to 0^{+}$ as

$\delta\rho_{y}(k)=\frac{8}{3\ln 2}\left(\frac{\pi^{2}}{8}\right)^{\frac{\ln
3}{\ln 2}}\sin\left[\frac{\ln
3}{\ln{2}}(2k+1)\pi\right]\varepsilon^{2\frac{\ln 3}{\ln
2}}+o\left(\varepsilon^{2\frac{\ln 3}{\ln 2}}\right).$ (33)

The numerical results for $\delta\rho_{y}(k)$ are depicted by open circles
$(k=0)$, squares $(k=1)$ and triangles $(k=2)$ in Figure 3. It is seen that
the numerical data fit perfectly the plots given by the asymptotic formula
(33), represented by dashed curves, for small values of $\varepsilon\leq
0.02$.

Figure 4: Real parts of the first three $(k=0,1,2)$ off-critical zeros
$E(\rho,\Delta)=0$. The deviation $\delta\rho_{x}(k)$ is defined by (34). The
numerical results are depicted by open circles $(k=0)$, squares $(k=1)$ and
triangles $(k=2)$. The plots of $\delta\rho_{x}(k)$ in the region of the small
anisotropy parameter $\varepsilon=1-\Delta$, anticipated to behave according
to the asymptotic formula (35), are represented by the corresponding dashed
curves.

As concerns the real parts of the first three $(k=0,1,2)$ off-critical zeros,
we define the deviations as follows

$\displaystyle\delta\rho_{x}(k)$ $\displaystyle:=$
$\displaystyle\rho_{x}(k)-\frac{2}{\ln{2}}\ln\varepsilon-\left(-3+\frac{2}{\ln
2}\ln\pi\right)-\frac{1}{\ln 2}\varepsilon$ (34) $\displaystyle-\frac{1}{\ln
2}\left(\frac{1}{4}+\frac{7\pi^{2}}{24}\right)\varepsilon^{2}-\frac{1}{\ln
2}\left(\frac{1}{12}+\frac{7\pi^{2}}{24}\right)\varepsilon^{3}.$

It is obvious from (22) that the deviations are anticipated to behave for
small values of anisotropy $\varepsilon\to 0^{+}$ as

$\delta\rho_{x}(k)=\frac{8}{3\ln 2}\left(\frac{\pi^{2}}{8}\right)^{\frac{\ln
3}{\ln 2}}\cos\left[\frac{\ln
3}{\ln{2}}(2k+1)\pi\right]\varepsilon^{2\frac{\ln 3}{\ln
2}}+o\left(\varepsilon^{2\frac{\ln 3}{\ln 2}}\right).$ (35)

The numerical results for $\delta\rho_{x}(k)$ are represented by open circles
$(k=0)$, squares $(k=1)$ and triangles $(k=2)$ in Figure 4. The numerical data
fit very well the plots deduced from the asymptotic formula (33) (dashed
curves).

## Acknowledgements

The support received from VEGA Grant No. 2/0092/21 and Project EXSES
APVV-20-0150 is acknowledged.

## References

  * [1] T.M. Apostol, Introduction to analytic number theory, Springer, New York,1976.
  * [2] E.W. Barnes, On the theory of the multiple gamma function, Trans. Camb. Philos. Soc. 19 (1904) 374–-425.
  * [3] J.M. Borwein, M.L. Glasser, R.C. McPhedran, J.G. Wan, J.L. Zucker, Lattice sums then and now, Cambridge University Press, Cambridge, 2013.
  * [4] J.S. Brauchart, Optimal discrete Riesz energy and discrepancy, Unif. Distrib. Theory 6 (2011) 207–220.
  * [5] S. Chowla, A. Selberg, On Epstein’s zeta function, Proc. Natl. Acad. Sci. USA 35 (1949) 371–374.
  * [6] H. Cohn, A. Kumar, Universally optimal distribution of points on spheres, J. Amer. Math. Soc. 20(1) (2007) 99–148.
  * [7] H.M. Edwards, Riemann’s Zeta function, Dover Publications, New York, 1974.
  * [8] E. Elizalde, Ten Physical Applications of Spectral Zeta Functions, Springer Verlag, Berlin, 2012.
  * [9] P. Epstein, Zur Theorie allgemeiner Zetafunctionen, Math. Ann. 56 (1903) 615–644.
  * [10] P. Epstein, Zur Theorie allgemeiner Zetafunctionen II, Math. Ann. 63 (1907) 205–216.
  * [11] N.J. Fine, Note on the Hurwitz zeta-function, Proc. Amer. Math. Soc. 2 (1951) 361-364.
  * [12] J. Hadamard, Étude sur les propriétés des fonction entières et un particulier d’une fonction considéré par Riemann, J. Math. Pure Appl. 9 (1893) 171–215.
  * [13] G.H. Hardy, Sur les zeros de la fonction $\zeta(s)$, Compt. Rend. Acad. Sci. 158 (1914) 1012–-1014.
  * [14] G.H. Hardy, J.E. Littlewood, The zeros of Riemann’s zeta-function on the critical line, Math. Z. 10 (1921) 283–317.
  * [15] A. Hurwitz, Einige Eigenschaften der Dirichletschen Fuctionen $F(s)=\sum\left(\frac{D}{n}\right)\cdot\frac{1}{n^{s}}$ die bei der Bestimmung der Klassenzahlen binärer quadratischer Formen auftreten, Z. Math. Phys. 27 (1882) 86–101.
  * [16] J.I. Hutchinson, On the Roots of the Riemann Zeta-Function. Trans. Amer. Math. Soc. 27 (1925) 49–60.
  * [17] A. Ivić, The Riemann Zeta Function, John Wiley & Sons, New York, 1985.
  * [18] A.M. Jaffe, The Millenium Grand Challenge in Mathematics, Notices of the AMS 53 (2006) 652–660.
  * [19] T. Nakamura, Real zeros of Hurwitz-Lerch zeta and Hurwitz-Lerch type of Euler-Zagierdouble zeta functions, Math. Proc. Cambridge Philos. Soc. 160 (2016) 39–50.
  * [20] B. Riemann, Über die Anzahl der Primzahlen unter einer gegebenen Grösse. Monats-berichte der Berliner Akademie (1859) 671–680.
  * [21] M. Riesz, Sur l’hypothèse de Riemann, Acta Math. 40 (1916) 185–-190.
  * [22] A. Selberg, Contributions to the theory of the Riemann zeta-function, Arch. Math. Naturvid. 48 (1946) 89–-155.
  * [23] R. Spira, Zeros of Hurwitz zeta functions, Mathematics of computations 30 (1976) 863–866.
  * [24] E.C. Titchmarsh, The Zeros of the Riemann Zeta-Function, Proc. Royal Soc. London A 151 (1935) 234–255.
  * [25] E.C. Titchmarsh, The Theory of The Riemann Zeta-function, 2nd ed., Clarendon Press, Oxford, 1988.
  * [26] I. Travěnec, L. Šamaj, Generation of off-critical zeros hypercubic Epstein zeta functions, Appl. Math. Comput. 413 (2022) 126611.
  * [27] L. Vepštas, An efficient algorithm for accelerating the convergence of oscillatory series, useful for computing the polylogarithm and Hurwitz zeta functions, Numer. Algor. 47 (2008) 211–252.
  * [28] W.J. Ventevogel, On the configuration of systems of interacting particles with minimum potential energy per particle, Physica A 92 (1978) 343–361.


11institutetext: Department of Computer Science, University of Copenhagen
22institutetext: Department of Neuroscience, University of Copenhagen
22email<EMAIL_ADDRESS>
# Operating critical machine learning models in resource constrained regimes

Raghavendra Selvan 1122 Julian Schön 11 Erik B Dam 11

###### Abstract

The accelerated development of machine learning methods, primarily deep
learning, are causal to the recent breakthroughs in medical image analysis and
computer aided intervention. The resource consumption of deep learning models
in terms of amount of training data, compute and energy costs are known to be
massive. These large resource costs can be barriers in deploying these models
in clinics, globally. To address this, there are cogent efforts within the
machine learning community to introduce notions of resource efficiency. For
instance, using quantisation to alleviate memory consumption. While most of
these methods are shown to reduce the resource utilisation, they could come at
a cost in performance. In this work, we probe into the trade-off between
resource consumption and performance, specifically, when dealing with models
that are used in critical settings such as in clinics.111Source Code:
https://github.com/raghavian/redl

###### Keywords:

Resource efficiency Image Classification Deep Learning

##  1 Introduction

Every third person on the planet does not yet have universal health coverage
according to the estimates from the United Nations[30]. And improving access
to universal health coverage is one of the UN Sustainable Development Goals.
The use of machine learning (ML) based methods could be instrumental in
achieving this objective. Reliable ML in clinical decision support could
multiply the capabilities of healthcare professionals by reducing the time
spent on several tedious and time-consuming steps. This can free up the
precious time and effort of healthcare professionals to cater to the people in
need[32].

The deployment of ML in critical settings like clinical environments is,
however, currently in a nascent state. In addition to the fundamental
challenges related to multi-site/multi-vendor validation, important issues
pertaining to fairness, bias and ethics of these methods are still being
investigated[24]. The demand for expensive material resources, in terms of
computational infrastructure and energy requirements along with the climate
impact of ML are also points of concern[25].

In this work, we focus on the question of improving the resource efficiency of
developing and deploying deep learning models for resource constrained
environments [28, 3]. Resource constraints in real world scenarios could be
manifested as the need for: low latency (corresponding to urgency), non-
specialized hardware (instead of requiring GPU/similar infrastructure) and
low-energy requirements (if models are deployed on portable settings). The
case for resource efficiency when deploying ML models is quite evident due to
these aforementioned factors. In this work, we also argue that the training
resource efficiency is important, as ML systems will be most useful in the
clinical settings if they are able to continually learn [4, 31]. We validate
the usefulness of several methods for improving resource efficiency of deep
learning models using detailed experiments. We show that resource efficiency –
in terms of reducing compute time, GPU memory usage, energy consumption – can
be achieved without degradation in performance for some classes of methods.

##  2 Methods for Resource Efficiency

In this work we mainly consider reducing the overall memory footprint of deep
learning models, which in turn could reduce computation time and the
corresponding energy costs. We do not explicitly study the influence of model
selection using efficient neural architecture search [29] or model compression
using pruning or similar methods [6], as these techniques require additional
design choices outside of the normal model development and deployment phases.
By focusing on strategies that can be easily incorporated, we also place
emphasis on the possible ease of clinical adoption.

Denote a neural network,
$f_{\mathbf{W}}(\cdot):{\mathbf{X}}\rightarrow{\mathbf{Y}}$ which acts on
inputs ${\mathbf{X}}$ to predict the output ${\mathbf{Y}}$ and has trainable
parameters ${\mathbf{W}}$. For recent deep neural networks, the number of
trainable parameters, $|{\mathbf{W}}|$, can be quite large [26] in
$\mathcal{O}(100M)$.

During training, the weights, ${\mathbf{W}}$, are updated using some form of a
gradient-based update rule. This requires the computation of the gradients of
the loss ${\mathcal{L}}(\cdot)$ with respect to ${\mathbf{W}}$ at iteration
$t$,
${\mathbf{g}}_{t}=\frac{\partial{{\mathcal{L}}}}{\partial{{\mathbf{W}}}}$. In
the case of stateful optimisers that use momentum-based optimisation, such as
Adam [16], the first and second order statistics of the gradient over time,
${\mathbf{m}}_{t},{\mathbf{v}}_{t}$, respectively, are maintained for improved
convergence. Thus, for a neural network with $|{\mathbf{W}}|$ parameters, at
any point during the training there are an additional $\approx
3\cdot|{\mathbf{W}}|$ variables stored in memory. Further, the intermediate
activations at layer $l$, ${\mathbf{h}}_{l}$ are also stored in memory to
efficiently perform backpropagation. Note that all the scalar entries of
${\mathbf{W}},{\mathbf{g}}_{t},{\mathbf{m}}_{t},{\mathbf{v}}_{t},{\mathbf{h}}_{l}\in{\mathbb{R}}$.
On most computers, these real numbers are discretised into floating point-32
(FP32) format; wherein each variable requires 32 bits.

Resource efficiency in this work is primarily addressed by reducing the
precision of these variables by quantisation [13, 28]. The implication of
using quantised or low-precision variables is that it not only reduces memory
footprint but can also reduce the compute cost of matrix multiplications [15]
at the expense of loss in precision. Note, however, that the overhead of
performing quantisation in some cases might outweigh the gains in computation.

In this work, we investigate a combination of the following three quantisation
strategies:

  1. 1.

Gradients and intermediate activations: Drastic quantisation of
${\mathbf{h}}_{l},{\mathbf{g}}_{t}$ have been studied extensively in
literature [15, 5]. For instance, in [15] the gradients and activations were
quantised to 1-bit and still yielded reasonable performance compared to the
32-bit models. However, not all operations are well suited for quantisation.
In frameworks like Pytorch [22], the precision of
${\mathbf{h}}_{l},{\mathbf{g}}_{t}$ are automatically cast into 16-bit
precision, and recast to 32-bit when necessary. This technique known as
automatic mixed precision (AMP) can yield considerable reduction of memory
costs [19].

  2. 2.

Optimiser states: Until recently, the main focus of quantisation in deep
learning has been on compressing
${\mathbf{g}}_{t},{\mathbf{h}}_{l},{\mathbf{W}}$. In a recent work in [10], it
was shown that the optimiser states, ${\mathbf{m}}_{t},{\mathbf{v}}_{t}$ can
use up to 75% of memory. By quantising optimiser states to 8-bits using a
dynamic quantisation scheme, this work reduces the memory consumption during
model training.

  3. 3.

Model weights: The resolution of the hypothesis space is controlled by the
precision of the weights, ${\mathbf{W}}$. Using lower precision weights for
models can reduce their expressive power. This is most evident when a model
trained in full precision is then cast into lower precision. In such cases, it
is a common strategy to fine tune the quantised model for a few epochs to
recover the loss in performance [20]. However, when a model is trained from
the outset in low-precision modes, they can achieve similar performances to
models with full precision weights [27].

##  3 Data & Experiments

Nation-wide screening for early detection of cancer is being attempted in many
countries. For instance, early breast cancer detection using mammography
screening has shown positive impact, both for the subjects and in allocating
healthcare resources [8]. To simulate a scenario where ML models are deployed
in such screening programs, we use a publicly available mammography screening
dataset for the task of recommending follow-up scans. We study the influence
of using resource optimised models and their impact on the predictions in the
mammography screening.
Datasets: The experiments in this work were primarily performed using a subset
of the data available as part of the Radiological Society of North America
(RSNA) Mammography Breast Cancer Detection
challenge[7]222https://www.kaggle.com/competitions/rsna-breast-cancer-
detection/overview. Accessed on 08/03/2023. Subjects under the age of $55$ yr
in the training set were used as the cohort of interest for early detection.
This yielded a dataset comprising $3199$ unique subjects. As each subject has
multiple scans (also from different views), the final dataset consists of
$11856$ mammography scans. While the dataset provides multiple markers which
could be used as labels, we focus on the BI-RADS333BI-RADS:Breast Imaging
Reporting & Data System score[21]. The dataset provides three grades of the
BI-RADS score: 0 if the breast required follow-up, 1 if the breast was rated
as negative for cancer, and 2 if the breast was rated as normal. We combine
the classes 1 and 2 to make a binary classification task of predicting if a
follow-up scan was required or not. This resulted in the following label
distribution: [$7903$, $3953$], corresponding to requiring follow-up (class-1)
and no follow-up (class-0). The dataset was randomly split into training,
validation and test sets based on subject identities. This resulted in a split
of [1919,640,640] subjects in the three sets, respectively. In terms of
mammography images, the resulting split was [7113,2371,2372] for the training,
validation and test sets. Sample images from the dataset and the age
distribution of the subjects is shown in Figure 1. We acknowledge that this
classification task is a simplification of the mammography reading done during
screening, but depending on the screening setup (e.g. number of readers), even
this simple task could supplement the workflow.

(A)

(B)

Figure 1: (A): Sample mammography images from the RSNA training set showing
the diverse acquisition and anatomical variations. (B): Age distribution of
the subjects in the dataset used. Maximum age of $55$ yr was used.

Additional experiments on a second dataset for lung nodule detection are also
reported for further validation. We use the LIDC-IDRI dataset which consists
of 1018 thoracic CT scans and predict the presence or absence of tumours [2].
Further details of the LIDC classification dataset are provided in the
Supplementary Material.
Classification models: To evaluate the impact of resource constraints,
different classes of deep learning models are used. Specifically, we use
Densenet [14], Efficientnet [29], Vision transformer [11] and Swin Transformer
[18]. All models were pretrained on ImageNet and fine tuned on the two
datasets444Pretrained weights for the models were obtained from
https://timm.fast.ai/.
Experimental design: All the models were trained for a maximum of 50 epochs,
with an early stopping patience of 5 epochs using batch size of 32. All the
experiments were performed on a desktop workstation with Intel i7 processor,
32GB memory and Nvidia RTX 3090 GPU with 24 GB memory. All experiments were
performed in Pytorch (v1.7.1) [22], and unless otherwise specified using the
Adam [16] optimizer with learning rate of $10^{-5}$. Measurements of run time,
energy and carbon consumption were performed using CarbonTracker [1]. All
models were trained three times with different random seeds. The training of
models in this work used 25.0 kWh of electricity contributing to 4.7 kg of
CO2eq.
Experiments: Resource constraints in real world scenarios could be due to:
compute time (corresponding to low latency), specialized hardware (access to
GPU infrastructure) and energy consumption (if models are deployed on portable
devices). Using these and the test set performance as the axes, we run all the
classification models with- and without- any resource optimisation. The three
main optimisation techniques we use are related to the methods described in
Section 2. Firstly, we perform automatic mixed precision training of the
models in order to obtain low compute and memory models [19]. This is
implemented using the automatic mixed precision (AMP) feature in Pytorch.
Secondly, we minimize the resource utilisation of the optimizer using an 8-bit
optimiser based on the implementation5558-bit optimiser from bitsandbytes from
[10]. Finally, we also investigate training of a low precision model by
casting all the model weights to half precision (float-16). We then compare
the performance of all the classification models without any explicit resource
optimisation and with the combination of the mixed precision training/half-
precision models and the 8-bit optimiser. Note that AMP cannot be combined
with Half Precision models, as the model itself is in low precision. We use
binary accuracy on the test set ($P_{T}$), average training time to
convergence ($T_{c}$), maximum GPU memory required, average inference time on
test set ($T$) and the total energy consumption ($E$) as the measures to
compare the models and different resource efficiency techniques.

Table 1: Quantitative comparison of the different classification methods reported over three random initialisations for RSNA and LIDC datasets. The use of 8-bit optimizer (8bit), automatic mixed precision (AMP), and half precision model (Half) are marked. Number of parameters: $|{\mathbf{W}}|$, mean accuracy with standard deviation over three seeds, average GPU memory required (in GB), average inference time, and average energy consumption of different settings reported when predicting on test set. Models that diverged at least during one of the random initialisations are marked with ∗ in $P_{T}$ column. Dataset | Method |  $|{\mathbf{W}}|$(M) | 8bit | AMP | Half | $\mathbf{P}_{T}$ | GPU |  $\mathbf{T}_{c}$(s) |  T(s) |  E(Wh)
---|---|---|---|---|---|---|---|---|---|---
RSNA[7] | Densenet[14] | 6.9 | ✗ | ✗ | ✗ |  0.738$\pm$0.005 | 6.6 | 361.3 | 2.7 | 104.2
✗ | ✓ | ✗ |  0.740$\pm$0.005 | 3.3 | 296.0 | 2.7 | 59.2
✓ | ✗ | ✗ |  0.743$\pm$0.003 | 6.5 | 240.0 | 2.7 | 58.1
✓ | ✗ | ✓ | 0.745$\pm$0.004 | 3.2 | 170.3 | 2.1 | 34.6
✓ | ✓ | ✗ | 0.744$\pm$0.006 | 3.2 | 191.3 | 2.7 | 42.3
Swin Trans.[18] | 86.9 | ✗ | ✗ | ✗ |  0.739$\pm$0.008 | 15.3 | 617.6 | 7.5 | 100.1
✗ | ✓ | ✗ |  0.735$\pm$0.006 | 10.7 | 497.0 | 7.5 | 109.8
✓ | ✗ | ✗ |  0.739$\pm$0.012 | 14.3 | 732.3 | 7.5 | 86.1
✓ | ✗ | ✓ |  0.733$\pm$0.013 | 7.5 | 316.0 | 4.3 | 118.5
✓ | ✓ | ✗ |  0.712$\pm$0.01∗ | 9.7 | 588.6 | 7.5 | 95.8
ViT[11] | 116.7 | ✗ | ✗ | ✗ |  0.728$\pm$0.001 | 5.6 | 322.3 | 1.9 | 52.0
✗ | ✓ | ✗ |  0.728$\pm$0.006 | 5.6 | 249.0 | 1.9 | 39.8
✓ | ✗ | ✗ |  0.716$\pm$0.009 | 4.1 | 497.0 | 1.9 | 81.0
✓ | ✗ | ✓ |  0.690$\pm$0.02∗ | 3.5 | 69.6 | 1.0 | 10.7
✓ | ✓ | ✗ |  0.665$\pm$0.02∗ | 4.7 | 132.6 | 1.9 | 21.0
LIDC[2] | Densenet[14] | 6.9 | ✗ | ✗ | ✗ |  0.656$\pm$0.005 | 6.6 | 674.3 | 3.4 | 104.2
✗ | ✓ | ✗ |  0.655$\pm$0.010 | 3.3 | 425.7 | 3.4 | 59.2
✓ | ✗ | ✗ |  0.677$\pm$0.008 | 6.5 | 370.7 | 3.4 | 58.1
✓ | ✗ | ✓ |  0.675$\pm$0.007 | 3.2 | 234.3 | 2.7 | 34.6
✓ | ✓ | ✗ |  0.676$\pm$0.005 | 3.2 | 291.3 | 2.7 | 42.3
Swin Trans.[18] | 86.9 | ✗ | ✗ | ✗ |  0.697$\pm$0.018 | 15.3 | 678.7 | 9.6 | 109.8
✗ | ✓ | ✗ |  0.696$\pm$0.018 | 10.8 | 544.0 | 9.6 | 86.1
✓ | ✗ | ✗ |  0.684$\pm$0.016 | 14.4 | 727.3 | 9.6 | 118.5
✓ | ✗ | ✓ | 0.704$\pm$0.007 | 9.9 | 401.0 | 5.6 | 64.1
✓ | ✓ | ✗ |  0.684$\pm$0.020 | 7.5 | 615.3 | 9.6 | 99.9
ViT[11] | 116.7 | ✗ | ✗ | ✗ |  0.610$\pm$0.041 | 5.6 | 231.7 | 2.4 | 37.2
✗ | ✓ | ✗ |  0.608$\pm$0.043 | 5.5 | 166.3 | 2.4 | 25.8
✓ | ✗ | ✗ |  0.636$\pm$0.012 | 4.3 | 415.3 | 2.4 | 67.1
✓ | ✗ | ✓ |  0.564$\pm$0.072 | 3.4 | 123.0 | 1.3 | 15.2
✓ | ✓ | ✗ |  0.561$\pm$0.076 | 4.5 | 99.0 | 2.4 | 19.1

Results: Experiments with the various resource optimisation settings and
methods for both datasets are reported comprehensively in Table 1. For each
method, five settings are reported. We treat the model with no additional
resource optimisation as the baseline, reported in the first row for each
method. For each model, we then report the following settings: [AMP, 8bit
optimizer, AMP+ 8bit optimizer, Half precision + 8bit optimiser] yielding five
settings for each model. The models are sorted based on the total number of
trainable parameters, reported in the third column.

At the outset, it is worth noting that the baseline model performance on the
test set of all Densenet and the transformer-based Swin- and Vision-
Transformers are comparable ($\approx 0.74$) in the RSNA mammography dataset.

The use of AMP during training does not degrade the performance of any of the
models, while reducing the maximum GPU memory utilised (up to 50% in most
cases). The run time, $T_{c}$ for all the models is also decreased when using
AMP.

The use of 8-bit optimiser reduces the GPU memory utilised and the convergence
time. The more interesting observation is that in almost all cases (except
ViT), it also converges to a better solution, yielding a small performance
improvement. This is consistent with the observations reported in [10], where
this behaviour is attributed to the elimination of outliers (which otherwise
could lead to exploding gradients) during quantisation and dequantisation of
the optimiser states.

Setting the model weights to half precision (16-bit instead of 32-bit)
immediately gives a reduction in the GPU memory utilisation. In our
experiments, however, there were stability issues when the half precision
model was trained using the standard Adam optimiser. Using the 8-bit optimiser
which is quantisation-aware helped convergence [10].

The best performing configuration across all the models for RSNA dataset is
the Densenet with 8-bit optimiser either when used with AMP or the half
precision model. This is an interesting result showing that these techniques
not only reduce resource utilisation but could also improve model convergence.

The observations about the influence of AMP, 8-bit optimiser and half-
precision hold within each model category for the LIDC dataset. The best
performing configuration across all the models for LIDC dataset is the Swin
Transformer with 8-bit optimiser when used with the half precision model.

The performance trends in Table 1 are captured visually in the radar plots
presented in Figure 2-A for Densenet. The mean performance measures of
Densenet over the three random initialisations, for the five settings, are
shown in different colours. The four axes are used to report four key
measures:[$P_{T},E,\text{GPU},T$]. The optimal model would have larger
footprint along $P_{T}$ and smaller values on the other axes. We notice that
the baseline model (in blue) spans a larger region over all axes, wheres the
half precision model trained with 8-bit optimiser (in red) has smaller span.
These trends are also visualised for the Vision transformer model in Figure
2-B, and for other models in the Supplementary Material.

(A): Densenet

(B): ViT

(C)

Figure 2: (A & B): Radar plots showing the mean metrics for performance
($P_{T},E,\text{GPU},T$) reported in Table 1 for the five settings, shown for
Densenet [14] and Vision transformer [11]. (C): Influence of different scales
of downsampling on the test performance for three configurations of Densenet
models on the RSNA dataset. For each, the mean and the standard deviation of
the three runs is shown as the curve and the shaded area, respectively.

##  4 Discussion

Data Quality: One common method to manage the resource constraints during
training or inference of deep learning models is the use of reduced quality
data. For instance, in image analysis, downsampling of images is a commonly
used strategy as this can have a direct impact on the size of the models. Such
reduction in data quality are undoubtedly lossy. We illustrate this for the
mammography screening task for the Densenet models in Figure 2-C. We reduce
the quality of the images both during training and inference by downsampling
the image at different scales: $S=[2,4,8]$. We observe a clear performance
degradation compared to the models trained with full resolution ($S=1$). While
this might be unsurprising, we argue, in lieu of the results in Table 1 that
optimising models instead of degrading data quality can offer a better trade-
off between performance and resource efficiency.
Resource Efficiency and Model Complexity: The most common measure of model
complexity is the number of parameters. The experiments and results reported
in Table 1 reveal slightly surprising trends. The best performing model for
the RSNA dataset (Densenet) had the fewest number of parameters (6.9M), and
ViT which was the largest model (116.7M) was also the fastest at inference
time ($<$2s). This is due to the fact that during training the resource
overhead for deep learning models is due to factors other than model weights,
such as the optimiser states or other derived data like gradients/activation
maps as discussed in Section 2. We have demonstrated that resource efficiency
can also be achieved by optimising these other axes.
Resource Efficiency of Vision Transformers: The use of AMP or 8-bit optimiser
settings degrades the performance of both the transformer-based models to some
extent in both datasets. This is in contrast to the CNN models, where there is
no large performance degradation. One reason for the degradation in
performance in transformers when using the 8-bit optimiser could be
instability during training. The authors in [10] suggest the use of a stable
embedding layer when training transformers for NLP tasks. For a fairer
comparison with the CNN models, we do not use the stable embedding layer for
the two transformer models. Our experiments indicate that CNN models could be
less sensitive when compared to transformer-based models in low/mixed
precision or quantised settings.
QALY-like metrics: Quality Adjusted Life Year (QALY) is a utilitarian index
used to estimate the cost-effectiveness of a health treatment; it attempts to
arithmetically relate quality and quantity of life [23]. In our experiments,
we did not notice degradation in performance for some of the model classes.
However, when there is a drop in performance, the tolerance levels in clinical
settings would require formulating a Pareto optimisation of performance and
resources [9]. It is in this context that a metric similar to QALY could be
useful, when balancing resource constraints with critical performance
guarantees.
Limitations: This work focused on optimising resource efficiency using only a
few of the techniques (mixed precision training, half precision models, 8-bit
optimisers). There are several other methods which could further optimise the
resource consumption. For instance, formulating model selection that can yield
inherently efficient models within neural architecture search (NAS) framework
could be a possibility. However, performing NAS in itself is prohibitively
expensive [12]. The other main limitation of the work is the gap in validating
these methods in actual, resource constrained setting which the authors did
not have access to. Testing these resource optimisation measures in diverse
settings such as on different edge devices, for instance, could be an
interesting direction for future work.

##  5 Conclusion

Large deep learning models are advancing medical image analysis and will
undoubtedly accelerate clinical workflows. For these models to make the most
impact, globally, resource efficiency is an important factor. In this work, we
have evaluated a subset of resource optimisation methods that can reduce
memory consumption, latency, and energy consumption, without any performance
degradation (Densenet in RSNA, Swin Trans. in LIDC). In contrast to our
initial hypothesis, surprisingly, there is no performance degradation across
the board; only systematic ones (like with Transformers which can be handled).
This prompts us to recommend these resource optimisation methods to be part of
the standard procedures when developing and deploying deep learning methods.

Acknowledgments The authors would like to thank Christian Igel (UCPH) and Lee
Lassen (Ambu A/S) for fruitful discussions, and Pedram Bakhtiarifard (UCPH)
for their public repository for making radar plots. RS and JS are partly
funded by European Union’s HorizonEurope research and innovation programme
under grant agreements No. 101070284 and No. 101070408.

## References

  * [1] Anthony, L.F.W., Kanding, B., Selvan, R.: Carbontracker: Tracking and predicting the carbon footprint of training deep learning models. ICML Workshop on Challenges in Deploying and monitoring Machine Learning Systems (July 2020), arXiv:2007.03051
  * [2] Armato III, S.G., McLennan, G., McNitt-Gray, M.F., Meyer, C.R., Yankelevitz, D., Aberle, D.R., Henschke, C.I., Hoffman, E.A., Kazerooni, E.A., MacMahon, H., et al.: Lung image database consortium: developing a resource for the medical imaging research community. Radiology 232(3), 739–748 (2004)
  * [3] Bartoldson, B.R., Kailkhura, B., Blalock, D.: Compute-efficient deep learning: Algorithmic trends and opportunities. arXiv preprint arXiv:2210.06640 (2022)
  * [4] Baweja, C., Glocker, B., Kamnitsas, K.: Towards continual learning in medical imaging. Workshop on Medical Imaging meets NeurIPS (2018), arXiv:1811.02496
  * [5] Chakrabarti, A., Moseley, B.: Backprop with approximate activations for memory-efficient network training. Advances in Neural Information Processing Systems 32 (2019)
  * [6] Cheng, Y., Wang, D., Zhou, P., Zhang, T.: Model compression and acceleration for deep neural networks: The principles, progress, and challenges. IEEE Signal Processing Magazine 35(1), 126–136 (2018)
  * [7] Chris Carr, FelipeKitamura, H.k.i.J.K.C.J.M.K.A.M.V.M.R.R.B.S.D.: Rsna screening mammography breast cancer detection (2022), https://kaggle.com/competitions/rsna-breast-cancer-detection
  * [8] Daysal, N.M., Mullainathan, S., Obermeyer, Z., Sarkar, S.K., Trandafir, M.: An economic approach to machine learning in health policy. Univ. of Copenhagen Dept. of Economics Discussion, CEBI Working Paper (24) (2022)
  * [9] Deb, K.: Multi-objective optimisation using evolutionary algorithms: an introduction. Springer (2011)
  * [10] Dettmers, T., Lewis, M., Shleifer, S., Zettlemoyer, L.: 8-bit optimizers via block-wise quantization. In: International Conference on Learning Representations (2022), https://openreview.net/forum?id=shpkpVXzo3h
  * [11] Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.: An image is worth 16x16 words: Transformers for image recognition at scale. In: International Conference on Learning Representations (2020)
  * [12] Elsken, T., Metzen, J.H., Hutter, F.: Neural architecture search: A survey. The Journal of Machine Learning Research 20(1), 1997–2017 (2019)
  * [13] Fiesler, E., Choudry, A., Caulfield, H.J.: Weight discretization paradigm for optical neural networks. In: Optical interconnections and networks. vol. 1281, pp. 164–173. SPIE (1990)
  * [14] Huang, G., Liu, Z., Van Der Maaten, L., Weinberger, K.Q.: Densely connected convolutional networks. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 4700–4708 (2017)
  * [15] Hubara, I., Courbariaux, M., Soudry, D., El-Yaniv, R., Bengio, Y.: Quantized neural networks: Training neural networks with low precision weights and activations. The Journal of Machine Learning Research 18(1), 6869–6898 (2017)
  * [16] Kingma, D., Ba, J.: Adam optimizer. arXiv preprint arXiv:1412.6980 pp. 1–15 (2014)
  * [17] Knegt, S.: A Probabilistic U-Net for segmentation of ambiguous images implemented in PyTorch. https://github.com/stefanknegt/Probabilistic-Unet-Pytorch (2018)
  * [18] Liu, Z., Lin, Y., Cao, Y., Hu, H., Wei, Y., Zhang, Z., Lin, S., Guo, B.: Swin transformer: Hierarchical vision transformer using shifted windows. In: Proceedings of the IEEE/CVF international conference on computer vision. pp. 10012–10022 (2021)
  * [19] Micikevicius, P., Narang, S., Alben, J., Diamos, G., Elsen, E., Garcia, D., Ginsburg, B., Houston, M., Kuchaiev, O., Venkatesh, G., Wu, H.: Mixed precision training. In: International Conference on Learning Representations (2018), https://openreview.net/forum?id=r1gs9JgRZ
  * [20] Nagel, M., Fournarakis, M., Amjad, R.A., Bondarenko, Y., Van Baalen, M., Blankevoort, T.: A white paper on neural network quantization. arXiv preprint arXiv:2106.08295 (2021)
  * [21] Orel, S.G., Kay, N., Reynolds, C., Sullivan, D.C.: Bi-rads categorization as a predictor of malignancy. Radiology 211(3), 845–850 (1999)
  * [22] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al.: Pytorch: An imperative style, high-performance deep learning library. In: Advances in Neural Information Processing Systems. pp. 8024–8035 (2019)
  * [23] Prieto, L., Sacristán, J.A.: Problems and solutions in calculating quality-adjusted life years (qalys). Health and quality of life outcomes 1, 1–8 (2003)
  * [24] Ricci Lara, M.A., Echeveste, R., Ferrante, E.: Addressing fairness in artificial intelligence for medical imaging. nature communications 13(1), 4581 (2022)
  * [25] Selvan, R., Bhagwat, N., Wolff Anthony, L.F., Kanding, B., Dam, E.B.: Carbon footprint of selecting and training deep learning models for medical image analysis. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2022: 25th International Conference, Singapore, September 18–22, 2022, Proceedings, Part V. pp. 506–516. Springer (2022)
  * [26] Sevilla, J., Heim, L., Ho, A., Besiroglu, T., Hobbhahn, M., Villalobos, P.: Compute trends across three eras of machine learning. arXiv preprint arXiv:2202.05924 (2022)
  * [27] Sun, X., Wang, N., Chen, C.Y., Ni, J., Agrawal, A., Cui, X., Venkataramani, S., El Maghraoui, K., Srinivasan, V.V., Gopalakrishnan, K.: Ultra-low precision 4-bit training of deep neural networks. Advances in Neural Information Processing Systems 33, 1796–1807 (2020)
  * [28] Sze, V., Chen, Y.H., Yang, T.J., Emer, J.S.: Efficient processing of deep neural networks: A tutorial and survey. Proceedings of the IEEE 105(12) (2017)
  * [29] Tan, M., Le, Q.: Efficientnet: Rethinking model scaling for convolutional neural networks. In: International conference on machine learning. pp. 6105–6114 (2019)
  * [30] United-Nations: Sustainable development goals report (2022)
  * [31] Vokinger, K.N., Feuerriegel, S., Kesselheim, A.S.: Continual learning in medical devices: Fda’s action plan and beyond. The Lancet Digital Health 3(6), e337–e338 (2021)
  * [32] Yu, K.H., Beam, A.L., Kohane, I.S.: Artificial intelligence in healthcare. Nature biomedical engineering 2(10), 719–731 (2018)

##  Appendix 0.A Additional Results

Figure 3: LIDC-IDRI dataset comprises 1018 thoracic CT images with lesions annotated by four radiologists [2]. Patches of 128x128 px are extracted from the 2D slices, yielding a total of 15096 patches [17]. Each patch has annotations from four raters marking the tumour regions. These segmentation masks were converted into binary labels indicating the presence (if $\geq 2$ radiologists marked a tumour) or absence of tumours (if $<2$ raters marked tumours), resulting in a fairly balanced dataset. All image intensities are normalised to be in [0, 1]. Training, validation and test splits are made following a [0.6,0.2,0.2] ratio. Table 2: Quantitative comparison of MLP and EfficientNet [28] reported over three random initialisations for the RSNA dataset. The use of 8-bit optimizer (8b_Opt), automatic mixed precision (AMP) and half precision model (Half) are marked. Number of parameters: $|\theta|$, mean accuracy with standard deviation over three seeds, average GPU memory required (in GB), average inference time and average energy consumption of different settings reported when predicting on test set. Method |  $|\theta|$(M) | 8b_Opt | AMP | Half | Acc. | GPU |  Tc(s) |  T(s) |  E (Wh)
---|---|---|---|---|---|---|---|---|---
MLP | 8.4 | ✗ | ✗ | ✗ |  0.685$\pm$0.010 | 0.4 | 85.0 | 0.2 | 7.3
✗ | ✓ | ✗ |  0.687$\pm$0.005 | 0.4 | 92.3 | 0.2 | 8.1
✓ | ✗ | ✗ |  0.693$\pm$0.007 | 0.2 | 41.3 | 0.2 | 3.4
✓ | ✓ | ✗ | 0.694$\pm$0.005 | 0.2 | 55.3 | 0.2 | 4.7
✓ | ✗ | ✓ |  0.680$\pm$0.005 | 0.1 | 27.3 | 0.2 | 2.4
Efficientnet[29] | 51.1 | ✗ | ✗ | ✗ |  0.679$\pm$0.026 | 14.6 | 666.7 | 4.9 | 104.3
✗ | ✓ | ✗ |  0.692$\pm$0.005 | 7.8 | 4027.3 | 4.9 | 618.4
✓ | ✗ | ✗ |  0.715$\pm$0.022 | 14.2 | 466.0 | 4.9 | 74.4
✓ | ✓ | ✗ | 0.717$\pm$0.008 | 7.4 | 2769.0 | 4.9 | 427.2
✓ | ✗ | ✓ | 0.717$\pm$0.008 | 7.4 | 2769.0 | 3.8 | 427.2

(A): MLP

(B): EfficientNet [28]

Figure 4: (A & B): Radar plots for multi-layered perceptron (MLP) and
EfficientNet [28] showing the mean metrics for performance
($P_{T},E,\text{GPU},T$) reported in Table 2 for the five settings.

(A)

(B)

(C)

Figure 5: Radar plots for DenseNet [14], Swin Transformer [18] and Vision
Transformer [11] on the LIDC dataset, showing the mean metrics for performance
($P_{T},E,\text{GPU},T$) reported in Table 1 for the five settings.


# Upper limits on the Polarized Isotropic Stochastic Gravitational-Wave
Background from Advanced LIGO-Virgo’s First Three Observing Runs

Yang Jiang  CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China  Qing-Guo Huang  Corresponding author:
<EMAIL_ADDRESS>CAS Key Laboratory of Theoretical Physics, Institute of
Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences, No.
19A Yuquan Road, Beijing 100049, China  School of Fundamental Physics and
Mathematical Sciences Hangzhou Institute for Advanced Study, UCAS, Hangzhou
310024, China

###### Abstract

Parity violation is expected to generate an asymmetry between the amplitude of
left and right-handed gravitational-wave modes which leads to a circularly
polarized stochastic gravitational-wave background (SGWB). Due to the three
independent baselines in the LIGO-Virgo network, we focus on the amplitude
difference in strain power characterized by Stokes’ parameters and do maximum-
likelihood estimation to constrain the polarization degree of SGWB. Our
results indicate that there is no evidence for the circularly polarized SGWB
in the data. Furthermore, by modeling the SGWB as a power-law spectrum, we
place upper limit on the normalized energy density
$\Omega_{\text{gw}}(25\,\text{Hz})<5.3\times 10^{-9}$ at $95\%$ confidence
level after marginalizing over the polarization degree and spectral index.

Introduction.  The stochastic gravitational-wave background (SGWB) is a
superposition of gravitational waves (GWs) from numerous unresolved and
uncorrelated sources. These GW sources can arise either from astrophysical
processes like compact binary coalescences (CBCs) Rosado (2011); Wu  _et al._
(2012); Zhu  _et al._ (2013), rotating neutron stars Ferrari  _et al._ (1999);
Regimbau and de Freitas Pacheco (2001); Lasky  _et al._ (2013), stellar core
collapses Crocker  _et al._ (2015), or cosmological contributions, such as
phase transitions (PT) Witten (1984); Kosowsky  _et al._ (1992); Dev and
Mazumdar (2016); Von Harling  _et al._ (2020), cosmic strings Damour and
Vilenkin (2000); Sarangi and Tye (2002); Siemens  _et al._ (2007); Abbott  _et
al._ (2018) and inflation models Turner (1997); Guzzetti  _et al._ (2016). The
detection of SGWB will provide a better understanding about the distribution
of astrophysical sources, the history of early Universe and testing the
theories of gravity. Until now, the terrestrial laser interferometers like
Advanced LIGO Aasi  _et al._ (2015) and Virgo Acernese  _et al._ (2015) have
been combined to search for SGWB in their output strain data. Based on current
noise power, the results show that there is no detectable correlation and thus
upper limits on the energy density of both isotropic and anisotropic SGWB have
been set in Abbott  _et al._ (2021a, b, 2022).

The SGWB is widely assumed to be unpolarized, implying that the different
polarized GW modes are treated as statistically identical and independent in
the analyses. This assumption is quite reliable for a superposition of signals
from clustered sources and parity is conserved in general relativity (GR).
However, once the parity is violated, the asymmetry between the amplitude of
left and right-handed GW modes leads to a circularly polarized SGWB. Some
processes in the early Universe, such as helical turbulence during a first-
order PT Kahniashvili  _et al._ (2005), can generate a polarized SGWB. When we
go beyond GR, various theories of gravity, including Chern-Simons gravity
Alexander and Yunes (2009); Lyth  _et al._ (2005); Satoh  _et al._ (2008);
Saito  _et al._ (2007), Hořava-Lifshitz gravity Horava (2009); Takahashi and
Soda (2009) and ghost-free scalar-tensor gravity Crisostomi  _et al._ (2018)
etc, can also give rise to parity violations. Accounting for these theories,
particular polarization mode can be enhanced or compressed through amplitude
birefringence effect Zhao  _et al._ (2020) during its propagation, resulting
in unequal left and right-handed components in SGWB. In all, detecting a
circularly polarized SGWB might yield a profound discovery for fundamental
physics.

In this letter, we perform the first multi-baseline search for the circularly
polarized isotropic SGWB in the data of Advanced LIGO–Virgo’s First Three
Observing Runs. The method for detecting the circularly polarized SGWB by
ground-based interferometers was proposed for the first time in Seto and
Taruya (2007, 2008). This method depends on the detector’s nontrivial reaction
to the Stokes V parameter. Even though SGWB can be detected in a one-baseline
analysis, the existence of polarization cannot be affirmed because one
baseline cannot distinguish the left-handed component from the right-handed
one in principle. In this sense, in Crowder  _et al._ (2013), the authors can
obtain the upper limit on parity violation indirectly by assuming a fiducial
model with a power-law spectrum of SGWB and then detecting the deviation,
because only the correlation of LIGO Hanford-LIGO Livingston pair was
available at that time. In principle, there are two different polarization
modes, at least one more baseline is needed for detecting the circularly
polarized SGWB. Fortunately, from O3 observing run, Virgo has been involved in
searching for SGWB. In this letter, the extra outputs of LIGO-Virgo pairs are
firstly adopted for distinguishing different polarization modes in the data.

Circularly Polarized SGWB.  The SGWB is formed as a superposition of plane
waves propagating along all possible directions with various frequencies:

$h_{ij}(t,\bm{x})=\sum_{A}\int_{-\infty}^{\infty}\mathrm{d}f\int\mathrm{d}^{2}\Omega\;h_{A}(f,\bm{\Omega})\mathrm{e}^{2\pi\mathrm{i}f(t-\bm{\Omega}\cdot\bm{x}/c)}\bm{e}_{ij}^{A}(\bm{\Omega}),$
(1)

where $\bm{e}_{ij}^{A}$ is the polarization basis tensor. As far as we are
concerned, the circularly polarization basis $A=\\{R,L\\}$ is a favorable
choice which is related to the linearly polarization basis by

$\bm{e}_{ij}^{R}=\frac{1}{\sqrt{2}}\left(\bm{e}^{+}_{ij}+\mathrm{i}\bm{e}^{\times}_{ij}\right),\;\bm{e}_{ij}^{L}=\frac{1}{\sqrt{2}}\left(\bm{e}^{+}_{ij}-\mathrm{i}\bm{e}^{\times}_{ij}\right).$
(2)

For isotropic SGWB, the quadratic expectation values of $h_{A}$ is

$\langle
h_{A}(f,\bm{\Omega})h^{*}_{A^{\prime}}(f^{\prime},\bm{\Omega}^{\prime})\rangle=\frac{1}{2}S_{h}^{A}(f)\delta_{AA^{\prime}}\delta(f-f^{\prime})\frac{\delta^{2}(\bm{\Omega},\bm{\Omega}^{\prime})}{4\pi},$
(3)

where

$S_{h}^{R}=\frac{1}{2}\left(I(f)+V(f)\right),\quad
S_{h}^{L}=\frac{1}{2}\left(I(f)-V(f)\right),$ (4)

and $I(f),V(f)$ are the so-called Stokes’ parameters. To reveal their
meanings, for a sinusoidal plane wave, the concise expressions for both of
them are given by

$I=|h_{R}|^{2}+|h_{L}|^{2},\quad V=|h_{R}|^{2}-|h_{L}|^{2}.$ (5)

In this sense, $I(f)$ and $V(f)$ denote the total intensity of SGWB and the
degree of parity violation, respectively. If $V\neq 0$, the parity is
violated. The normalized energy density of SGWB is defined by

$\Omega_{\text{gw}}(f)=\frac{1}{\rho}_{c}\frac{\mathrm{d}\rho_{\text{gw}}}{\mathrm{d}\ln
f},$ (6)

where $\rho_{c}=3H_{0}^{2}c^{2}/(8\pi G)$ is the critical energy density and
$\Omega_{\text{gw}}$ is related to the Stokes’ parameter $I(f)$ by

$\Omega_{\text{gw}}(f)=\frac{2\pi^{2}f^{3}}{3H_{0}^{2}}I(f).$ (7)

Separating different circularly polarized modes.  The detection of SGWB relies
on the correlation analysis Allen and Romano (1999) of strain data. If two
detectors couple together to form a baseline, during a certain period $T$, a
cross-correlating statistic is defined as

$\hat{C}(f)=\frac{2}{T}\mathrm{Re}\left[\tilde{s}_{1}^{*}(f)\tilde{s}_{2}(f)\right],$
(8)

where $\tilde{s}_{1,2}(f)$ are the Fourier transform of the time-series strain
outputs, $\hat{C}(f)$ is expected to be Gaussian distributed with expectation

$\langle\hat{C}(f)\rangle=\Gamma_{I}(f)I(f)+\Gamma_{V}(f)V(f),$ (9)

and variance (small signal-to-noise ratio limit)

$\sigma^{2}(f)\approx\frac{1}{2T\Delta f}P_{1}(f)P_{2}(f).$ (10)

Here $\Delta f$ is the width of frequency bins and $P_{1,2}(f)$ are the one-
side noise power spectrum of the detectors. $\Gamma(f)$ is the overlap
reduction function of the baseline. The overlap reduction function can be
calculated by the antenna pattern of the detector Crowder  _et al._ (2013) as
follows

$\begin{split}\Gamma_{I}(f)&=\frac{1}{8\pi}\int\mathrm{d}^{2}\Omega\left(F_{1}^{+}F_{2}^{+*}+F_{1}^{\times}F_{2}^{\times*}\right)\mathrm{e}^{2\pi\mathrm{i}f\bm{\Omega}\cdot\Delta\bm{x}},\\\
\Gamma_{V}(f)&=-\frac{\mathrm{i}}{8\pi}\int\mathrm{d}^{2}\Omega\left(F_{1}^{+}F_{2}^{\times*}-F_{1}^{\times}F_{2}^{+*}\right)\mathrm{e}^{2\pi\mathrm{i}f\bm{\Omega}\cdot\Delta\bm{x}}.\end{split}$
(11)

Note that there are three baselines denoted by H-L (LIGO Handford-LIGO
Livingston), H-V (LIGO Hanford-Virgo) and L-V (LIGO Livingston-Virgo),
respectively. The overlap functions of the interferometer pairs involved in
our analysis are plotted in Fig. 1.

Figure 1: Overlap reduction functions of H-L (LIGO Hanford-LIGO Livingston),
H-V (LIGO Hanford-Virgo) and L-V (LIGO Livingston-Virgo) pairs, respectively.

In this letter we adopt the data released by LVK Collaborations during
O1$\sim$O3 observing runs Abbott  _et al._ . The analyzed frequency band is
$20\sim 1726$ Hz with a resolution of $1/32$ Hz. For the case with more than
one baseline, the combined likelihood can be written as

$p(\bm{\hat{C}}|S)\propto\prod_{f}\exp\left[-\frac{1}{2}\left(\bm{\hat{C}}-\Gamma
S\right)^{\dagger}N^{-1}\left(\bm{\hat{C}}-\Gamma S\right)\right].$ (12)

The noise correlation matrix $N$ is diagonal with elements in Eq. (10) and

$\bm{\hat{C}}=\begin{bmatrix}\hat{C}_{1}\\\ \hat{C}_{2}\\\
\vdots\end{bmatrix},\>\Gamma=\begin{bmatrix}\Gamma_{1I}&\Gamma_{1V}\\\
\Gamma_{2I}&\Gamma_{2V}\\\ \vdots&\vdots\end{bmatrix},\>S=\begin{bmatrix}I\\\
V\end{bmatrix}.$ (13)

Maximizing the likelihood is equivalent to minimizing

$\|\bm{\hat{C}}^{\prime}-\Gamma^{\prime}S\|,$ (14)

where $\bm{\hat{C}}^{\prime}=\sqrt{N^{-1}}\bm{\hat{C}}$ and
$\Gamma^{\prime}=\sqrt{N^{-1}}\Gamma$. This least–squares problem can be
solved by applying singular value decomposition Thrane  _et al._ (2009); Press
_et al._ (2007) to $\Gamma^{\prime}$

$\Gamma^{\prime}=U\Sigma W^{\dagger},$ (15)

here $U$ and $W$ are unitary matrices, $\Sigma$ is diagonal with singular
values $w_{i}$ arranged from large to small. This decomposition defines the
_Moore-Penrose inverse_

$\Gamma^{\prime-1}=W\cdot\mathrm{diag}(1/w_{i})\cdot U^{\dagger}.$ (16)

The maximum likelihood estimator is

$\hat{S}=\Gamma^{\prime-1}\hat{\bm{C}}^{\prime},$ (17)

and the covariance becomes

$\mathrm{Cov}_{ij}=\sum_{k}\frac{W_{ik}W_{jk}}{w_{k}^{2}}.$ (18)

Here we need to remind that the magnitude of $w_{i}$ is a criterion of matrix
singularity. When $w_{i}\to 0$, the covariance matrix tends to be infinite.
This happens quite often due to the degeneration of overlap functions and
glitches in detector’s frequency band. Actually, these insensitive frequency
bins not only make little contribution to the constraint on energy spectrum,
but also make trouble to numerical calculation. Therefore, in practice, the
frequency bins with $w_{i}/w_{\text{max}}<10^{-10}$ will be removed in our
analysis. This criteria cuts off about $18.9\%$ frequency bins of
$\bm{\hat{C}}(f)$.

The estimations of Stokes’ parameters $I$ and $V$ in the frequency band
$20\sim 100$ Hz corresponding to the most sensitive band of LIGO-Virgo is
plotted in Fig. 2. These results indicate that there is no evidence for the
circularly polarized SGWB in the data of LIGO-Virgo’s first three observing
runs.

Figure 2: The estimation of Stokes’ parameters $I(f)$ and $V(f)$ with LIGO-
Virgo O1$\sim$ O3 runs. The black lines denote $1\sigma$ standard deviation
which is equal to $\sqrt{\mathrm{Cov}_{ii}(f)}$.

Power-law models.  By adopting Bayesian inference technique according to Eq.
(12), we can place constraints on various GW spectra. In the LIGO-Virgo
frequency band, for most theoretical models, $\Omega_{\text{gw}}(f)$ can be
approximated as a power-law:

$\Omega_{\text{gw}}(f)=\Omega_{\alpha}\left(\frac{f}{f_{\text{ref}}}\right)^{\alpha},$
(19)

where $\alpha$ is an index which is assumed to be a constant and the reference
frequency $f_{\text{ref}}$ is taken to be 25 Hz in this letter. In particular,
cosmic string and slow roll inflation are well approximated by $\alpha=0$ in
LIGO-Virgo frequency band Damour and Vilenkin (2000); Turner (1997); Bar-Kana
(1994), CBCs produce a spectrum with $\alpha=2/3$ Regimbau (2011), and
$\alpha=3$ is a fiducial choice because it denotes a flat spectrum of $I(f)$
and some astrophysical processes like supernovae can produce such a kind of
signal Sandick  _et al._ (2006). In addition, we introduced a new parameter

$\Pi(f)=V(f)/I(f),$ (20)

which encodes the parity violation. If $\Pi(f)\neq 0$, the parity is violated.
The range of $\Pi(f)$ is $[-1,1]$ in which the lower and upper bounds
correspond to full left or right-handed polarizations, respectively. For
simplicity, $\Pi(f)$ is taken as a constant in our analysis.

First of all, $\bm{\Theta}=(\Omega_{\alpha},\alpha)$ and $\Pi$ are taken to be
free parameters and $I(f)$ can be derived using Eq. (7). According to Bayes
theorem, the posterior is given by

$p(\bm{\Theta},\Pi|\bm{\hat{C}})\propto
p(\bm{\hat{C}}|S[\Omega(f;\bm{\Theta}),\Pi])p(\bm{\Theta},\Pi).$ (21)

The ratio of evidence, so-called Bayes factor, is a factor measuring the
relative possibility of hypotheses. Similar to Abbott  _et al._ (2019, 2021a),
we take a log-uniform prior for $\Omega_{\alpha}$ and choose the lower bound
to be $10^{-13}$. See Table 1 for the priors for other parameters. Our results
of the posterior distribution for the parameters are shown in Fig. 3. The
Bayes factor is $\log\mathcal{B}=-0.2$ between signal and pure noise
hypothesis, which indicates that there is no evidence for claiming the
existence of such signal. Besides, we have not found any significant
restriction on the polarization parameter $\Pi$. At $95\%$ confidence level
(CL), the 1D marginalized posterior of $\Omega_{\alpha}$ gives the upper limit
of $5.3\times 10^{-9}$ on the strength of SGWB at $25$ Hz.

Parameter |  Prior
---|---
$\Omega_{\alpha}$ |  $\text{LogUniform}[10^{-13},\,10^{-7}]$
$\Pi$ |  $\text{Uniform}[-1,\,1]$
$\alpha$ |  $\text{Gaussian}(0,\,3.5^{2})$
|
Table 1: Prior distributions for free parameters in the analysis. We set a
Gaussian prior with zero mean and standard deviation 3.5 for $\alpha$. Figure
3: Posterior distributions of $\alpha$, $\log\Omega_{\alpha}$ and $\Pi$. The
lower left panels are the 2D posterior density with $68\%$ and $95\%$ CL while
the diagonal subplots show the marginalized 1D posteriors. The gray dotted
lines denote priors adopted.

For some given theoretical models, $\alpha$ can be fixed. Therefore we also
provide the constraints on the models with $\alpha=0,2/3,3$, respectively. Our
results are given in Fig. 4. The left panel illustrates $95\%$ exclusion
contours in $\Omega_{\alpha}$–$\Pi$ plane. The 1D probability distributions of
$\Omega_{\alpha}$ and $\Pi$ are shown in center and right panels. To
demonstrate how the constraints on SGWB strength depend on the GW
polarization, the upper limits of $\Omega_{\alpha}$ for the unpolarized
$(\Pi=0)$ and fully polarized $(\Pi=\pm 1)$ cases are listed in Table 2. Due
to the correlation contribution from $V(f)$, we are able to put tighter
constraints on polarized SGWB, in particular in the circumstance of
$\alpha=3$. After marginalizing over the polarization degree, the upper limits
for the power-law spectra with $\alpha=0,2/3,3$ are $3.6\times 10^{-9}$,
$2.3\times 10^{-9}$ and $4.5\times 10^{-11}$ at $95\%$ CL, respectively.

Figure 4: Left: Exclusion contours of the strength $\Omega_{\alpha}$ and polarization parameter $\Pi$ for three different power-law models at $95\%$ confidence levels. Center & Right: Marginalized probability distributions of $\Omega_{\alpha}$ and $\Pi$. |  $\Pi=-1$ |  $\Pi=0$ |  $\Pi=1$ |  $\Pi$ Marg.
---|---|---|---|---
$\alpha=0$ |  $3.7\times 10^{-9}$ |  $4.0\times 10^{-9}$ |  $3.6\times 10^{-9}$ |  $3.6\times 10^{-9}$
$\alpha=2/3$ |  $1.7\times 10^{-9}$ |  $2.7\times 10^{-9}$ |  $1.7\times 10^{-9}$ |  $2.3\times 10^{-9}$
$\alpha=3$ |  $1.1\times 10^{-11}$ |  $3.7\times 10^{-10}$ |  $1.6\times 10^{-11}$ |  $4.5\times 10^{-11}$
|  |  |  |
Table 2: Upper limits of $\Omega_{\alpha}$ at $95\%$ CL for fixed values of
$\alpha$ and $\Pi$. Besides, we also list the result after marginalizing
$\Pi$.

Conclusions.  In this work, the three independent baselines in the LIGO-Virgo
network allow us to provide the first constraint on the circularly polarized
isotropic SGWB by adopting maximum likelihood estimation and Bayesian
statistics. Our results indicate that there is no evidence for the polarized
SGWB in the data of Advanced LIGO-Virgo’s first three observing runs.

By now, the ground-based detectors are not sensitive enough to claim a
detection of SGWB. Moreover, the sensitivities of LIGO-Virgo pairs are worse
than LIGO Hanford–Livingston pair. This further reduces the efficiency of
detector networks and we are unable to obtain a significant constraint on the
parity parameter $\Pi(f)$. But the extra baselines from LIGO-Virgo pairs do
help to reduce the degeneracy between $\Omega_{\alpha}$ and $\Pi$. Therefore a
lower upper limit for the strength of polarized SGWB compared to the
unpolarized SGWB Abbott  _et al._ (2021a) is obtained in our analysis.

Acknowledgments.  We acknowledge the use of HPC Cluster of ITP-CAS. This work
is supported by the National Key Research and Development Program of China
Grant No.2020YFC2201502, grants from NSFC (grant No. 11975019, 11991052,
12047503), Key Research Program of Frontier Sciences, CAS, Grant NO. ZDBS-
LY-7009, CAS Project for Young Scientists in Basic Research YSBR-006, the Key
Research Program of the Chinese Academy of Sciences (Grant NO. XDPB15).

## References

  * Rosado (2011) Pablo A. Rosado, “Gravitational wave background from binary systems,” Phys. Rev. D 84, 084004 (2011).
  * Wu  _et al._ (2012) C. Wu, V. Mandic, and T. Regimbau, “Accessibility of the gravitational-wave background due to binary coalescences to second and third generation gravitational-wave detectors,” Phys. Rev. D 85, 104024 (2012).
  * Zhu  _et al._ (2013) Xing-Jiang Zhu, Eric J. Howell, David G. Blair, and Zong-Hong Zhu, “On the gravitational wave background from compact binary coalescences in the band of ground-based interferometers,” Mon. Not. Roy. Astron. Soc. 431, 882–899 (2013), arXiv:1209.0595 [gr-qc] .
  * Ferrari  _et al._ (1999) Valeria Ferrari, Sabino Matarrese, and Raffaella Schneider, “Stochastic background of gravitational waves generated by a cosmological population of young, rapidly rotating neutron stars,” Mon. Not. Roy. Astron. Soc. 303, 258 (1999), arXiv:astro-ph/9806357 .
  * Regimbau and de Freitas Pacheco (2001) T. Regimbau and Jose A. de Freitas Pacheco, “Cosmic background of gravitational waves from rotating neutron stars,” Astron. Astrophys. 376, 381 (2001), arXiv:astro-ph/0105260 .
  * Lasky  _et al._ (2013) Paul D. Lasky, Mark F. Bennett, and Andrew Melatos, “Stochastic gravitational wave background from hydrodynamic turbulence in differentially rotating neutron stars,” Phys. Rev. D 87, 063004 (2013).
  * Crocker  _et al._ (2015) K. Crocker, V. Mandic, T. Regimbau, K. Belczynski, W. Gladysz, K. Olive, T. Prestegard, and E. Vangioni, “Model of the stochastic gravitational-wave background due to core collapse to black holes,” Phys. Rev. D 92, 063005 (2015), arXiv:1506.02631 [gr-qc] .
  * Witten (1984) Edward Witten, “Cosmic separation of phases,” Phys. Rev. D 30, 272–285 (1984).
  * Kosowsky  _et al._ (1992) Arthur Kosowsky, Michael S. Turner, and Richard Watkins, “Gravitational waves from first order cosmological phase transitions,” Phys. Rev. Lett. 69, 2026–2029 (1992).
  * Dev and Mazumdar (2016) P. S. Bhupal Dev and A. Mazumdar, “Probing the Scale of New Physics by Advanced LIGO/VIRGO,” Phys. Rev. D 93, 104001 (2016), arXiv:1602.04203 [hep-ph] .
  * Von Harling  _et al._ (2020) Benedict Von Harling, Alex Pomarol, Oriol Pujolàs, and Fabrizio Rompineve, “Peccei-Quinn Phase Transition at LIGO,” JHEP 04, 195 (2020), arXiv:1912.07587 [hep-ph] .
  * Damour and Vilenkin (2000) Thibault Damour and Alexander Vilenkin, “Gravitational wave bursts from cosmic strings,” Phys. Rev. Lett. 85, 3761–3764 (2000).
  * Sarangi and Tye (2002) Saswat Sarangi and S. H. Henry Tye, “Cosmic string production towards the end of brane inflation,” Phys. Lett. B 536, 185–192 (2002), arXiv:hep-th/0204074 .
  * Siemens  _et al._ (2007) Xavier Siemens, Vuk Mandic, and Jolien Creighton, “Gravitational-wave stochastic background from cosmic strings,” Phys. Rev. Lett. 98, 111101 (2007).
  * Abbott  _et al._ (2018) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Constraints on cosmic strings using data from the first Advanced LIGO observing run,” Phys. Rev. D 97, 102002 (2018), arXiv:1712.01168 [gr-qc] .
  * Turner (1997) Michael S. Turner, “Detectability of inflation-produced gravitational waves,” Phys. Rev. D 55, R435–R439 (1997).
  * Guzzetti  _et al._ (2016) M. C. Guzzetti, N. Bartolo, M. Liguori, and S. Matarrese, “Gravitational waves from inflation,” Riv. Nuovo Cim. 39, 399–495 (2016), arXiv:1605.01615 [astro-ph.CO] .
  * Aasi  _et al._ (2015) J Aasi, B Abbott, R Abbott, Temeka Abbott, Matthew Abernathy, K Ackley, C Adams, Teneisha Adams, Paolo Addesso, R Adhikari, Vaishali Adya, C Affeldt, Nishu Aggarwal, Odylio Aguiar, Anirban Ain, P Ajith, A Alemic, Bruce Allen, D Amariutei, and John Zweizig, “Advanced ligo,” Classical and Quantum Gravity 32 (2015), 10.1088/0264-9381/32/7/074001.
  * Acernese  _et al._ (2015) F. Acernese _et al._ (VIRGO), “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (KAGRA, Virgo, LIGO Scientific), “Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run,” Phys. Rev. D 104, 022004 (2021a), arXiv:2101.12130 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (KAGRA, Virgo, LIGO Scientific), “Search for anisotropic gravitational-wave backgrounds using data from Advanced LIGO and Advanced Virgo’s first three observing runs,” Phys. Rev. D 104, 022005 (2021b), arXiv:2103.08520 [gr-qc] .
  * Abbott  _et al._ (2022) R. Abbott _et al._ (LIGO Scientific, Virgo, KAGRA), “All-sky, all-frequency directional search for persistent gravitational waves from Advanced LIGO’s and Advanced Virgo’s first three observing runs,” Phys. Rev. D 105, 122001 (2022), arXiv:2110.09834 [gr-qc] .
  * Kahniashvili  _et al._ (2005) Tina Kahniashvili, Grigol Gogoberidze, and Bharat Ratra, “Polarized cosmological gravitational waves from primordial helical turbulence,” Phys. Rev. Lett. 95, 151301 (2005), arXiv:astro-ph/0505628 .
  * Alexander and Yunes (2009) Stephon Alexander and Nicolas Yunes, “Chern-Simons Modified General Relativity,” Phys. Rept. 480, 1–55 (2009), arXiv:0907.2562 [hep-th] .
  * Lyth  _et al._ (2005) David H. Lyth, Carlos Quimbay, and Yeinzon Rodriguez, “Leptogenesis and tensor polarisation from a gravitational Chern-Simons term,” JHEP 03, 016 (2005), arXiv:hep-th/0501153 .
  * Satoh  _et al._ (2008) Masaki Satoh, Sugumi Kanno, and Jiro Soda, “Circular Polarization of Primordial Gravitational Waves in String-inspired Inflationary Cosmology,” Phys. Rev. D 77, 023526 (2008), arXiv:0706.3585 [astro-ph] .
  * Saito  _et al._ (2007) Shun Saito, Kiyotomo Ichiki, and Atsushi Taruya, “Probing polarization states of primordial gravitational waves with CMB anisotropies,” JCAP 09, 002 (2007), arXiv:0705.3701 [astro-ph] .
  * Horava (2009) Petr Horava, “Quantum Gravity at a Lifshitz Point,” Phys. Rev. D 79, 084008 (2009), arXiv:0901.3775 [hep-th] .
  * Takahashi and Soda (2009) Tomohiro Takahashi and Jiro Soda, “Chiral primordial gravitational waves from a lifshitz point,” Phys. Rev. Lett. 102, 231301 (2009).
  * Crisostomi  _et al._ (2018) Marco Crisostomi, Karim Noui, Christos Charmousis, and David Langlois, “Beyond Lovelock gravity: Higher derivative metric theories,” Phys. Rev. D 97, 044034 (2018), arXiv:1710.04531 [hep-th] .
  * Zhao  _et al._ (2020) Wen Zhao, Tao Zhu, Jin Qiao, and Anzhong Wang, “Waveform of gravitational waves in the general parity-violating gravities,” Phys. Rev. D 101, 024002 (2020), arXiv:1909.10887 [gr-qc] .
  * Seto and Taruya (2007) Naoki Seto and Atsushi Taruya, “Measuring a Parity Violation Signature in the Early Universe via Ground-based Laser Interferometers,” Phys. Rev. Lett. 99, 121101 (2007), arXiv:0707.0535 [astro-ph] .
  * Seto and Taruya (2008) Naoki Seto and Atsushi Taruya, “Polarization analysis of gravitational-wave backgrounds from the correlation signals of ground-based interferometers: Measuring a circular-polarization mode,” Phys. Rev. D 77, 103001 (2008), arXiv:0801.4185 [astro-ph] .
  * Crowder  _et al._ (2013) S. G. Crowder, R. Namba, V. Mandic, S. Mukohyama, and M. Peloso, “Measurement of Parity Violation in the Early Universe using Gravitational-wave Detectors,” Phys. Lett. B 726, 66–71 (2013), arXiv:1212.4165 [astro-ph.CO] .
  * Allen and Romano (1999) Bruce Allen and Joseph D. Romano, “Detecting a stochastic background of gravitational radiation: Signal processing strategies and sensitivities,” Phys. Rev. D 59, 102001 (1999), arXiv:gr-qc/9710117 .
  * (36) R. Abbott _et al._ (LIGO Scientific, Virgo, KAGRA), https://dcc.ligo.org/G2001287/public.
  * Thrane  _et al._ (2009) Eric Thrane, Stefan Ballmer, Joseph D. Romano, Sanjit Mitra, Dipongkar Talukder, Sukanta Bose, and Vuk Mandic, “Probing the anisotropies of a stochastic gravitational-wave background using a network of ground-based laser interferometers,” Phys. Rev. D 80, 122002 (2009), arXiv:0910.0858 [astro-ph.IM] .
  * Press  _et al._ (2007) William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery, “Numerical recipes: the art of scientific computing, 3rd edition,” (2007).
  * Bar-Kana (1994) Rennan Bar-Kana, “Limits on direct detection of gravitational waves,” Phys. Rev. D 50, 1157–1160 (1994).
  * Regimbau (2011) Tania Regimbau, “The astrophysical gravitational wave stochastic background,” Res. Astron. Astrophys. 11, 369–390 (2011), arXiv:1101.2762 [astro-ph.CO] .
  * Sandick  _et al._ (2006) Pearl Sandick, Keith A. Olive, Frédéric Daigne, and Elisabeth Vangioni, “Gravitational waves from the first stars,” Phys. Rev. D 73, 104024 (2006).
  * Abbott  _et al._ (2019) B. P. Abbott _et al._ (LIGO Scientific, Virgo), “Search for the isotropic stochastic background using data from Advanced LIGO’s second observing run,” Phys. Rev. D 100, 061101 (2019), arXiv:1903.02886 [gr-qc] .