# Moments of derivatives of the Riemann zeta function: Characteristic
polynomials and the hybrid formula

Christopher Hughes  Department of Mathematics, University of York, York, YO10
5DD, United Kingdom<EMAIL_ADDRESS>and  Andrew Pearce-Crump
Department of Mathematics, University of York, York, YO10 5DD, United Kingdom
<EMAIL_ADDRESS>
###### Abstract.

We conjecture results about the moments of mixed derivatives of the Riemann
zeta function, evaluated at the non-trivial zeros of the Riemann zeta
function. We do this in two different ways, both giving us the same
conjecture. In the first, we find asymptotics for the moments of derivatives
of the characteristic polynomials of matrices in the Circular Unitary
Ensemble. In the second, we consider the hybrid model approach first proposed
by Gonek, Hughes and Keating.

##  1\. Introduction

Let $\zeta(s)$ be the Riemann zeta function. Throughout this paper we assume
the Riemann Hypothesis, that is, that all non-trivial zeros can be written as
$\rho=\tfrac{1}{2}+i\gamma$ with real $\gamma$. We will develop the philosophy
of Keating and Snaith [11] and Hughes, Keating and O’Connell [10] from random
matrix theory to study the discrete moments of the $n$th derivative of
$\zeta(s)$, denoted $\zeta^{(n)}(s)$. That is, we study

$\frac{1}{N(T)}\sum_{0<\gamma\leq
T}\zeta^{(n)}\left(\tfrac{1}{2}+i\gamma\right)^{k},$ (1.1)

where

$\displaystyle N(T)$ $\displaystyle=\\#\\{\text{nontrivial zeros of zeta with
}0<\gamma\leq T\\}$ $\displaystyle=\frac{T}{2\pi}\log\frac{T}{2\pi e}+O(\log
T).$ (1.2)

In fact, we will study the more general problem of moments of mixed
derivatives,

$\frac{1}{N(T)}\sum_{0<\gamma\leq
T}\zeta^{(n_{1})}\left(\tfrac{1}{2}+i\gamma\right)\ldots\zeta^{(n_{k})}\left(\tfrac{1}{2}+i\gamma\right).$

Clearly if $n_{1}=\dots=n_{k}=n$, then we recover (1.1).

The model proposed by Keating and Snaith [11] is that the characteristic
polynomial of a large random unitary matrix can be used to model the value
distribution of the Riemann zeta function near a large height T. They used the
characteristic polynomial

$Z(\theta)=\det(I-Ue^{-i\theta})=\prod_{m=1}^{N}(1-e^{i(\theta_{m}-\theta)}),$
(1.3)

where the $\theta_{m}$ are the eigenangles of an $N\times N$ random unitary
matrix $U$, to model $\zeta(s)$. The motivation for studying this is that it
has been conjectured that the limiting distribution of the non-trivial zeros
of $\zeta(s)$, on the scale of their mean spacing, is asymptotically the same
as that of the eigenangles $\theta_{m}$ of matrices chosen according to Haar
measure, in the limit as $N\to\infty$. Equating the mean densities yields the
connection between matrix size and height up the critical line

$N=\log\frac{T}{2\pi}.$

Averaging over all Haar-distributed $N\times N$ unitary matrices, Keating and
Snaith [11] calculated the moments of $|Z(\theta)|$ and found that as
$N\to\infty$

$\mathbb{E}_{N}\left[|Z(\theta)|^{2k}\right]\sim\frac{G^{2}(k+1)}{G(2k+1)}N^{k^{2}}$

where $\mathbb{E}_{N}$ denotes expectation with respect to Haar measure, and
$G(z)$ is the Barnes G-function.

Comparing with all known and conjectured results for continuous moments of
zeta, they were lead to conjecture the following result.

######  Conjecture 1 (Keating–Snaith).

For any fixed $k$ with $\Re(k)>-1/2$,

$\frac{1}{T}\int_{0}^{T}\left|\zeta\left(\tfrac{1}{2}+it\right)\right|^{2k}\
dt\sim a(k)\frac{G^{2}(k+1)}{G(2k+1)}\left(\log\frac{T}{2\pi}\right)^{k^{2}},$

as $T\rightarrow\infty$, and

$a(k)=\prod_{p\text{
prime}}\left(1-\frac{1}{p}\right)^{k^{2}}\sum_{m=0}^{\infty}\left(\frac{\Gamma(m+k)}{m!\,\Gamma(k)}\right)^{2}p^{-m}$
(1.4)

is an arithmetic factor which is discussed below.

The Keating and Snaith result is for a continuous moment, but random matrix
theory has also been used to predict discrete moments of $\zeta(s)$. Hughes,
Keating and O’Connell [10] followed a similar approach and calculated the
discrete moments of $|Z^{\prime}(\theta)|$ and found that as $N\to\infty$

$\mathbb{E}_{N}\left[|Z^{\prime}(\theta)|^{2k}\right]\sim\frac{G^{2}(k+2)}{G(2k+3)}N^{k(k+2)}.$

As in the continuous case above, motivated by this result they made a
conjecture which agreed with all known and previously conjectured results and
is believed to hold in general.

######  Conjecture 2 (Hughes–Keating–O’Connell).

For any fixed $k$ with $\Re(k)>-3/2$,

$\frac{1}{N(T)}\sum_{0<\gamma\leq
T}\left|\zeta^{\prime}\left(\tfrac{1}{2}+i\gamma\right)\right|^{2k}\sim
a(k)\frac{G^{2}(k+2)}{G(2k+3)}\left(\log\frac{T}{2\pi}\right)^{k(k+2)},$

where $a(k)$ is the same arithmetic factor as in the previous conjecture,
given in (1.4).

What wasn’t clear at the time was why the same arithmetic factor (1.4) should
appear in these two conjectures. Indeed, the arithmetic factor itself was
added in an ad hoc manner after the random matrix theory calculations were
computed. This additional factor has been shown to be included naturally in
the hybrid model developed by Gonek, Hughes and Keating [7] for continuous
moments and of Bui, Gonek and Milinovich [3] for discrete moments. This hybrid
approach gives equal weighting to $Z_{X}$, a partial Hadamard product over the
zeros of zeta, and $P_{X}$, a partial Euler product over the primes. It is
stated exactly in Theorem 7 in Section 3. Another strength of the hybrid
approach is that it explained exactly why the same factor should appear in
both conjectures.

The previous two problems where random matrix theory was used both concerned
moments of a real function, where an absolute value of zeta was taken. For our
problem, we consider the moments of the complex function. In this paper we
attack our problem in two ways. First, by direct calculation of the
characteristic polynomial, and then by modelling zeta using the hybrid
approach.

In Section 2 we will calculate

$\mathbb{E}_{N}\left[\frac{1}{N}\sum_{m=1}^{N}Z(\theta_{m}+\alpha_{1})\ldots
Z(\theta_{m}+\alpha_{k})\right]$

where, as before, $\mathbb{E}_{N}$ denotes expectation with respect to Haar
measure. By Taylor expanding each characteristic polynomial, we will be able
to prove the following theorem.

######  Theorem 3.

For fixed, positive integers $n_{1},\ldots,n_{k}$, we have

$\mathbb{E}_{N}\left[\frac{1}{N}\sum_{m=1}^{N}Z^{(n_{1})}(\theta_{m})\ldots
Z^{(n_{k})}(\theta_{m})\right]\\\ \sim(-1)^{n_{1}+\ldots+n_{k}+k}\
i^{n_{1}+\ldots+n_{k}}\frac{n_{1}!\ldots
n_{k}!}{(n_{1}+\ldots+n_{k}+1)!}N^{n_{1}+\ldots+n_{k}},$

as $N\rightarrow\infty$.

Setting $n_{1}=\dots=n_{k}=n$ gives the following Corollary.

######  Corollary 3.1.

For $n,k\in\mathbb{N}$, as $N\rightarrow\infty$

$\mathbb{E}_{N}\left[\frac{1}{N}\sum_{m=1}^{N}Z^{(n)}(\theta_{m})^{k}\right]\sim(-1)^{k(n+1)}\
i^{kn}\frac{(n!)^{k}}{(kn+1)!}N^{kn}.$

In order to use this to create a conjecture for zeta, we need to know what the
arithmetic term equals. To do that, we take the hybrid approach, discussed in
Section 3. In Section 4 (Theorems 8 and 9) we calculate the discrete moments
of $P_{X}^{(n_{1})}\dots P_{X}^{(n_{k})}$ and show that they asymptotically
vanish unless all the $n_{j}$ are $0$ (that is, no derivatives are taken), in
which case the mean is 1. Furthermore, we predict that the discrete moments of
$Z_{X}$ are accurately modelled by the characteristic polynomial, and we will
verify this in the case of the first derivative in Theorem 13 in Section 5).
That is, we conjecture that there is no arithmetic term at the leading order
for these moments of the complex Riemann zeta function.

Note that the derivative of $Z(\theta)$ is with respect to $\theta$, and this
acts like the derivative of $\zeta(1/2+it)$ with respect to $t$. This is
obviously equal to $i\zeta^{\prime}(1/2+it)$, and so we expect there to be no
factors of $i$ for the Riemann zeta function. To be precise, we make the
following prediction:

######  Conjecture 4.

For fixed, positive integers $n_{1},\ldots,n_{k}$, the moments of mixed
derivatives of $\zeta(s)$, evaluated at the non-trivial zeros of $\zeta(s)$,
are given by

$\frac{1}{N(T)}\sum_{0<\gamma\leq
T}\zeta^{(n_{1})}\left(\tfrac{1}{2}+i\gamma\right)\ldots\zeta^{(n_{k})}\left(\tfrac{1}{2}+i\gamma\right)\\\
\sim(-1)^{n_{1}+\ldots+n_{k}+k}\frac{n_{1}!\ldots
n_{k}!}{(n_{1}+\ldots+n_{k}+1)!}\left(\log\frac{T}{2\pi}\right)^{n_{1}+\ldots+n_{k}}$

as $T\rightarrow\infty$.

Setting $n_{1}=\dots=n_{k}=n$ gives the following conjecture.

######  Conjecture 5.

For $n,k\in\mathbb{N}$, the $k$th moment of the $n$th derivative of $\zeta(s)$
is given by

$\frac{1}{N(T)}\sum_{0<\gamma\leq
T}\zeta^{(n)}\left(\tfrac{1}{2}+i\gamma\right)^{k}\sim(-1)^{k(n+1)}\frac{(n!)^{k}}{(kn+1)!}\left(\log\frac{T}{2\pi}\right)^{kn}$

as $T\rightarrow\infty$.

Very little is rigorously known about these types of results. For $k=0$ the
result is trivial, as the sum is just $N(T)$, defined in equation (1). For
$k=1$, they have been previously studied in relation to Shanks’ Conjecture, a
full history of which can be found in [8], and the results agree with
Conjecture 5 in this case.

These results are for leading order only. In a forthcoming paper [9] we use a
different method to conjecture the full asymptotic for zeta.

These results and conjectures are stated for $k$ a positive integer. However,
in the case of the first derivative, when $n=1$, Section 5 applies for
$k\in\mathbb{C}$ with $\Re(k)>-3$, and some of our results in Section 4
(namely Theorem 9) also apply for non-integer $k$. Therefore it is plausible
that in the case of the first derivative, Conjecture 5 also holds for complex
$k$.

######  Conjecture 6.

For $\Re(k)>-3$ as $T\to\infty$

$\frac{1}{N(T)}\sum_{0<\gamma\leq
T}\zeta^{\prime}\left(\tfrac{1}{2}+i\gamma\right)^{k}\sim\frac{1}{\Gamma(k+2)}\left(\log\frac{T}{2\pi}\right)^{k}.$

In the case $k=-1$, this result is already known and can be found as equation
(3.9) in a paper of Garaev and Sankaranarayanan, [5].

When $k=-2$, the right hand side of Conjecture 6 is zero. A separate
calculation shows that

$\mathbb{E}_{N}\left[\frac{1}{N}\sum_{m=1}^{N}Z^{\prime}(\theta_{m})^{-2}\right]=0$

exactly. For the case of zeta, what we mean is the left hand side is bounded
by some error terms with no leading-order asymptotic. An educated guess
suggests the error terms will be $O(T^{-1/2+\epsilon})$.

##  2\. Proof of Theorem 3

In this section we prove Theorem 3 concerning moments of characteristic
polynomials. We start with the product representation of the characteristic
polynomial of an $N\times N$ random unitary matrix $Z(\theta)$ in (1.3), given
by

$Z(\theta)=\prod_{m=1}^{N}(1-e^{i(\theta_{m}-\theta)}).$

The Haar probability density of $U(N)$ equals [13]

$d\mu_{N}=\frac{1}{N!(2\pi)^{N}}\prod_{1\leq j<\ell\leq
N}|e^{i\theta_{j}}-e^{i\theta_{\ell}}|^{2}\ d\theta_{1}\ldots d\theta_{N}.$
(2.1)

To extract the moments that we wish to calculate, we use shifts and calculate

$\mathbb{E}_{N}\left[\frac{1}{N}\sum_{m=1}^{N}Z(\theta_{m}+\alpha_{1})\ldots
Z(\theta_{m}+\alpha_{k})\right]$

By invariance of Haar measure there are no distinguished eigenvalues, so this
expectation equals

$\mathbb{E}_{N}\left[Z(\theta_{N}+\alpha_{1})\ldots
Z(\theta_{N}+\alpha_{k})\right]$

Since we can Taylor expand assuming the $\alpha_{j}$ are very small, we obtain

$\mathbb{E}_{N}\left[\left(Z^{\prime}(\theta_{N})\frac{\alpha_{1}}{1!}+Z^{\prime\prime}(\theta_{N})\frac{\alpha_{1}^{2}}{2!}+\ldots+Z^{(n_{1})}(\theta_{N})\frac{\alpha_{1}^{n_{1}}}{n_{1}!}+\ldots\right)\right.\\\
\times\ldots\times\left.\left(Z^{\prime}(\theta_{N})\frac{\alpha_{k}}{1!}+\ldots+Z^{(n_{k})}(\theta_{N})\frac{\alpha_{k}^{n_{k}}}{n_{k}!}+\ldots\right)\right],$

where obviously $Z(\theta_{N})=0$. By calculating the
$\alpha_{1}^{n_{1}}\ldots\alpha_{k}^{n_{k}}$-th coefficient, we can find
moments of the mixed derivatives of the characteristic polynomial.

Writing out the expectation in terms of the Weyl density (2.1), we have

$\mathbb{E}_{N}\left[\prod_{r=1}^{k}Z(\theta_{N}+\alpha_{r})\right]=\\\
\frac{1}{N!(2\pi)^{N}}\idotsint_{0}^{2\pi}\prod_{1\leq j<\ell\leq
N}|e^{i\theta_{j}}-e^{i\theta_{\ell}}|^{2}\prod_{m=1}^{N}\prod_{r=1}^{k}(1-e^{i\theta_{m}}e^{-i\theta_{N}}e^{-i\alpha_{r}})\
d\theta_{m}.$

We split off the terms in the first product with $\ell=N$, and all the terms
in the second product with $m=N$ (which can then be pulled out of the
integral) to give

$\frac{1}{N!(2\pi)^{N}}\prod_{r=1}^{k}(1-e^{-i\alpha_{r}})\idotsint_{0}^{2\pi}\prod_{1\leq
j<\ell\leq N-1}|e^{i\theta_{j}}-e^{i\theta_{\ell}}|^{2}\\\
\times\prod_{m=1}^{N-1}\left|e^{i\theta_{m}}-e^{i\theta_{N}}\right|^{2}\prod_{r=1}^{k}(1-e^{i\theta_{m}}e^{-i\theta_{N}}e^{-i\alpha_{r}})\
d\theta_{m}$

We may make a change of variables $\phi_{m}=\theta_{m}-\theta_{N}$ for
$m=1,\dots,N-1$ and trivially perform the $\theta_{N}$ integral to obtain

$\frac{1}{N!(2\pi)^{N-1}}\prod_{r=1}^{k}(1-e^{-i\alpha_{r}})\idotsint_{0}^{2\pi}\prod_{1\leq
j<\ell\leq N-1}\left|e^{i\phi_{j}}-e^{i\phi_{\ell}}\right|^{2}\\\
\times\prod_{m=1}^{N-1}(1-e^{i\phi_{m}})(1-e^{-i\phi_{m}})\prod_{r=1}^{k}(1-e^{i\phi_{m}}e^{-i\alpha_{r}})\
d\phi_{m}$

Notice then that we can rewrite the integrals back in terms of Haar measure
for $(N-1)\times(N-1)$ unitary matrices, obtaining

$\frac{1}{N}\prod_{r=1}^{k}(1-e^{-i\alpha_{r}})\times\mathbb{E}_{N-1}\left[\prod_{m=1}^{N-1}(1-e^{i\phi_{m}})(1-e^{-i\phi_{m}})\prod_{r=1}^{k}(1-e^{i\phi_{m}}e^{-i\alpha_{r}})\right]\\\
=\frac{(-1)^{(k+1)(N-1)}}{N}\prod_{r=1}^{k}A_{r}^{-N}\left(A_{r}-1\right)\times\mathbb{E}_{N-1}\left[\prod_{m=1}^{N-1}e^{-i\phi_{m}}(e^{i\phi_{m}}-1)^{2}\prod_{r=1}^{k}(e^{i\phi_{m}}-A_{r})\right]$

where we have done some simple algebraic manipulations to put it into a
suitable form for later, and written $A_{r}=e^{i\alpha_{r}}$. We may then use
Heine’s Lemma [12]111Though the identity was first written down in 1939 by
Szegő, he gave the credit to Heine (who lived from 1821 to 1881). to write the
expectation as a Toeplitz determinant, and thus we have

$\mathbb{E}_{N-1}\left[\prod_{m=1}^{N-1}e^{-i\phi_{m}}(e^{i\phi_{m}}-1)^{2}\prod_{r=1}^{k}(e^{i\phi_{m}}-A_{r})\right]=D_{N-1}[f]$

where

$D_{N-1}[f]:=\det\left|\hat{f}_{j-\ell}\right|_{1\leq j,\ell\leq N-1}$

with

$\hat{f}_{j-\ell}=\frac{1}{2\pi}\int_{0}^{2\pi}f(\vartheta)e^{-i(j-\ell)\vartheta}\mathrm{d}\vartheta$

for the symbol

$f(z)=z^{-1}(z-1)^{2}\prod_{r=1}^{k}\left(z-A_{r}\right).$ (2.2)

Putting this together, we have

$\mathbb{E}_{N}\left[\prod_{r=1}^{k}Z(\theta_{N}+\alpha_{r})\right]=\frac{(-1)^{(k+1)(N-1)}}{N}\prod_{r=1}^{k}A_{r}^{-N}\left(A_{r}-1\right)\times
D_{N-1}[f]$ (2.3)

Basor and Forrester [2] evaluate $D_{N-1}[f]$ with symbols of the form (2.2)
as a $(k+2)\times(k+2)$ determinant,

$D_{N-1}[f]=\frac{(-1)^{(k+1)(N-1)}\prod_{r=1}^{k+2}(-A_{r})}{\prod_{1\leq
j<\ell\leq k+2}(A_{\ell}-A_{j})}\\\
\times\det\begin{bmatrix}A_{1}^{N-1}&A_{1}^{N}&A_{1}^{N+1}&\dots&A_{1}^{N-1+k}&(-A_{1})^{-1}\\\
A_{2}^{N-1}&A_{2}^{N}&A_{2}^{N+1}&\dots&A_{2}^{N-1+k}&(-A_{2})^{-1}\\\
A_{3}^{N-1}&A_{3}^{N}&A_{3}^{N+1}&\dots&A_{3}^{N-1+k}&(-A_{3})^{-1}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
A_{k}^{N-1}&A_{k}^{N}&A_{k}^{N+1}&\dots&A_{k}^{N-1+k}&(-A_{k})^{-1}\\\
A_{k+1}^{N-1}&A_{k+1}^{N}&A_{k+1}^{N+1}&\dots&A_{k+1}^{N-1+k}&(-A_{k+1})^{-1}\\\
A_{k+2}^{N-1}&A_{k+2}^{N}&A_{k+2}^{N+1}&\dots&A_{k+2}^{N-1+k}&(-A_{k+2})^{-1}\\\
\end{bmatrix},$

where we write $A_{k+1}$ and $A_{k+2}$ for now, but will let them both tend to
$1$ later.

We can manipulate this determinant through elementary row operations. We begin
by factoring $A_{r}^{N-1}$ from each row, and factor the $-1$ from the last
column. This determinant becomes

$D_{N-1}[f]=\frac{(-1)^{N(k+1)}\prod_{r=1}^{k+2}A_{r}^{N}}{\prod_{1\leq
j<\ell\leq
k+2}(A_{\ell}-A_{j})}\det\begin{bmatrix}1&A_{1}&A_{1}^{2}&\dots&A_{1}^{k}&A_{1}^{-N}\\\
1&A_{2}&A_{2}^{2}&\dots&A_{2}^{k}&A_{2}^{-N}\\\
1&A_{3}&A_{3}^{2}&\dots&A_{3}^{k}&A_{3}^{-N}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&A_{k}&A_{k}^{2}&\dots&A_{k}^{N-1+k}&A_{k}^{-N}\\\
1&A_{k+1}&A_{k+1}^{2}&\dots&A_{k+1}^{k}&A_{k+1}^{-N}\\\
1&A_{k+2}&A_{k+2}^{2}&\dots&A_{k+2}^{k}&A_{k+2}^{-N}\\\ \end{bmatrix}.$

Subtract the first row from each other row, and bring the factor of
$1/(A_{r}-A_{1})$ inside the determinant for each row $r=2,\dots k+2$ to give

$D_{N-1}[f]=\frac{(-1)^{N(k+1)}\prod_{r=1}^{k+2}A_{r}^{N}}{\prod_{2\leq
j<\ell\leq
k+2}(A_{\ell}-A_{j})}\det\begin{bmatrix}1&A_{1}&A_{1}^{2}&\dots&A_{1}^{k}&A_{1}^{-N}\\\
0&1&A_{1}+A_{2}&\dots&\displaystyle\sum_{j_{1}=0}^{k-1}A_{1}^{j_{1}}A_{2}^{k-1-j_{1}}&C_{2,2}\\\
0&1&A_{1}+A_{3}&\dots&\displaystyle\sum_{j_{1}=0}^{k-1}A_{1}^{j_{1}}A_{3}^{k-1-j_{1}}&C_{3,2}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&1&A_{1}+A_{k}&\dots&\displaystyle\sum_{j_{1}=0}^{k-1}A_{1}^{j_{1}}A_{k}^{k-1-j_{1}}&C_{k,2}\\\
0&1&A_{1}+A_{k+1}&\dots&\displaystyle\sum_{j_{1}=0}^{k-1}A_{1}^{j_{1}}A_{k+1}^{k-1-j_{1}}&C_{k+1,2}\\\
0&1&A_{1}+A_{k+2}&\dots&\displaystyle\sum_{j_{1}=0}^{k-1}A_{1}^{j_{1}}A_{k+2}^{k-1-j_{1}}&C_{k+2,2}\\\
\end{bmatrix}$

where in the $r$th row and $m$th column we make use of the following identity

$\frac{A_{r}^{m-1}-A_{1}^{m-1}}{A_{r}-A_{1}}=\sum_{j_{1}=0}^{m-2}A_{1}^{j_{1}}A_{r}^{m-2-j_{1}}.$
(2.4)

for $m$ from $2$ to $k+1$, and for the $(k+2)$nd column we use the identity

$C_{r,2}:=\frac{A_{r}^{-N}-A_{1}^{-N}}{A_{r}-A_{1}}=-\sum_{j_{1}=1}^{N}A_{1}^{-j_{1}}A_{r}^{-(N+1-j_{1})}.$
(2.5)

Now subtract the second row from each other row below it and bring the factor
of $1/(A_{r}-A_{2})$ inside the determinant for each row $r=3,\dots k+2$ to
give $D_{N-1}[f]$ as

$\frac{(-1)^{N(k+1)}\prod_{r=1}^{k+2}A_{r}^{N}}{\prod_{3\leq j<\ell\leq
k+2}(A_{\ell}-A_{j})}\det\begin{bmatrix}1&A_{1}&A_{1}^{2}&\dots&A_{1}^{k}&A_{1}^{-N}\\\
0&1&A_{1}+A_{2}&\dots&\displaystyle\sum_{j_{1}=0}^{k-1}A_{1}^{j_{1}}A_{2}^{k-1-j_{1}}&C_{2,2}\\\
0&1&1&\dots&\displaystyle\sum_{j_{1}=0}^{k-2}\sum_{j_{2}=0}^{k-2-j_{1}}A_{1}^{j_{1}}A_{2}^{j_{2}}A_{3}^{k-2-j_{1}-j_{2}}&C_{3,3}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&1&1&\dots&\displaystyle\sum_{j_{1}=0}^{k-2}\sum_{j_{2}=0}^{k-2-j_{1}}A_{1}^{j_{1}}A_{2}^{j_{2}}A_{k}^{k-2-j_{1}-j_{2}}&C_{k,3}\\\
0&1&1&\dots&\displaystyle\sum_{j_{1}=0}^{k-2}\sum_{j_{2}=0}^{k-2-j_{1}}A_{1}^{j_{1}}A_{2}^{j_{2}}A_{k+1}^{k-2-j_{1}-j_{2}}&C_{k+1,3}\\\
0&1&1&\dots&\displaystyle\sum_{j_{1}=0}^{k-2}\sum_{j_{2}=0}^{k-2-j_{1}}A_{1}^{j_{1}}A_{2}^{j_{2}}A_{k+2}^{k-2-j_{1}-j_{2}}&C_{k+2,3}\\\
\end{bmatrix}$

where we have used (2.4) for the $r$th row and $m$th column for $m=3,...,k+1$
and where the last column uses the identity (2.5) with $A_{2}$ instead of
$A_{1}$ with

$\displaystyle C_{r,3}$
$\displaystyle=-\sum_{j_{1}=1}^{N}A_{1}^{-j_{1}}\left(\frac{A_{r}^{-(N+1-j_{1})}-A_{2}^{-(N+1-j_{1})}}{A_{r}-A_{2}}\right)$
$\displaystyle=(-1)^{2}\sum_{j_{1}=1}^{N}\sum_{j_{2}=1}^{N+1-j_{1}}A_{1}^{-j_{1}}A_{2}^{-j_{2}}A_{r}^{-(N+2-j_{1}-j_{2})}.$

Continuing in this manner, it is clear that $D_{N-1}[f]$ becomes

$D_{N-1}[f]=(-1)^{N(k+1)}\prod_{r=1}^{k+2}A_{r}^{N}\\\
\times\det\begin{bmatrix}1&A_{1}&A_{1}^{2}&\dots&A_{1}^{k}&A_{1}^{-N}\\\
0&1&A_{1}+A_{2}&\dots&\displaystyle\sum_{j_{1}=0}^{k-1}A_{1}^{j_{1}}A_{2}^{k-1-j_{1}}&C_{2,2}\\\
0&0&1&\dots&\displaystyle\sum_{j_{1}=0}^{k-2}\sum_{j_{2}=0}^{k-2-j_{1}}A_{1}^{j_{1}}A_{2}^{j_{2}}A_{3}^{k-2-j_{1}-j_{2}}&C_{3,3}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\dots&A_{1}+A_{2}+\dots+A_{k}&C_{k,k}\\\ 0&0&0&\dots&1&C_{k+1,k+1}\\\
0&0&0&\dots&0&C_{k+2,k+2}\\\ \end{bmatrix},$

that is, we have an upper-triangular matrix with $1$’s on the diagonal and
zeros underneath. This means that the determinant equals the bottom right-most
element, which is

$C_{k+2,k+2}=(-1)^{(k+1)}\sum_{j_{1}=1}^{N}\sum_{j_{2}=1}^{N+1-j_{1}}\dots\sum_{j_{k}=1}^{N+k-1-j_{1}-\dots-
j_{k-1}}A_{1}^{-j_{1}}\dots A_{k}^{-j_{k}}\\\
\times\sum_{j_{k+1}=1}^{N+k-j_{1}-\dots-
j_{k}}A_{k+1}^{j_{k+1}}A_{k+2}^{-(N+k+2-j_{1}-\dots-j_{k+1})}$

We now combine $D_{N-1}[f]$ with the other terms in (2.3), and let
$A_{k+1},A_{k+2}\rightarrow 1$. This gives

$\mathbb{E}_{N}\left[\prod_{r=1}^{k}Z(\theta_{N}+\alpha_{r})\right]\\\
=\frac{1}{N}\prod_{r=1}^{k}\left(A_{r}-1\right)\sum_{j_{1}=1}^{N}\sum_{j_{2}=1}^{N+1-j_{1}}\dots\sum_{j_{k}=1}^{N+k-1-j_{1}-\dots-
j_{k-1}}A_{1}^{-j_{1}}\dots A_{k}^{-j_{k}}\sum_{j_{k+1}=1}^{N+k-j_{1}-\dots-
j_{k}}1.$

Re-indexing the summations to start from $0$ instead of $1$ and evaluating the
last sum, we have

$\mathbb{E}_{N}\left[\prod_{r=1}^{k}Z(\theta_{N}+\alpha_{r})\right]\\\
=\frac{1}{N}\prod_{r=1}^{k}A_{r}^{-1}\left(A_{r}-1\right)\sum_{j_{1}=0}^{N-1}\sum_{j_{2}=0}^{N-1-j_{1}}\dots\sum_{j_{k}=0}^{N-1-j_{1}-\dots-
j_{k-1}}A_{1}^{-j_{1}}\dots A_{k}^{-j_{k}}(N-j_{1}-\dots-j_{k}).$ (2.6)

Recall that the coefficient of $\alpha_{1}^{n_{1}}\dots\alpha_{k}^{n_{k}}$,
which we denote by $[\alpha_{1}^{n_{1}}\dots\alpha_{k}^{n_{k}}]$, from (2.6)
will yield the desired mixed derivative, and we wish to find this to leading
order as $N\to\infty$. It is clear that to do this, we want to obtain the
highest possible power of $N$ from the summations as they depend upon $N$, and
the smallest possible contribution from the product
$\prod_{r=1}^{k}A_{r}^{-1}\left(A_{r}-1\right)$ as there is no $N$ dependence
in that product.

Recalling that $A_{r}=e^{i\alpha_{r}}$ and Taylor expanding gives

$\prod_{r=1}^{k}A_{r}^{-1}\left(A_{r}-1\right)=i^{k}\alpha_{1}\dots\alpha_{k}+\text{
higher order terms.}$ (2.7)

Therefore we need to find the large $N$ behaviour of
$[\alpha_{1}^{n_{1}-1}\dots\alpha_{k}^{n_{k}-1}]$ from the summations in
(2.6), that is from

$\sum_{j_{1}=0}^{N-1}\sum_{j_{2}=0}^{N-1-j_{1}}\dots\sum_{j_{k}=0}^{N-1-j_{1}-\dots-
j_{k-1}}e^{-ij_{1}\alpha_{1}}\dots e^{-ij_{k}\alpha_{k}}(N-j_{1}-\dots-
j_{k}).$

The coefficient of $\alpha_{1}^{n_{1}-1}$ in the Taylor expansion of
$e^{-ij_{1}\alpha_{1}}$ equals $(-ij_{1})^{n_{1}-1}/(n_{1}-1)!$, so the
coefficient $[\alpha_{1}^{n_{1}-1}...\alpha_{k}^{n_{k}-1}]$ equals

$\sum_{j_{1}=0}^{N-1}\dots\sum_{j_{k}=0}^{N-1-j_{1}-\dots-
j_{k-1}}\frac{(-ij_{1})^{n_{1}-1}}{(n_{1}-1)!}\dots\frac{(-ij_{k})^{n_{k}-1}}{(n_{k}-1)!}(N-j_{1}-\dots-
j_{k})\\\
=\frac{(-i)^{n_{1}+\dots+n_{k}-k}}{(n_{1}-1)!\dots(n_{k}-1)!}\sum_{j_{1}=0}^{N-1}\dots\sum_{j_{k}=0}^{N-1-j_{1}-\dots-
j_{k-1}}j_{1}^{n_{1}-1}\dots j_{k}^{n_{k}-1}(N-j_{1}-\dots-j_{k}).$

To evaluate the leading order coefficient of these sums, we will turn them
into Riemann sums and hence into an integral. The sums equal

$N^{n_{1}+\dots+n_{k}+1}\sum_{j_{1}=0}^{N-1}\sum_{j_{2}=0}^{N-1-j_{1}}\dots\sum_{j_{k}=0}^{N-1-j_{1}-\dots-
j_{k-1}}\frac{1}{N^{k}}\left(\frac{j_{1}}{N}\right)^{n_{1}-1}\dots\left(\frac{j_{k}}{N}\right)^{n_{k}-1}\left(1-\frac{j_{1}}{N}-\dots-\frac{j_{k}}{N}\right)$

and now each summation can be considered as a Riemann sum and so we can
rewrite the summation in the above line as the $k$-fold integral

$\int_{0}^{1}\int_{0}^{1-x_{1}}\dots\int_{0}^{1-x_{1}-\dots-
x_{k-1}}x_{1}^{n_{1}-1}\dots x_{k}^{n_{k}-1}(1-x_{1}-\dots-x_{k})\
dx_{1}dx_{2}\dots dx_{k}.$

By Theorem 1.8.6 in [1], this integral equals

$\frac{(n_{1}-1)!\dots(n_{k}-1)!}{(n_{1}+\ldots+n_{k}+1)!}$

and so we have

$\displaystyle[\alpha_{1}^{n_{1}-1}...\alpha_{k}^{n_{k}-1}]$
$\displaystyle\sim\frac{(-i)^{n_{1}+\dots+n_{k}-k}}{(n_{1}-1)!\dots(n_{k}-1)!}\frac{(n_{1}-1)!\ldots(n_{k}-1)!}{(n_{1}+\ldots+n_{k}+1)!}N^{n_{1}+\ldots+n_{k}+1}$
(2.8)

Combining (2.7) with (2.8) and the term $1/N$ in (2.6), we obtain the
coefficient $[\alpha_{1}^{n_{1}}...\alpha_{k}^{n_{k}}]$ of
$\mathbb{E}_{N}\left[Z(\theta_{N}+\alpha_{1})\dots
Z(\theta_{N}+\alpha_{k})\right]$ is asymptotic to

$\frac{(-1)^{n_{1}+\dots+n_{k}+k}\
i^{n_{1}+\dots+n_{k}}}{(n_{1}+\ldots+n_{k}+1)!}N^{n_{1}+\ldots+n_{k}}.$

Recalling that the coefficient of $\alpha_{r}^{n_{r}}$ in the Taylor expansion
of $Z(\theta_{N}+\alpha_{r})$ equals $Z^{(n_{r})}(\theta_{N})/n_{r}!$, we see
that

$\mathbb{E}_{N}\left[Z^{(n_{1})}(\theta_{N})\dots
Z^{(n_{k})}(\theta_{N})\right]\sim(-1)^{n_{1}+\dots+n_{k}+k}\
i^{n_{1}+\dots+n_{k}}\frac{n_{1}!\dots
n_{k}!}{(n_{1}+\ldots+n_{k}+1)!}N^{n_{1}+\ldots+n_{k}},$

and this completes the proof of Theorem 3.

##  3\. The hybrid formula

In 2007 Gonek, Hughes and Keating [7], proved a hybrid formula for the zeta
function, expressing it as a partial product over primes times a rapidly
decaying product over non-trivial zeros. We state a result with slightly more
control over the smoothing (in [7] and other later papers, $X=Y$ is taken).

######  Theorem 7.

Let $s=\sigma+it$ with $\sigma\geq 0$ and $|t|\geq 2$, let $X,Y\geq 2$ be real
parameters, and let $K$ be any fixed positive integer. Let $f(x)$ be a
nonnegative $C^{\infty}$ function of mass 1 supported on $[0,1]$ and set
$u(x)=Yf(Y\log\frac{x}{e}+1)/x$. Thus, $u(x)$ is a function of mass 1
supported on $[e^{1-1/Y},e]$. Set

$U(z)=\int_{1}^{e}u(y)E_{1}(z\log y)\;\mathrm{d}y\,,$ (3.1)

where $E_{1}(z)$ is the exponential integral
$\int_{z}^{\infty}e^{-w}/w\;\mathrm{d}w$. Then

$\zeta(s)=P_{X}(s)Z_{X}(s)\left(1+O\left(\frac{Y^{K+1}X^{\max(1-\sigma,0)}}{(|s|\log
X)^{K}}\right)+O\left(\frac{X^{1-\sigma}\log X}{Y}\right)\right)\;,$

where

$P_{X}(s)=\exp\left(\sum_{n\leq X}\frac{\Lambda(n)}{\log
n}\frac{1}{n^{s}}\right)\,,$ (3.2)

$\Lambda(n)$ is von Mangoldt’s function, and

$Z_{X}(s)=\exp\left(-\sum_{\rho_{n}}U\big{(}(s-\rho_{n})\log
X\big{)}\right)\,.$ (3.3)

The constants implied by the $O$ terms depend only on $f$ and $K$.

######  Remark.

The proof of this result follows precisely the method in [7], with obvious
changes for the different smoothing.

Note that $U(z)$ is not an analytic function; it has a logarithmic singularity
at $z=0$. However, recall the formula

$E_{1}(z)=-\log z-\gamma-\sum_{m=1}^{\infty}\frac{(-1)^{m}z^{m}}{m!\,m}\,,$
(3.4)

where $|\arg z|<\pi$, $\log z$ denotes the principal branch of the logarithm,
and $\gamma$ is Euler’s constant. (It is clear that the sum is an entire
function of $z$ since it is absolutely convergence for all $z\in\mathbb{C}$).
From this and (3.1) we observe that we may interpret $\exp(-U(z))$ to be an
analytic function, asymptotic to $Cz$ as $z\to 0$, for some constant $C$ which
depends upon the smoothing function $u$.

Just as in work of Bui, Gonek and Milinovich [3], we can differentiate whilst
still obtaining the asymptotic, obtaining under the Riemann Hypothesis

$\zeta^{\prime}(\rho)=P_{X}(\rho)Z_{X}^{\prime}(\rho)\left(1+O\left(\frac{Y^{K+1}X^{1/2}}{(|\rho|\log
X)^{K}}\right)+O\left(\frac{X^{1/2}\log X}{Y}\right)\right)$

It is believed that $P_{X}$ and $Z_{X}$ operate pseudo-independently. Hence
the moments of zeta are products of moments of $P_{X}$ and $Z_{X}$. This is
known as the Splitting Conjecture.

In the next section we will evaluate the moments of mixed derivatives of
$P_{X}$, when averaged over zeros of zeta. We will show that the moments only
really contribute when no derivatives are taken. In the last section we will
model the moments of the first derivative of $Z_{X}$ when averaged over zeros
of zeta, using a method that allows non-integer moments to be taken. Taken
collectively, these results give weight to Conjectures 5 and 6.

##  4\. The moments of $P_{X}$

In this section we will show that the mixed moments of $P_{X}$ vanish if any
derivatives are taken, whilst the mean of $P_{X}^{k}$ is asymptotically 1,
when averaged over zeros of the zeta function.

######  Theorem 8.

Under the Riemann Hypothesis, for $X>2$ with $X=O(\log T)$, for any non-
negative integers $n_{1},\dots,n_{k}$ we have,

$\sum_{0<\gamma\leq T}P_{X}^{(n_{1})}(\rho)\dots
P_{X}^{(n_{k})}(\rho)=\begin{cases}N(T)+O(T\log\log T)&\text{ if
}n_{1}=\dots=n_{k}=0\\\ O\left(T(\log\log T)^{1+n_{1}}\right)&\text{ if
}n_{1}>0\text{ and }n_{2}=\dots=n_{k}=0\\\ O\left(T\right)&\text{ if more than
one $n_{i}>0$}\end{cases}$

as $T\to\infty$, where $\rho=\tfrac{1}{2}+i\gamma$ is a non-trivial zero of
the Riemann zeta function, $P_{X}(s)$ is given in (3.2) and $N(T)$ is the
number of such zeros up to height $T$, given by (1).

In the case that no derivatives are taken, one can permit $k$ to be real and
positive.

######  Theorem 9.

Under the Riemann Hypothesis, for $X>2$ with $X=O(\log T)$, for $k$ a fixed
positive real number,

$\sum_{0<\gamma\leq T}P_{X}(\rho)^{k}=N(T)+O(T\log\log T)$

as $T\to\infty$.

This should be compared to a result of Bui, Gonek and Milinovich (Theorem 2.2
of [3]) which stated that, under the Riemann Hypothesis, for $\epsilon>0$ and
$X,T\to\infty$ with $X=O\left((\log T)^{2-\epsilon}\right)$ then for any
$k\in\mathbb{R}$ one has

$\frac{1}{N(T)}\sum_{0<\gamma\leq
T}\left|P_{X}(\rho)\right|^{2k}=a(k)\left(e^{\gamma_{0}}\log
X\right)^{k^{2}}\left(1+O_{k}\left((\log X)^{-1}\right)\right)$

where $a(k)$ is an explicit product over primes given by (1.4).

Before we prove Theorem 8, we need to state and prove some preliminary lemmas.

######  Lemma 10.

For $k$ a fixed positive real number, the function $P_{X}(s)^{k}$ can be
written as a non-vanishing absolutely convergent Dirichlet series

$P_{X}(s)^{k}=\sum_{m=1}^{\infty}\frac{a_{k}(m)}{m^{s}}$

with $0\leq a_{k}(m)\leq d_{k}(m)$ where $d_{k}(m)$ is the $k$th divisor
function.

###### Proof.

Since $P_{X}(s)$ given in (3.2) is the exponential of a finite Dirichlet
polynomial, so its $k$th power is an entire non-vanishing function, so can be
written as a Dirichlet series

$\displaystyle P_{X}(s)^{k}$ $\displaystyle=\exp\left(k\sum_{n\leq
X}\frac{\Lambda(n)}{\log n}\frac{1}{n^{s}}\right)$
$\displaystyle=\sum_{m=1}^{\infty}\frac{a_{k}(m)}{m^{s}}$

where the sum converges absolutely for any $s$, and $a_{k}(m)=0$ if $m$ has
any prime divisor greater than $X$. If $p\leq X$ is a prime, then
$a_{k}(p)=k$. Since these coefficients come from the exponential of positive
terms, the $a_{k}(m)$ are positive. In the case when $m\leq X$, these
coefficients will equal $d_{k}(m)$ and in the case when $m>X$ they will be
less than $d_{k}(m)$ due to missing contributions from terms with $n>X$ in the
exponent. ∎

In the proof of the two theorems, it will be important to truncate the
infinite sum.

######  Lemma 11 (Gonek-Hughes-Keating).

For $2\leq X\ll(\log T)^{2-\epsilon}$ with $\epsilon>0$, and $k$ a fixed
positive real number, and $\theta$ a small positive number, we have

$P_{X}(\tfrac{1}{2}+it)^{k}=\sum_{m\leq
T^{\theta}}\frac{a_{k}(m)}{m^{1/2+it}}+O_{k}\left(T^{-\epsilon\theta/2}\right)$

###### Proof.

This follows immediately from the proof of Lemma 2 in [7]. ∎

In [6] Gonek proved unconditionally a uniform version of Landau’s formula
concerning sums over the non-trivial zeros of the Riemann zeta function,
including the following corollary which requires the assumption of the Riemann
Hypothesis.

######  Lemma 12 (Gonek).

Under the Riemann Hypothesis, for $T>1,m\in\mathbb{N}$ with $m\geq 2$ and
$\rho$ a non-trivial zero of the Riemann zeta function $\zeta(s)$,

$\sum_{0<\gamma\leq
T}m^{-\rho}=-\frac{T}{2\pi}\frac{\Lambda(m)}{m}+O(\log(2mT)\log\log(3m)),$

where $\Lambda(m)$ is the von Mangoldt function.

######  Remark.

Note this result does not apply in the case when $m=1$, when the sum over
zeros trivially equals $N(T)$.

###### Proof of Theorem 8.

To differentiate $P_{X}(\sigma)$ $j$ times, we return to the exponential form
(3.2). The Faà di Bruno formula on multiple derivatives of $e^{f}$ states

$\frac{d^{j}}{d\sigma^{j}}e^{f(\sigma)}=\left(f^{(j)}+\dots+(f^{\prime})^{j}\right)e^{f}$

where the sum ranges over all possible combinations of derivatives of $f$ so
that total derivative is $j$, with certain multinomial coefficients that don’t
matter for us.

In our case, one can check that the longest sum comes from the first
derivative raised to the $j$th power, essentially giving a sum of length
$X^{j}$. That is

$P_{X}^{(j)}(s)=\sum_{n\leq
X^{j}}\frac{\tilde{\beta}_{j}(n)}{n^{s}}\exp\left(\sum_{n\leq
X}\frac{\Lambda(n)}{\log n}\frac{1}{n^{s}}\right)$

for coefficients $\tilde{\beta}_{j}(n)$ that come from appropriate
combinations of derivatives of the sum in the exponential, whose precise
nature does not concern us, other than noting that

$\tilde{\beta}_{j}(1)=\begin{cases}0&\text{ if }j\neq 0\\\ 1&\text{ if
}j=0\end{cases}$

and for prime $p\leq X$,

$\tilde{\beta}_{j}(p)=\begin{cases}(-\log p)^{j}&\text{ if }j\neq 0\\\
0&\text{ if }j=0\end{cases}$

(this comes from differentiating the sum in the exponential $j$ times; other
terms in the Faà di Bruno formula contribute to coefficients of higher powers
of $p$). If $n$ is divisible by any prime greater than $X$,
$\tilde{\beta}_{j}(n)=0$.

Combining $k$ of these expressions together (replacing $j$ with derivatives
$n_{1},n_{2},\dots,n_{k}$, we have

$P_{X}^{(n_{1})}(s)\dots P_{X}^{(n_{k})}(s)=\sum_{n\leq
X^{n_{1}+n_{2}+\dots+n_{k}}}\frac{\beta(n)}{n^{s}}P_{X}(s)^{k}$ (4.1)

where

$\beta(n)=(\tilde{\beta}_{n_{1}}\ast\tilde{\beta}_{n_{2}}\ast\dots\ast\tilde{\beta}_{n_{k}})(n)$

is the convolution of the appropriate $\tilde{\beta}_{j}$.

In the case when $\Re(s)=1/2$ and $X=(\log T)^{2-\epsilon}$, using Lemma 11 we
have

$\displaystyle P_{X}^{(n_{1})}(s)\dots P_{X}^{(n_{k})}(s)$
$\displaystyle=\sum_{n\leq
X^{n_{1}+n_{2}+\dots+n_{k}}}\frac{\beta(n)}{n^{s}}\left(\sum_{m\leq
T^{\theta}}\frac{a_{k}(m)}{m^{s}}+O_{k}\left(T^{-\epsilon\theta/2}\right)\right)$
$\displaystyle=\sum_{m\leq
T^{\theta^{\prime}}}\frac{b(m)}{m^{s}}+O\left(T^{-\epsilon^{\prime}}\right)$
(4.2)

for some $\epsilon^{\prime}>0$ and $\theta^{\prime}>\theta$, which both depend
upon $\epsilon$, $k$, $\theta$, and $n_{1},\dots,n_{k}$. Essentially, since
the first sum’s length is only a power of $\log T$, the two sums together are
truncated at $T^{\theta^{\prime}}$ for some $\theta^{\prime}>\theta$, as $T$
gets large.

Shortly we will use some facts about the $b(m)$, which we record here:

  * •

Only if $n_{1}=n_{2}=\dots=n_{k}=0$ does $b(1)=1$, otherwise $b(1)=0$
(essentially because $\tilde{\beta}_{j}(1)=0$ for any $j\neq 0$ since
differentiation kills constants).

  * •

If $n_{1}=n_{2}=\dots=n_{k}=0$, then $b(p)=a_{k}(p)$ and for prime $p\leq X$,
$a_{k}(p)=k$.

  * •

If exactly all bar one of the $n_{i}$ vanish (say $n_{1}\neq 0$, and
$n_{2}=n_{3}=\dots=n_{k}=0$) then
$b(p)=\beta(p)=\tilde{\beta}_{n_{1}}(p)=(-\log p)^{n_{1}}$. To see this, note
that $\tilde{\beta}_{n_{1}}(1)=0$, which kills all other contributions in the
convolution other than $\tilde{\beta}_{n_{1}}(p)$.

  * •

If more than one $n_{i}\neq 0$, then $b(p)=0$ since there will be at least two
pieces in the convolution that have $\tilde{\beta}_{n_{1}}(1)=0$ and
$\tilde{\beta}_{n_{2}}(1)=0$, so every term vanishes.

Assuming the Riemann Hypothesis, we can sum (4.1) over the zeros of zeta,
obtaining

$\sum_{0<\gamma\leq T}P_{X}^{(n_{1})}(\rho)\dots
P_{X}^{(n_{k})}(\rho)=b(1)N(T)+\sum_{0<\gamma\leq T}\sum_{2\leq m\leq
T^{\theta^{\prime}}}\frac{b(m)}{m^{\rho}}+O\left(N(T)T^{-\epsilon^{\prime}}\right)$
(4.3)

and applying Lemma 12 to the middle sum yields

$\sum_{0<\gamma\leq T}\sum_{2\leq m\leq
T^{\theta}}\frac{b(m)}{m^{\rho}}=-\frac{T}{2\pi}\sum_{2\leq m\leq
T^{\theta^{\prime}}}\frac{b(m)\Lambda(m)}{m}+\sum_{2\leq m\leq
T^{\theta}}b(m)O(\log(2mT)\log\log(3m)).$ (4.4)

To bound the first term on the RHS of (4.4), note that the $\Lambda(m)$ forces
the sum to be over primes and prime powers only, and the $b(m)$ ensures we
only need primes $\leq X$. Since $b(m)\ll m^{\epsilon}$, the square- and
higher-powers of primes form a convergent sum, so

$-\frac{T}{2\pi}\sum_{2\leq m\leq
T^{\theta^{\prime}}}\frac{b(m)\Lambda(m)}{m}=O\left(T\sum_{\begin{subarray}{c}p\leq
X\\\ p\text{ prime}\end{subarray}}\frac{b(p)\log p}{p}\right)+O(T)$

The value of $b(p)$ depends on how many $n_{i}=0$

$b(p)=\begin{cases}k&\text{ if }n_{1}=\dots=n_{k}=0\\\ (-\log
p)^{n_{1}}&\text{ if }n_{1}>0\text{ and }n_{i}=0\text{ otherwise}\\\ 0&\text{
if more than one $n_{i}>0$}.\end{cases}$

Hence, summing over all primes up to $X$ and using the Prime Number Theorem,

$\sum_{p\leq x}\frac{b(p)\log p}{p}\ll\begin{cases}\log X&\text{ if
}n_{1}=\dots=n_{k}=0\\\ (\log X)^{1+n_{1}}&\text{ if }n_{1}>0\text{ and
}n_{i}=0\text{ otherwise}\\\ 0&\text{ if more than one $n_{i}>0$}.\end{cases}$

For the second error term in (4.4) since the $b(m)$ all have the same sign
(for any given fixed $n_{1},\dots,n_{k}$), we have this equals

$O\left(\log T\log\log T\sum_{m=1}^{\infty}b(m)\right)$

By (4.2) and (4.1) the sum is simply

$\sum_{m=1}^{\infty}b(m)=\sum_{n\leq
X^{n_{1}+n_{2}+\dots+n_{k}}}\beta(n)P_{X}(0)^{k}$

To bound the first sum on the right-hand side, recall that $\beta(n)\ll
n^{\epsilon}$, and to bound the $P_{X}(0)^{k}$ note that since

$\displaystyle\log P_{X}(0)$ $\displaystyle=\sum_{n\leq
X}\frac{\Lambda(n)}{\log n}$ $\displaystyle=\frac{X}{\log
X}+O\left(\frac{X}{(\log X)^{2}}\right)$

by the Prime Number Theorem we have

$\sum_{m=1}^{\infty}b(m)=O\left(\exp\left(c^{\prime}\frac{X}{\log
X}\right)\right)$

for any fixed $c^{\prime}>k$.

Plugging these two error terms back into (4.4), we have shown

$\sum_{0<\gamma\leq T}\sum_{2\leq m\leq
T^{\theta}}\frac{b(m)}{m^{\rho}}=O\left(T(\log
X)^{1+n_{1}}\right)+O(T)+O\left(\log T\log\log
T\exp\left(c^{\prime}\frac{X}{\log X}\right)\right)$

where the first error term is only present if $n_{2}=\dots=n_{k}=0$.

To complete the proof we now choose $X$ to balance the various error terms
that appear when we substitute this into (4.3). If we insist $X=O(\log T)$
then we have

$\sum_{0<\gamma<T}P_{X}^{(n_{1})}(\rho)\dots
P_{X}^{(n_{k})}(\rho)=\frac{T}{2\pi}\log T+O\left(T(\log\log T)\right)$

in the case when $n_{1}=\dots=n_{k}=0$,

$\sum_{0<\gamma<T}P_{X}^{(n_{1})}(\rho)\dots
P_{X}^{(n_{k})}(\rho)=O\left(T(\log\log T)^{1+n_{1}}\right)$

if $n_{1}>0$ and $n_{2}=\dots=n_{k}=0$, and

$\sum_{0<\gamma<T}P_{X}^{(n_{1})}(\rho)\dots
P_{X}^{(n_{k})}(\rho)=O\left(T\right)$

if more than one $n_{i}>0$.

This completes the proof of Theorem 8. ∎

###### Proof of Theorem 9.

In the case when no derivatives are taken, we can exploit the multiplicative
nature of $P_{X}$ to allow $k$ to be a fixed real number, not just an integer.
Much of the proof goes through unchanged.

By Lemma 11 and Lemma 12 we have

$\displaystyle\sum_{0<\gamma\leq T}P_{X}(\rho)^{k}$
$\displaystyle=N(T)+\sum_{0<\gamma\leq T}\sum_{2\leq m\leq
T^{\theta}}\frac{a_{k}(m)}{m^{\rho}}+O\left(N(T)T^{-\epsilon^{\prime}}\right)$
$\displaystyle=N(T)-\frac{T}{2\pi}\sum_{2\leq m\leq
T^{\theta}}\frac{a_{k}(m)\Lambda(m)}{m}+O\left(\log T\log\log T\sum_{2\leq
m\leq T^{\theta}}a_{k}(m)\right)$

since $a_{k}(1)=1$.

As in the case of Theorem 8, assuming $X<T^{\theta}$ we can bound

$\displaystyle\sum_{2\leq m\leq T^{\theta}}\frac{a_{k}(m)\Lambda(m)}{m}$
$\displaystyle=\sum_{p<X}\frac{k\log p}{p}+O(1)$ $\displaystyle\ll\log X$

and, exactly as before,

$\sum_{m=1}^{\infty}a_{k}(m)=P_{X}(0)^{k}\ll\exp\left(c^{\prime}\frac{X}{\log
X}\right)$

for any fixed $c^{\prime}>k$.

Therefore

$\sum_{0<\gamma\leq T}P_{X}(\rho)^{k}=N(T)+O\left(T\log X\right)+O\left(\log
T\log\log T\exp\left(c^{\prime}\frac{X}{\log X}\right)\right)$

As before, the error terms can be made to approximately match when $X=O(\log
T)$. ∎

######  Remark.

In the previous proofs we can take $X$ to be anything up to $o(\log T\log\log
T)$ before the second error term dominates the main term of
$\frac{1}{2\pi}T\log T$. Since the error term as we have written it comes from
a subsidiary main term, we cannot reduce the error down to $O(T)$.

##  5\. The moments of $Z_{X}^{\prime}$

In this section we will concentrate on the random matrix equivalent for
$Z_{X}$, given in (3.3), and for simplicity just look at its first derivative.
Again in this section, we will not require $k$ to be an integer. Taking an
$N\times N$ unitary matrix with eigenvalues
$e^{i\theta_{1}},\dots,e^{i\theta_{N}}$, we wish to calculate

$\mathbb{E}_{N}\left[Z_{N,X}^{\prime}(\theta_{N})^{k}\right]$

where

$Z_{N,X}(\theta)=\prod_{m=1}^{N}\exp\left(-\sum_{j=-\infty}^{\infty}U\big{(}i(\theta-\theta_{m}+2\pi
j)\log X\big{)}\right)$

plays the role of $Z_{X}(s)$ in Theorem 7, and where $U$ is given in (3.1).

######  Remark.

The sum over $j$ makes the arguments $2\pi$-periodic, thus permitting the move
from $e^{i\theta}$ to $\theta$ without any branch-cut ambiguities. The decay
of $U$ for large inputs means that the infinite sum over $j$ is absolutely
summable.

######  Theorem 13.

If $k\in\mathbb{C}$ is such that $k\not\in\\{-3,-4,-5,\dots\\}$ then

$\mathbb{E}_{N}\left[\frac{1}{N}\sum_{m=1}^{N}Z_{N,X}^{\prime}(\theta_{m})^{k}\right]\sim\frac{e^{ik\pi/2}}{\Gamma(k+2)}N^{k}$

as $N\to\infty$.

######  Remark.

The proof deals with the analytic continuation of the objects under
consideration. However, for our purposes as a model for zeta, we must stop
before the first singularity; that is, in Conjecture 6 we take $\Re(k)>-3$.

######  Remark.

Restricting to $k\in\mathbb{N}$, we see this Theorem has the same right hand
side as Corollary 3.1 in the case $n=1$. That is, $Z^{\prime}_{N,X}$ and the
first derivative of the characteristic polynomial asymptotically have the same
moments.

###### Proof.

It is convenient to factor out the behaviour around the eigenvalues (where
$Z_{N,X}$ vanishes) as follows

$Z_{N,X}(\theta)=\prod_{m=1}^{N}\left(1-e^{-i(\theta-\theta_{m})}\right)e^{F_{X}(\theta-\theta_{m})}$

where

$F_{X}(\vartheta)=-\log\left(1-e^{-i\vartheta}\right)-\sum_{j=-\infty}^{\infty}U\big{(}i(\vartheta+2\pi
j)\log X\big{)}$ (5.1)

is a $2\pi$-periodic function which is continuously differentiable for all
real $\vartheta$ — the logarithmic singularities cancel out. The first
statement is obvious. The second statement can be seen from the definition of
$U(z)$ given in (3.1) and the series expansion of $E_{1}(z)$ around $z=0$
given in (3.4) (and indeed, this fact will be made clear when we look at the
decay of its Fourier coefficients in Lemma 14).

Differentiating,

$Z_{N,X}^{\prime}(\theta)=\sum_{n=1}^{N}ie^{-i(\theta-\theta_{n})}e^{F_{X}(\theta-\theta_{n})}\prod_{\begin{subarray}{c}m=1\\\
m\neq
n\end{subarray}}^{N}\left(1-e^{-i(\theta-\theta_{m})}\right)e^{F_{X}(\theta-\theta_{m})}\\\
+\sum_{n=1}^{N}\left(1-e^{-i(\theta-\theta_{n})}\right)F_{X}^{\prime}(\theta-\theta_{n})e^{F_{X}(\theta-\theta_{n})}\prod_{\begin{subarray}{c}m=1\\\
m\neq
n\end{subarray}}^{N}\left(1-e^{-i(\theta-\theta_{m})}\right)e^{F_{X}(\theta-\theta_{m})}$

and substituting $\theta=\theta_{N}$, we see that for every $n\neq N$ there is
a term in the product that vanishes (so only the $n=N$ terms survives), and
every term in the second term vanishes, including the $n=N$ term. That is,

$Z_{N,X}^{\prime}(\theta_{N})=ie^{F_{X}(0)}\prod_{m=1}^{N-1}\left(1-e^{-i(\theta_{N}-\theta_{m})}\right)e^{F_{X}(\theta_{N}-\theta_{m})}$

######  Remark.

There is nothing special about picking $\theta_{N}$ over any of the other
eigenvalues. This is made for notational convenience. The true calculation
would be to evaluate

$\frac{1}{N}\sum_{m=1}^{N}Z_{N,X}^{\prime}(\theta_{m})$

similar to what we do in Section 2, but by rotation invariance the result will
be the same.

We now calculate the expected value of $Z_{N,X}^{\prime}(\theta_{N})^{k}$ for
$k\in\mathbb{C}$ when averaged over all $N\times N$ unitary matrices
distributed according to Haar measure, employing a trick found in [10].

$\displaystyle\mathbb{E}_{N}\left[Z_{N,X}^{\prime}(\theta_{N})^{k}\right]$
$\displaystyle=e^{ik\pi/2}e^{kF_{X}(0)}\mathbb{E}_{N}\left[\prod_{m=1}^{N-1}\left(1-e^{-i(\theta_{N}-\theta_{m})}\right)^{k}e^{kF_{X}(\theta_{N}-\theta_{m})}\right]$
$\displaystyle=e^{ik\pi/2}e^{kF_{X}(0)}\frac{1}{N}\mathbb{E}_{N-1}\left[\prod_{m=1}^{N-1}\left|1-e^{i\vartheta_{m}}\right|^{2}\left(1-e^{i\vartheta_{m}}\right)^{k}e^{kF_{X}(-\vartheta_{m})}\right]$
(5.2)

where we interpret $e^{i\vartheta_{m}}$ as the eigenvalues of an
$(N-1)\times(N-1)$ unitary matrix. (Briefly, what we have done is write the
first expectation out as a $N$-dimensional Weyl integral, then change
variables to $\vartheta_{m}=\theta_{m}-\theta_{N}$ for $m=1,\dots,N-1$. Making
use of the fact the integrand is $2\pi$-periodic, we turn it back into a
($N-1$)-dimension Weyl integral, plus an extra trivial integral over
$\theta_{N}$).

Because $k$ is complex, we need to be precise about the choice of branch
taken. We define $\left(1-e^{i\vartheta}\right)^{k}$ to be the value obtained
from continuous variation of $\left(1-e^{i\vartheta}e^{-\epsilon}\right)^{k}$
for $\epsilon>0$, converging to the value $1$ as $\epsilon\to\infty$.

As in Section 2, the inner expectation can be calculated by Heine’s identity
to be equal to the Toeplitz determinant

$D_{N-1}[f]:=\det\left|\hat{f}_{j-\ell}\right|_{1\leq j,\ell\leq N-1}$

with symbol

$f(\vartheta)=\left|1-e^{i\vartheta}\right|^{2}\left(1-e^{i\vartheta}\right)^{k}e^{kF_{X}(-\vartheta)}$

Note that $f(\vartheta)$ has only one singularity at $\vartheta=0$. Toeplitz
determinants with such symbols have been calculated asymptotically by Ehrhardt
and Silbermann in [4] (building on extensive previous work by several other
authors).

They show (in their Theorem 2.5) that if $b(\vartheta)$ is a suitably smooth
(defined in their paper) $2\pi$-periodic function with winding number zero and
whose logarithm has the Fourier series expansion

$\log b(\vartheta)=\sum_{m=-\infty}^{\infty}s_{m}e^{im\vartheta}$

then for symbols of the form

$f(\vartheta)=b(\vartheta)\left(1-e^{i\vartheta}\right)^{\gamma}\left(1-e^{-i\vartheta}\right)^{\delta}$

subject to $\gamma+\delta\not\in\\{-1,-2,-3,\dots\\}$, then for any
$\epsilon>0$

$D_{N-1}[f]=C_{1}N^{\gamma\delta}C_{2}^{N-1}\left(1+O\left(\frac{1}{N^{1-\epsilon}}\right)\right)$

where

$C_{1}=\frac{G(1+\gamma)G(1+\delta)}{G(1+\gamma+\delta)}\exp\left(\sum_{m=1}^{\infty}ms_{m}s_{-m}\right)\exp\left(-\delta\sum_{m=1}^{\infty}s_{m}\right)\exp\left(-\gamma\sum_{m=1}^{\infty}s_{-m}\right)$

where $G$ is the Barnes $G$-function, and where

$C_{2}=\exp(s_{0}).$

In our case, $\gamma=k+1$ and $\delta=1$ and
$b(\vartheta)=\exp(kF_{X}(-\vartheta))$. To complete the proof, we simply need
to calculate the Fourier coefficients of $kF_{X}(-\vartheta)$ (those
coefficients being the desired $s_{m}$).

We will, in fact, delay the exact calculation of $s_{m}$ to Lemma 14 to
quickly complete the proof of the theorem. All we need from the lemma is that
$s_{m}=0$ if $m\leq 0$. From (5.2) and using Ehrhardt and Silbermann’s formula
above, with the values $\gamma=k+1$ and $\delta=1$, we have

$\mathbb{E}_{N}\left[Z_{N,X}^{\prime}(\theta_{N})^{k}\right]\sim
e^{ik\pi/2}e^{kF_{X}(0)}\frac{1}{N}\frac{G(k+2)G(2)}{G(k+3)}\exp\left(-\sum_{m=1}^{\infty}s_{m}\right)N^{k+1}$

Finally, note that $G(2)=1$, $G(k+3)=\Gamma(k+2)G(k+2)$ and

$e^{kF_{X}(0)}=\exp\left(\sum_{m=-\infty}^{\infty}s_{m}\right)=\exp\left(\sum_{m=1}^{\infty}s_{m}\right)$

(the first equality following from setting $\vartheta=0$ in the Fourier series
expansion; the second equality comes from knowing $s_{m}=0$ for $m\leq 0$).
Therefore,

$\mathbb{E}_{N}\left[Z_{N,X}^{\prime}(\theta_{N})^{k}\right]\sim\frac{e^{ik\pi/2}}{\Gamma(k+2)}N^{k}$

as $N\to\infty$, as required. ∎

######  Lemma 14.

Let $F_{X}(\vartheta)$ be defined in (5.1), then

$kF_{X}(-\vartheta)=\sum_{m=-\infty}^{\infty}s_{m}e^{im\vartheta}$

where

$s_{m}=\begin{cases}0&\text{ if }m\leq 0\\\
\frac{k}{m}\left(1-\int_{1}^{\exp(m/\log X)}u(y)\mathrm{d}y\right)&\text{ if
}1\leq m<\log X\\\ 0&\text{ if }m\geq\log X\end{cases}$

######  Remark.

Since $u$ has total mass 1, for $m>0$ we can write

$1-\int_{1}^{\exp(m/\log X)}u(y)\mathrm{d}y=\int_{\exp(m/\log
X)}^{\infty}u(y)\mathrm{d}y$

although due to the support condition on $u$, if $m\geq\log X$ then
$\exp(m/\log X)>e$ and so the integral vanishes.

###### Proof of Lemma 14.

Note that

$s_{m}=\frac{1}{2\pi}\int_{-\pi}^{\pi}kF_{X}(-\vartheta)e^{-im\vartheta}\mathrm{d}\vartheta.$

From (5.1) we have

$kF_{X}(-\vartheta)=-k\log\left(1-e^{i\vartheta}\right)-k\sum_{j=-\infty}^{\infty}U\big{(}i(-\vartheta+2\pi
j)\log X\big{)}.$

For $m$ an integer, we have

$\frac{1}{2\pi}\int_{-\pi}^{\pi}-k\log\left(1-e^{i\vartheta}\right)e^{-im\vartheta}\mathrm{d}\vartheta=\begin{cases}0&\text{
if }m\leq 0\\\ \frac{k}{m}&\text{ if }m\geq 1\end{cases}$ (5.3)

(To non-rigorously see why, note that
$\log\left(1-e^{i\vartheta}\right)=-\sum_{\ell=1}^{\infty}e^{i\ell\vartheta}/\ell$
and the only term in the sum that survives the integration is when $\ell=m$
for positive integers $m$).

Furthermore, we have

$\frac{1}{2\pi}\int_{-\pi}^{\pi}-k\sum_{j=-\infty}^{\infty}U\left(i(-\vartheta+2\pi
j)\log
X\right)e^{-im\vartheta}\mathrm{d}\vartheta=\frac{-k}{2\pi}\int_{-\infty}^{\infty}U(-i\vartheta\log
X)e^{-im\vartheta}\mathrm{d}\vartheta\\\
=\frac{-k}{2\pi}\int_{-\infty}^{\infty}\int_{1}^{e}u(y)E_{1}(-i\vartheta\log
X\log y)e^{-im\vartheta}\;\mathrm{d}y\mathrm{d}\vartheta$

where the final equality comes from inserting the definition of $U$ from (3.1)
and the support of $u$. Swapping the order of integration, this is

$\frac{-k}{2\pi}\int_{1}^{e}u(y)\int_{-\infty}^{\infty}E_{1}(-i\vartheta\log
X\log y)e^{-im\vartheta}\;\mathrm{d}\vartheta\mathrm{d}y.$

The Fourier transform of the exponential integral is

$\frac{1}{2\pi}\int_{-\infty}^{\infty}E_{1}(-i\vartheta\log X\log
y)e^{-im\vartheta}\;\mathrm{d}\vartheta=\begin{cases}0&\text{ if }m<\log X\log
y\\\ \frac{1}{2m}&\text{ if }m=\log X\log y\\\ \frac{1}{m}&\text{ if }m>\log
X\log y\end{cases}$

Therefore, making use of the fact that $u$ is supported on $[1,e]$ and has
total mass 1, we have

$\frac{-k}{2\pi}\int_{1}^{e}u(y)\int_{-\infty}^{\infty}E_{1}(-i\vartheta\log
X\log
y)e^{-im\vartheta}\;\mathrm{d}\vartheta\mathrm{d}y=\frac{-k}{m}\int_{1}^{\max(1,\exp(m/\log
X))}u(y)\mathrm{d}y\\\ =\begin{cases}0&\text{ if }m\leq 0\\\
-\frac{k}{m}\int_{1}^{\exp(m/\log X)}u(y)\mathrm{d}y&\text{ if }1\leq m<\log
X\\\ -\frac{k}{m}&\text{ if }m\geq\log X\end{cases}$

The Fourier coefficient, $s_{m}$ is the sum of this and (5.3). Note that both
terms are zero for $m\leq 0$ and the two terms perfectly cancel each other if
$m\geq\log X$. ∎

## References

  * [1] G.E. Andrews, R. Askey, and R. Roy, _Special functions_ , Encyclopedia of Mathematics and its Applications, vol. 71, Cambridge University Press, 1999.
  * [2] E. L. Basor and P. J. Forrester, _Formulas for the evaluation of Toeplitz determinants with rational generating functions_ , Math. Nachr. 170 (1994), 5–18.
  * [3] H. M. Bui, Steven M. Gonek, and Micah B. Milinovich, _A hybrid Euler-Hadamard product and moments of $\zeta^{\prime}(\rho)$_, Forum Mathematicum 27 (2015), 1799–1828.
  * [4] T. Ehrhardt and B. Silbermann, _Toeplitz determinants with one Fisher–Hartwig singularity_ , Journal of Functional Analysis 148 (1997), 229–256.
  * [5] M. Z. Garaev and A. Sankaranarayanan, _The sum involving derivative of $\zeta(s)$ over simple zeros_, J. Number Theory 117 (2006), 122–130.
  * [6] S. M. Gonek, _An explicit formula of Landau and its applications to the theory of the zeta-function_ , A tribute to Emil Grosswald: number theory and related analysis, Contemp. Math., vol. 143, Amer. Math. Soc., Providence, RI, 1993, pp. 395–413.
  * [7] S. M. Gonek, C. P. Hughes, and J. P. Keating, _A hybrid Euler-Hadamard product for the Riemann zeta function_ , Duke Math. J. 136 (2007), 507–549.
  * [8] C. Hughes, G. Martin, and A. Pearce-Crump, _A heuristic for discrete mean values of the derivative of the Riemann zeta function_ , Preprint (arXiv:2305.14253), 2023.
  * [9] C. Hughes and A. Pearce-Crump, _Moments of derivatives of the Riemann zeta function: The Ratios Conjecture_ , (in preparation).
  * [10] C. P. Hughes, J. P. Keating, and N. O’Connell, _Random matrix theory and the derivative of the Riemann zeta function_ , Proc. R. Soc. Lond. A 456 (2000), 2611–2627.
  * [11] J. P. Keating and N. C. Snaith, _Random matrix theory and $\zeta(1/2+it)$_, Comm. Math. Phys. 214 (2000), 57–89.
  * [12] G. Szegő, _Orthogonal Polynomials_ , AMS Colloquium Publications XXII, 1939.
  * [13] H. Weyl, _Classical Groups_ , Princeton University Press, 1946.


# Phase control of escapes in the fractional damped Helmholtz oscillator

Mattia Coccolo  Nonlinear Dynamics, Chaos and Complex Systems Group,
Departamento de Física, Universidad Rey Juan Carlos, Tulipán s/n, 28933
Móstoles, Madrid, Spain  Jesús M. Seoane  Nonlinear Dynamics, Chaos and
Complex Systems Group, Departamento de Física, Universidad Rey Juan Carlos,
Tulipán s/n, 28933 Móstoles, Madrid, Spain  Stefano Lenci  Department of Civil
and Building Engineering, and Architecture, Polytechnic University of Marche,
60131 Ancona, Italy  Miguel A.F. Sanjuán  Nonlinear Dynamics, Chaos and
Complex Systems Group, Departamento de Física, Universidad Rey Juan Carlos,
Tulipán s/n, 28933 Móstoles, Madrid, Spain

###### Abstract

We analyze the nonlinear Helmholtz oscillator in the presence of fractional
damping, a characteristic feature in several physical situations. In our
specific scenario, as well as in the non-fractional case, for large enough
excitation amplitudes, all initial conditions are escaping from the potential
well. To address this, we incorporate the phase control technique into a
parametric term, a feature commonly encountered in real-world situations. In
the non-fractional case it has been shown that, a phase difference of
$\phi_{OPT}\approx\pi$, is the optimal value to avoid the escapes of the
particles from the potential well. Here, our investigation focuses on
understanding when particles escape, considering both the phase difference
$\phi$ and the fractional parameter $\alpha$ as control parameters. Our
findings unveil the robustness of phase control, as evidenced by the
consistent oscillation of the optimal $\phi$ value around its non-fractional
counterpart when varying the fractional parameter. Additionally, our results
underscore the pivotal role of the fractional parameter in governing the
proportion of bounded particles, even when utilizing the optimal phase.

##  I Introduction

Phase control is a crucial concept for complex dynamical systems. Nonlinear
oscillators exhibit behaviors that defy simple periodicity, often showcasing
rich and intricate dynamics. The control of their dynamical behaviors and,
therefore, obtaining a desirable situation or avoiding an undesirable
situation, is a relevant topic in the context of Nonlinear Dynamics. There are
different ways of applying a control on a given dynamical system, that
typically can be classified in two main types depending on how they interact
with the system: feedback and non-feedback methods BOCCALETTI:PR . Among the
non-feedback methods, the phase control technique Zambrano:2006 ; Seoane  is a
typical scheme for controlling driven dynamical systems. Nonlinear oscillators
described by a second-order differential equation with a fixed damping
parameter, a nonlinear potential, and driven by an external periodic forcing
are well-suited for the application of this method.

In such systems, phase control becomes a fundamental technique to analyze and
tame oscillations. It involves defining and tracking the phase of the
oscillator relative to a reference point, which is typically expressed as an
angle or dimensionless quantity. This phase information allows researchers and
engineers to understand the temporal relationships within the oscillator’s
dynamics, making it possible to predict, control, and synchronize its
behavior. Phase control has applications in diverse fields, from physics and
biology to engineering and information technology, offering insights into
phenomena ranging from cardiac rhythms and brain activity to synchronized
lasers and chaos-based secure communications. By harnessing the concept of
phase control, we can better understand the complex, nonlinear dynamics of
oscillatory systems, enabling advancements in science and technology. A
notable example of the application of phase control can be found in the field
of neuroscience. When studying brain rhythms, researchers have employed phase
control techniques to investigate neural synchronization and communication.
For instance, Buzsaki and Draguhn Buzsaki  demonstrated the significance of
phase relationships in the coordination of neural ensembles during cognitive
processes, shedding light on the brain’s dynamic functioning. In the realm of
physics, phase control is paramount for manipulating the synchronization of
lasers, as exemplified by the work of Strogatz and Mirollo Strogatz1 . They
discussed the concept of ”phase-locking” in coupled oscillators, which is
fundamental in laser physics, where the synchronization of laser light waves
is essential for applications in communication and precision measurement.
Furthermore, the study of heart rhythms and cardiac arrhythmias relies heavily
on phase control techniques. Recent research Namasivayam ; Monga  explored the
use of phase control to investigate the onset of dangerous cardiac
arrhythmias, providing insights into potential therapies for heart-related
disorders. The application of this method has also been utilized to control
the spread of epidemics DUARTE:2021 . All this shows that it remains a topic
of active research with promising implications for diverse scientific and
technological advancements.

In general, all these problems have been faced in the context of non-
fractional derivatives. However, recent works show that, in several physical
situations, the damping term can be fractional and it creates an important
influence in the dynamical behavior of the system Coccolo_fr . In this
specific case, a fractional damping term introduce non-integer-order
derivatives to the equations, allowing for the description of systems with
memory effects, anomalous diffusion, or other non-standard behaviors. These
fractional systems are of interest in various scientific fields, including
physics, engineering, and applied mathematics Boroviec ; Dafermos ; Lu ;
Chellaboina ; De ; Elliott ; Horr ; Ding ; Li ; Strogatz .

Here, we study the phase control technique in the nonlinear Helhmoltz
oscillator with fractional damping. In Ref. Seoane , the authors studied and
analyzed the phase control in the non-fractional Helmholtz oscillator. In this
work, the authors show the optimal value of the phase for which most
trajectories keep inside the well and they do not escape. Now, in our current
situation, the dynamics of a Helmholtz oscillator with fractional damping can
exhibit intriguing and rich behavior Coccolo_fr  since the fractional
parameter $\alpha$ can also play the role of a control parameter and therefore
is relevant for the escaping dynamics. The fractional damping term influences
the rate at which energy is dissipated from the system, and the fractional
order determines the memory effect or the long-term correlation in the damping
process. These features of the fractional derivatives can affect the phase
control. The potential of the Helmholtz oscillator is considered a prototype
for transient chaos and escape phenomena. Therefore, we extend what have been
showed in Seoane  and study the effect of the fractional parameter and the
phase control on the escape of the particles. In particular, we analyze the
attractor and the escape time of the particle, i.e., the final state of the
trajectories, inside or outside the potential well, and the time to reach such
state, respectively.

Our objective is to effectively manage the impact of minor disturbances as a
control. In particular, we assume that the strength of our control is smaller
than the external forcing. As a consequence, apart from optimizing energy
consumption for controlling the trajectory, it shows a relevant characteristic
of a control technique.

This paper is organized as follows. In Sec. II, we present the Helmholtz
oscillator with the fractional damping. The phase control scheme and the
consequences of the dynamical behavior of the system are presented in Sec.
III. The effects of the fractional parameter when we vary the phase and the
amplitude of the control forcing are described in Sec. IV. Likewise, Sec. V
shows how the fractional parameter affects the control of escapes from the
potential well for different initial conditions. Finally, the main conclusions
and a discussion of the results are presented in Sec. VI.

##  II Model description

Figure 1: This plot represents the potential of the Helmholtz oscillator, the
effect of the modulation term $[1+\epsilon\cos{(t+\phi)}]$ is to make the
bottom of the well oscillate vertically.

The Helmholtz oscillator is represented by a nonlinear second order
differential equation considering the presence of the potential defined as
$V(x)=-\frac{1}{2}x^{2}-\frac{1}{3}x^{3}$, as shown in Fig. 1. It represents
the equation of motion of a unit mass particle in a cubic potential under the
influence of both a periodic forcing and a dissipative force and is given by
equation:

$\displaystyle\ddot{x}+\mu\dot{x}-x-x^{2}=F\cos(\omega t),$ (1)

where $\mu$ is the damping parameter, $F$ and $\omega$ the forcing amplitude,
and the forcing frequency, respectively, taking all of them positive values.
This represents a simple paradigmatic example of a dynamical system with
escapes. In this system, normally, the damping term, considered like the first
derivative in Eq. (1), has an impact proportional to the constant $\mu$ in the
system. For our numerical simulations we fix $\mu=0.1$.

In what follows, the first derivative of Eq. (1) is replaced by an order
$\alpha$ fractional derivative, i.e, $\dot{x}(t)\rightarrow
d^{\alpha}x(t)/dx^{\alpha}$, for which we have the Helmholtz Oscillator with a
fractional damping term $\alpha$. This reads as follows:

$\displaystyle\ddot{x}+\mu\frac{d^{\alpha}x}{dx^{\alpha}}-x-x^{2}=F\cos(\omega
t).$ (2)

Adding a fractional damping term to the Helmholtz oscillator equation involves
replacing the integer-order damping term with a fractional derivative term.
Mathematically, if $v(t)=\frac{d}{dt}x(t)$ is the velocity, then the
fractional derivative of order $\alpha$ could be represented as
$v^{\alpha}(t)=\frac{d^{\alpha}x(t)}{dt^{\alpha}}$. This fractional derivative
does not have a straightforward physical interpretation like the first-order
derivative does, because it blends the effects of all past states of the
system into the current state. In fact, if the first-order derivative of space
with respect to time represents velocity, then a fractional derivative of
space with respect to time can be thought of as a “fractional velocity,” which
incorporates the memory effect of the velocity. This means that the fractional
velocity at a given moment is influenced not only by the current rate of
change of position but also by the entire history of the object’s motion. In
practical terms, the fractional derivative introduces memory effects or non-
local influences into the damping mechanism. Physically, the fractional
damping term could represent phenomena such as viscoelasticity in materials,
where the damping force depends not only on the current velocity of the
oscillator but also on its past velocities. This means that the velocity has a
history-dependent nature. It could also model systems with complex
interactions between different parts, where the damping effect is influenced
by the history of the motion Suzuki .

To integrate the system, Eq. 2 can be written as a set of three fractional
differential equations as follows:

$\displaystyle\frac{d^{\alpha}x}{dt^{\alpha}}$ $\displaystyle=$ $\displaystyle
y$ (3) $\displaystyle\frac{d^{1-\alpha}y}{dy^{1-\alpha}}$ $\displaystyle=$
$\displaystyle z$ $\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle
F\cos(\omega t)+x^{2}+x-\mu y,$

where $z$ is a mathematical component coming from the transformation of the
model into a fractional order system. To obtain the solution of Eqs. (3), we
use the Grünwald-Letnikov ALGORITHM:1  fractional derivative, for which the
algorithm to numerically solve this system is given by

$\displaystyle x(t_{k})$ $\displaystyle=$ $\displaystyle
y(t_{k-1})h^{\alpha}-\sum_{j=\upsilon}^{k}c_{j}^{(\alpha)}x(t_{k-j})$ (4)
$\displaystyle y(t_{k})$ $\displaystyle=$ $\displaystyle
z(t_{k-1})h^{1-\alpha}-\sum_{j=\upsilon}^{k}c_{j}^{(1-\alpha)}y(t_{k-j})$
$\displaystyle z(t_{k})$ $\displaystyle=$ $\displaystyle\Psi
h-\sum_{j=\upsilon}^{k}c_{j}^{(1)}z(t_{k-j}),$

where $\Psi=F\cos(\omega t_{k})-x^{2}(t_{k})-x(t_{k})-\mu y(t_{k})$ and $h$ is
the discrete-time step. The coefficients $c_{j}^{\alpha}$ are the binomial
coefficients derived in the numerical scheme implemented, $c_{0}^{\alpha}=1$
and

$\displaystyle c_{j}^{\alpha}=(1-\frac{\alpha+1}{j})c_{j-1}^{\alpha}.$ (5)

The model equation that we will analyze once the phase control is introduced
parametrically is the following

$\displaystyle\ddot{x}+\mu\frac{d^{\alpha}x}{dx^{\alpha}}-x-[1+\epsilon\cos{(t+\phi)}]x^{2}=F\cos(\omega
t).$ (6)

The $\epsilon$ parameter is the modulation amplitude and the $\phi$ parameter
is the phase. So, the potential becomes:

$V(x,t)=-\frac{1}{2}x^{2}-\frac{1}{3}[1+\epsilon\cos{(t+\phi)}]x^{3}.$ (7)

If we fix $\epsilon$ and $\phi$ the bottom of the potential well oscillate in
function of $t$ periodically in the vertical direction. This control
implementation was used, in the non-fractional case, in Refs. Zambrano:2006 ;
Seoane  and it was also tested in an electronic circuit which mimics the
dynamics of the system.

Figure 2: Here, we show some trajectories of the Helmholtz oscillator. The
effect of the interaction between parameter $\alpha$ and $\phi$ is evident in
panels (b) and (d). The initial conditions are
$(x_{0},\dot{x}_{0})=(-0.5,-0.1)$ and $F=0.2$. Figure 3: scape parameter set
$(\omega-\alpha)$ for $F=0.2$, $\epsilon=0.1$ and $\phi=\pi$. The color bar
shows the ratio over $1$ of the trajectories that do not leave the potential
well by varying in every point the initial conditions in the region
$[x_{0}]\times[y_{0}]=[-1.5,0]\times[-0.7,0.7]$. We can observe that, in our
case and for our choice of parameters values, in which we have taken resonant
frequencies, $\omega=1$ is an optimal value to avoid escapes as already stated
in Seoane .

##  III Phase control and the influence of the fractional derivative parameter

Here, we implement the phase control scheme in the fractional Helmholtz
oscillator when all initial conditions escape from the potential well. Also,
we work with resonance frequencies between the periodic forcing and the
modulation term. Therefore, we fix $F=0.2$, $\mu=0.1$ and $\omega=1$. This
choice of the frequency value takes into account resonant frequencies between
the main forcing $F\cos\omega t$ and the control term $[1+\epsilon\cos(\omega
t+\phi)]$. The value $\omega=1$ has been identified as optimal for preventing
escapes from the well, as observed in the non-fractional case Seoane .
Nevertheless, we have computed Fig. 3, wherein we demonstrate that, in the
fractional case and for different values of the forcing frequency, $\omega=1$
locally yields the best performance in retaining most of the particles within
the well. Therefore, the equation of motion reads:

$\displaystyle\ddot{x}+0.1\frac{d^{\alpha}x}{dx^{\alpha}}-x-[1+\epsilon\cos{(t+\phi)}]x^{2}=0.2\cos
t.$ (8)

After switching on the phase control, that means setting the parameter
$\epsilon\neq 0$, we plot trajectories for different values of the parameters
$\alpha$ and $\phi$ in Fig. 2. Here, we fix $\epsilon=0.1$ and initial
conditions $(x_{0},\dot{x}_{0})=(-0.5,-0.1)$ and discard the transient
dynamics. As we can see in Fig. 2(a) the particle is controlled for $\phi=\pi$
in the non-fractional case ($\alpha=1$) and in the fractional case
$\alpha=0.73$ as shown in Fig. 2(c). However, when we change to $\phi=0.7\pi$
the non-fractional case trajectory escapes from the potential well in Fig.
2(b), while in the fractional case it is trapped in it (see Fig. 2(d)). So the
fractional parameter plays a major role to influence the escape or the control
of the particle inside the well. For a better understanding of the effect of
the fractional parameter, we show in Fig. 4 an attractive region ($A$) and the
escape times ($T$) distribution as a function of $\alpha$ for different choice
of parameters $F,\epsilon$ and $\phi$. The attractive region, denoted by the
blue line, assumes the value $1$, when the particle is controlled inside the
well and its escape time is infinite, or $0$ when it escapes. The escape
times, denoted by the red line, have been normalized to $1$ instead of the
maximum time of integration $t=125$ to compare them with the attractive region
in the same conditions. We can see in Fig. 4(a-d) that the escape time and the
attractive region follow similar patterns. Indeed, the information that we
extract from the two curves is different. The blue curve $A$ tells us the
value of the parameter $\alpha$ that control the particle inside the well. The
orange curve $T$ indicates the impact of the parameter on the transient
dynamics of the system. In fact, in Fig. 4(e) and Fig. 4(f), there are no
controlled particles ($A=0$ for all $\alpha$ values), but the escape times
change by varying the $\alpha$ parameter. Moreover, our numerical experiments
have shown that the particles can be controlled with a smaller value of
$\epsilon$ when $F=0.16$. This is the reason why we have decided to keep on
our analysis only with the case $F=0.2$ for which the particles always escape
in the non-fractional case ($\alpha=1$) except for $\epsilon\approx 0.1$ and
$\phi\approx\pi$.

Figure 4: The attractive region $A$ (denoted by the blue line) and the escape
times distribution $T$ (red line, it is normalized at the maximum $1$ instead
of $125$ for a better comparison with the attractor) as a function of the
fractional parameter $\alpha$ are shown in the panels (a-f). When the
attractive region is equal to $1$ the particle does not escape, while when it
is $0$ the particle escapes. We can see that depending on the choice of the
parameters $F,\epsilon$ and $\phi$, the value of $\alpha$ plays a relevant
role to control the trajectories or to obtain different escape times. Figure
5: The attractive region, in the left column, and the escape times
distribution, in the right column, in the parameter set $(\phi-\alpha)$, for
$F=0.2$ and $\epsilon=0.06$ (a) and (b), and for $F=0.2$ and $\epsilon=0.1$
(c) and (d). The color code denotes the escape time of every trajectory. The
yellow structures indicate the regions for which the particles spend more time
to escape from the potential well. Figure 6: The attractive region, in the
left column, and the escape times distribution, in the right column, in the
parameter set $(\epsilon-\alpha)$, for $F=0.2$ and $\phi=0.7\pi$ (a) and (b),
and for $F=0.2$ and $\phi=\pi$ (c) and (d). The color code denotes the escape
time of every trajectory. The yellow structures indicate the regions for which
the particles spend more time to escape from the potential well. Figure 7: The
basins of attraction, in the left column, and the escape times distribution,
in the right column, in the initial conditions space $(x_{0},y_{0})$, for
$F=0.2$ and $\phi=\pi$ (a) and (b), and for $F=0.2$ and $\phi=0.7\pi$ (c) and
(d). In all the panels $\epsilon=0.1,\alpha=0.73$ and the initial conditions
are set in the region $[-1.5,0]\times[-0.7,0.7]$. The percentages inside the
figures of the left column represent the ratio of the initial conditions for
which the trajectories do not escape from the potential well. Figure 8: Escape
parameter set $(\phi-\alpha)$, for $F=0.2$ and $\epsilon=0.06$ (a) and
$\epsilon=0.1$ (b). In all panels, the color bar shows the ratio over $1$ of
the trajectories that do not leave the potential well by varying in every
point the initial conditions in the region
$[x_{0}]\times[y_{0}]=[-1.5,0]\times[-0.7,0.7]$. Figure 9: Escape parameter
set $(\epsilon-\alpha)$, for $F=0.2$ and $\phi=0.7\pi$ (a) and $\phi=\pi$ (b).
In all panels, the color bar shows the ratio over $1$ of the trajectories that
do not leave the potential well by varying in every point the initial
conditions in the region $[x_{0}]\times[y_{0}]=[-1.5,0]\times[-0.7,0.7]$.
Figure 10: The optimal $\phi$ as a function of $\alpha$ (blue curve), for
$F=0.2$ and $\epsilon=0.06$ (a) and $\epsilon=0.1$ (b). The panels show the
$\phi$ value for which the ratio of the controlled trajectories is maximum for
the initial conditions in the region
$[x_{0}]\times[y_{0}]=[-1.5,0]\times[-0.7,0.7]$. The red solid line shows the
value of $\pi$. The red curve, beneath the blue one, represents the percentage
normalized to $1$ of the controlled trajectory for the represented
$\phi_{OPT}$ along the $\alpha$ axis.

##  IV Effect of the fractional parameter when the phase control parameters
change

The effect of the fractional parameter in function of the phase $\phi$ is
evident in Fig. 5. Here, we show in Fig. 5(a) and Fig. 5(b) the attractive
region (the yellow points) and the escape times for $F=0.2$ and
$\epsilon=0.06$. Then, we show in Fig. 5(c) and 5(d) the attractive region and
the escape times for $F=0.2$ and $\epsilon=0.1$ as a function of the phase
$\phi$ and the fractional parameter $\alpha$. We have denoted a color code bar
to indicate the escape time of every trajectory as a function of the
corresponding values of $\alpha$ and $\phi$. It is important to remark that
the attractive regions give us at first glance the parameter values for which
the particles remain bounded inside the well. On the other hand, the escape
times show the parameters values for which we can modify the system dynamics
to obtain longer or smaller time for the trajectory to escape. Of course, both
are related since the region of maximum time delay are the yellow attractive
region. Therefore, it can be useful to see them separately. We have
experienced along all the numerical simulations that, generally, the higher
the intensity of the parameter $\epsilon$, the larger the area of the
attractive region. In fact, in Fig. 5(a) the controlled trajectories are
$3.5\%$ of the total, while in Fig. 5(c) we can control the $7.8\%$. Also, the
escape times are strongly related with the attractive region topology. Similar
behaviors are reported in Fig. 6 where the attractive region and the escape
times gradient plot are depicted in the parameter set $(\epsilon-\alpha)$. The
range of $\epsilon$ values has been selected to guarantee that the modulation
term is smaller than the amplitude of the external periodic forcing. This
ensures that the control is obtained through minimizing the impact of
perturbations on the system and optimizing the energy used to bound the
trajectories. We mean $\epsilon\in(0,0.1]$ which is smaller than $F=0.2$. So,
in Fig. 6(a) and in Fig. 6(b) we have set $\phi=0.7\pi$ and in Fig. 6(c) and
Fig. 6(d) $\phi=\pi$. The first $\phi$ value has been chosen because it is a
value of $\phi$ that shows lots of escaping trajectories in Fig. 5, but also
intersects briefly the attractive region. The second one is the value of
$\phi$ that cuts the attractive region in the middle. The difference in
controlled trajectories is evident. In fact, in Fig. 6(a) they are the $3.5\%$
of the total, while in Fig. 6(c) they are $13.8\%$.

The proposed figures express the deep impact of the $\alpha$ parameter on the
phase control method. In fact, on the border of the attractive region a small
change of the parameter can be pivotal to decide if the trajectory is
controlled inside the well or not.

##  V Control of escapes on the initial conditions

The importance of the fractional parameter on the phase control is well
described by the figures presented until now. Nevertheless, to deepen the
understanding of its impact, we decided to carry on the study by changing the
initial conditions in order to complete the previous study carried out in the
last sections. Therefore, we analyze the basins of attraction in the space of
initial conditions. We define a basin of attraction as the set of initial
conditions that lead to a certain attractor or fixed point. In our case, our
attractors are the infinity (when the particle escapes) and the region inside
the potential well (when the particles do not escape). For this purpose, we
have depicted in Fig. 7, the basin of attractions in the left column and the
escape time in the right column in the initial conditions space
$[-1.5,0]\times[-0.7,0.7]$. Here, we have fixed $\epsilon=0.1$ and $F=0.2$ in
Fig. 7, where the phase is $\phi=\pi$ in the first two panels and
$\phi=0.7\pi$ in the second two panels. We indicate the percentages of the
number of the controlled trajectories (denoted in yellow) thinking similarly
to the framework called dynamical integrity Thompson ; Lenci . This can be
defined as the analysis that can provide valuable information about the
expected dynamics under realistic conditions. It permits to safely operate a
system with the desired behavior, depending on expected disturbances (for us
the modulation term and the periodic forcing). We can see the impact of the
phase $\phi$ in the dynamic of the system when the initial conditions change.
Therefore, we have decided to calculate how many trajectories can be
controlled if the initial conditions change in the region stated above for
different values of $\phi$ and $\alpha$. In Fig. 8, we plot the Escape
parameter set, that shows the percentage of controlled trajectories with
initial conditions $[x_{0}\times y_{0}]=[-1.5,0]\times[-0.7,0.7]$, defined in
a grid of $40\times 40$, in the parameter set $(\phi-\alpha)$. Here, the
percentage values are normalized to $1$ for a better display of the results.
The parameters values are $F=0.2$ in both panels and $\epsilon=0.06$ in Fig.
8(a) and $\epsilon=0.1$ in Fig. 8(b). It is possible to spot the difference in
the percentage values between the two figures due to the different amplitude
of the modulation term. More important is to stress out the complexity of the
effect of the fractional parameter on the escape ratio of the particles. For a
better understanding of the importance of the fractional parameter, we have
carried out in Fig. 9 a similar analysis in the parameter set
$(\epsilon-\alpha)$. In Fig. 9(a) we have fixed $\phi=0.7\pi$ and in Fig. 9(b)
$\phi=\pi$. It is possible to appreciate the difference in the percentage
value between the two figures and how the effect of the interaction between
$\epsilon$ and $\alpha$ is not simple. In fact, although the numerical
simulations show that, generally, the trend should be the bigger the value of
$\epsilon$ the bigger the percentage of the controlled trajectories. Notice
that locally it is not always respected and we can observe there some
fluctuations. Moreover, Fig. 9(a) shows a different trend and we can see that
fact very clearly.

Finally, from the data of Fig. 8 we extrapolate the optimal $\phi$ values that
control the larger number of trajectories in function of the fractional
parameter. These values are depicted in Fig. 10. Here, it is possible to
appreciate a fluctuation around $\phi=\pi$ in both cases $\epsilon=0.06$ and
$\epsilon=0.1$. In the case of $\epsilon=0.06$, where the control amplitude is
smaller, the average value of the optimal phase is
$\overline{\phi_{op}}=3.0913$, with a standard deviation of $\sigma=0.45$. On
then other hand, for the case of $\epsilon=0.1$, the average value of the
optimal phase is $\overline{\phi_{op}}=3.1479$, with a value of $\sigma=0.21$.
For $\epsilon=0.06$, the control amplitude is not so large and, therefore, the
fluctuations around $\phi=\pi$ are more visible although the general trend is
quite close to $\pi$, as corroborated by the value of the $\sigma$. In the
case of $\epsilon=0.1$, in which the control amplitude is larger than the
previous one, the fluctuations around $\phi=\pi$ are not so visible as in the
previous case since the forcing amplitude is predominant in the effects
between $\epsilon$ and $\alpha$. Again, the value of $\sigma$ confirms the
analysis. It means that, in this last situation, the role of the fractional
parameter $\alpha$ is negligible on the optimal value of $\phi$ and we mainly
recover the non-fractional situation in which $\phi=\pi\simeq\phi_{op}$. This
means that the smaller the $\epsilon$ value the stronger the influence of the
fractional parameter on the optimal value of the phase $\phi$ and consequently
on the control of the escapes. Nevertheless, we notice that the percentage of
controlled trajectories differs significantly for each optimal phase value in
the two cases. As illustrated by the red curves in Fig. 10, when
$\epsilon=0.06$, the maximum number of trajectories that the optimal $\phi$
can control is in percentage less then $50\%$, while it is higher in the other
case. Then, for the sake of clarity, we can integrate the $\phi_{OPT}$ curve
along the $\alpha$ axis and calculate the percentage

$n=\frac{1}{N}\int_{\alpha_{m}}^{\alpha_{M}}\phi_{OPT}(\alpha)d\alpha,$ (9)

where $\alpha_{m}=0.5$, $\alpha_{M}=1.5$ and $N$ is the number of launched
trajectories for each $\alpha$ value. So, in the first case, using the optimal
$\phi$ values, we are able to control the $11.5\%$ of all the calculated
trajectories, while in the second case we can control the $24.7\%$. Moreover,
it is essential to highlight that the red curves vary significantly in
function of the $\alpha$ parameter. This means that the proportion of the
bounded trajectories is directly related to the value of $\alpha$. Therefore,
this study underscores the impact of the fractional parameter on the
percentage of controlled trajectories, even when utilizing the values of
$\phi=\phi_{OPT}$.

To summarize, the main findings of this section are as follows: the phase
control technique demonstrates robustness when varying the fractional
parameter value, with the optimal $\phi$ value consistently oscillating around
$\phi_{OPT}=\pi$. Conversely, altering $\alpha$ significantly affects the
percentage of controlled trajectories when initial conditions vary.

##  VI Conclusions and discussion

By using the Helmholtz oscillator with a fractional damping, we have analyzed
the effects of a parametric control with a phase difference $\phi$ between the
main forcing and the control. Afterwards, we have studied the influence of the
fractional parameter $\alpha$ on the control. In the situation in which all
particles are escaping from the potential well, $F=0.2$, in which $\alpha=1$
(non-fractional case), we have found the optimal values of the phase $\phi$
for different values of the fractional parameter. This optimal value varies
with $\alpha$ but its average is close to the ones corresponding to the non-
fractional case, $\phi\approx\pi$. The fluctuations of the optimal value of
$\phi$ with respect the variation of $\alpha$ are illustrated in detail and
they are due to the decrease or increase or the damping term. Although the
fractional parameter may not significantly influence the main optimal value of
$\phi$, it does play a crucial role in determining the proportion of
controlled trajectories, even when utilizing an optimal value of
$\phi=\phi_{OPT}$. Finally, and since the optimal value of the phase
controlling the escapes oscillates around the value in the non-fractional
case, we conjecture that the phase control technique is a robust way to
control the escapes in an open dynamical systems.

##  VII Acknowledgment

This work has been supported by the Spanish State Research Agency (AEI) and
the European Regional Development Fund (ERDF, EU) under Project No.
PID2019-105554GB-I0 (MCIN/AEI/10.13039/501100011033).

## References

  * (1) Boccaletti S, Grebogi C, Lai YC, Mancini H, Maza D. The control of chaos: theory and applications. Phys. Rep. 2000; 329: 103-197.
  * (2) Zambrano S, Allaria E, Brugioni S, Leyva I, Meucci R, Sanjuán MAF, Arecchi FT. Numerical and experimental exploration of phase control of chaos. Chaos 2006; 16; 013111.
  * (3) Seoane JM, Zambrano S, Euzzor S, Meucci R, Arecchi FT, Sanjuán MAF. Avoiding escapes in open dynamical systems using phase control. Phys. Rev. E 2008; 78: 016205.
  * (4) Buzsaki G, Draguhn A. Neuronal oscillations in cortical networks. Science 2004; 304: 1926-1929.
  * (5) Strogatz SH, Mirollo RE. Phase-locking and critical phenomena in lattices of coupled nonlinear oscillators with random intrinsic frequencies. Physica D 1988; 31: 143-168.
  * (6) Namasivayam M, McCready M, Adji A, O’Rourke MF. The role of heart rate in diastolic coronary perfusion and subclinical myocardial ischemia. J. Am. Coll. Cardiol. 2017; 69: 1647-1647.
  * (7) Monga B, Moehlis J. Optimal phase control of biological oscillators using augmented phase reduction. Biol. Cybern. 2019; 113: 161-178.
  * (8) Duarte J, Januário C, Martins N, Seoane, JM, Sanjuán MAF. Controlling Infectious Diseases: The Decisive Phase Effect on a Seasonal Vaccination Strategy. Int J Bifurc Chaos 2021; 31: 2130044.
  * (9) Ortiz A, Yang J, Coccolo M, Seoane JM, Sanjuán MAF. Fractional damping enhances chaos in the nonlinear Helmholtz oscillator. Nonlinear Dyn. 2020; 102: 2323-2337.
  * (10) Borowiec M, Litak G, Syta A. Vibration of the Duffing oscillator: effect of fractional damping. Shock Vib. 2007; 14: 29-36.
  * (11) Dafermos C. Dissipation in materials with memory. In: Lodge AS, Renardy M, Nohel JA, editors. Viscoelasticity and Rheology. Proceedings of a Symposium Conducted by the Mathematics Research Center, Wisconsin–Madison: Academic Press; 1985, p. 221-234.
  * (12) Lu Z, Wang Z, Zhou Y, Lu X. Nonlinear dissipative devices in structural vibration control: A review. J. Sound Vib. 2018; 423: 18-49.
  * (13) Chellaboina V, Haddad WM. Exponentially dissipative nonlinear dynamical systems: a nonlinear extension of strict positive realness. In: Jeltsema D, Scherpen J, editors. Proceedings of the 2000 American Control Conference, Chicago: IEEE Xplore; 2000, p. 3123-3127.
  * (14) De S, Kunal K, Aluru NR. Nonlinear intrinsic dissipation in single layer MoS 2 resonators. RSC Adv. 2017; 7: 6403-6410.
  * (15) Elliott SJ, Tehrani MG, Langley RS. Nonlinear damping and quasi-linear modelling. Philos. Trans. R. Soc. A 2015; 373: 20140402.
  * (16) Horr AM, Schmidt LC. A fractional-spectral method for vibration of damped space structures. Eng. Struct. 1996; 18: 947-956.
  * (17) Ding Y, Liu X, Chen P, Luo X, Luo Y. Fractional-Order Impedance Control for Robot Manipulator. Fractal fract. 2022; 6: 684.
  * (18) Li X, Liu Z, Li J, Tisdell C. Existence and controllability for nonlinear fractional control systems with damping in Hilbert spaces. Acta Math. Sci. 2019; 39; 229-242.
  * (19) Strogatz SH. Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering. Boca Raton: CRC press; 2018.
  * (20) Suzuki JL, Gulian M, Zayernouri M, D’Elia M. Fractional modeling in action: A survey of nonlocal models for subsurface transport, turbulent flows, and anomalous materials. J. Peridynam. Nonlocal Model. 2023;5(3):392-459.
  * (21) Diethelm, K, Ford, NJ, Freed, AD, Luchko, Y. Algorithms for the fractional calculus: A selection of numerical methods. Comput. Methods Appl. Mech. Engrg. 2005; 194; 743–773.
  * (22) Thompson JMT.: Dynamical integrity: three decades of progress from macro to nanomechanics. In: Global Nonlinear Dynamics for Engineering Design and System Safety, pp. 1–26. Springer (2019)
  * (23) Lenci S, Rega G. Optimal control of homoclinic bifurcation: Theoretical treatment and practical reduction of safe basin erosion in the Helmholtz oscillator. J. Vib. Control 2003; 9; 281–316.