Let $(u_{k+1}, w_{k+1})_{k\ge 0}$ and $(u'_{k+1}, w'_{k+1})_{k\ge 0}$ be two trajectories of the NISPP method (Algorithm <ref>) for any adjacent datasets $\bS \simeq \bS'$ whose batches are denotes by $\bB_{k+1}, \bB'_{k+1}$ respectively. Moreover, denote $\mathbf{a}_{k+1}  :=  \gnorm{F_{\bB_{k+1}}(u_{k+1}) - F_{\bB'_{k+1}}(u_{k+1})}{}{}$ and $\delta_{k+1} := \gnorm{u_{k+1}-u'_{k+1}}{}{} (= \gnorm{w_{k+1}-w'_{k+1}}{}{})$ for  $k$-th iteration of Algorithm <ref>. Then, if $i = \inf\{k : \bB_{k+1} \ne \bB'_{k+1}\}$,

\begin{equation}
\delta_{j+1} \begin{cases}
=0 &\text{if } j+1 \le i\\
\le \sum_{k = i}^j\frac{2\mathbf{a}_{k+1}}{\lambda_k} + \sqrt{\frac{4\nu}{\lambda_k}} &\text{otherwise}.
\end{cases}\label{eq:int_rel_prox_28}
\end{equation}

It is clear from the definition $i$ that $\bB_j = \bB_j'$ for all $j \le i$. This implies $u_j = u'_j$ and $w_j = w_j'$ for all $j \le i$. Hence, we conclude first case of (<ref>).

Using (<ref>) for $\nu$-approximate strong VI solution, we have,

\begin{align}
\inprod{F_{\bB_{k+1}}(u_{k+1}) + \lambda_k(u_{k+1}-w_k)}{w - u_{k+1}} &\ge -\nu\label{eq:int_rel_prox_26},\\
\inprod{F_{\bB'_{k+1}}(u'_{k+1}) + \lambda_k(u'_{k+1}-w'_k)}{w - u'_{k+1}} &\ge -\nu\label{eq:int_rel_prox_27}.
\end{align}
Then, adding (<ref>) with $w = u'_{k+1}$ and (<ref>) with $w = u_{k+1}$, we have

\begin{align}
\inprod{F_{\bbeta_{k+1}}(u_{k+1}) - F_{\bbeta'_{k+1}}(u'_{k+1})}{u_{k+1}-u'_{k+1}} &\le 2\nu - \lambda_k\inprod{u_{k+1} -u'_{k+1}}{(u_{k+1}-w_k) - (u'_{k+1}-w'_k)}\nonumber\\
&=2\nu - \lambda_k \delta_{k+1}^2  + \lambda_k\inprod{u_{k+1} -u'_{k+1}}{w_k - w_k'} \nonumber\\
&\le 2\nu - \lambda_k\delta_{k+1}^2 + \frac{\lambda_k}{2}[\delta_{k}^2 + \delta_{k+1}^2]\nonumber\\
&\le 2\nu  -\frac{\lambda_k}{2}\delta_{k+1}^2 + \frac{\lambda_k}{2}\delta_{k}^2 \label{eq:int_rel_prox_25}.
\end{align}

Also note that

\begin{align*}
\inprod{F_{\bB_{k+1}}(u_{k+1}) - F_{\bB'_{k+1}}(u'_{k+1})}{u_{k+1}-u'_{k+1}} &= \inprod{F_{\bB'_{k+1}}(u_{k+1}) - F_{\bB'_{k+1}}(u'_{k+1})}{u_{k+1}-u'_{k+1}} \\
&\quad+ \inprod{F_{\bB_{k+1}}(u_{k+1}) - F_{\bB'_{k+1}}(u_{k+1})}{u_{k+1}-u'_{k+1}}\\
&\ge \inprod{F_{\bB_{k+1}}(u_{k+1}) - F_{\bB'_{k+1}}(u_{k+1})}{u_{k+1}-u'_{k+1}}
\end{align*}

where the last inequality follows from monotonicity of $F_{\bB'_{k+1}}$. Using above relation along with (<ref>), we obtain

\begin{align*}
\frac{\lambda_k}{2} \delta_{k+1}^2 &\le \frac{\lambda_k}{2}\delta_{k}^2 + 2\nu + \inprod{F_{\bB'_{k+1}}(u_{k+1}) - F_{\bB_{k+1}}(u_{k+1})}{u_{k+1}-u'_{k+1}}\\
\Rightarrow \delta_{k+1}^2 &\le \delta_{k}^2 + \frac{4\nu}{\lambda_k} + \frac{2}{\lambda_k}\mathbf{a}_{k+1}\delta_{k+1},
\end{align*}
where we used the definition $a_k$ along with Cauchy-Schwarz inequality. Solving for the quadratic inequality in $\delta_{k+1}$, we obtain the following recursion

\[ \delta_{k+1} \le \frac{\mathbf{a}_{k+1}}{\lambda_k} + \sqrt{\frac{\mathbf{a}_{k+1}^2}{\lambda_k^2} + \delta_{k}^2 + \frac{4\nu}{\lambda_k}}\]

which can be further simplified to

\[ \delta_{k+1} \le \delta_{k} + \frac{2\mathbf{a}_{k+1}}{\lambda_k} + \sqrt{\frac{4\nu}{\lambda_k}}.\]

Solving this recursion and noting the base case that $\delta_i = 0$, we obtain (<ref>).

A direct consequence of the previous analysis are in-expectation and high probability uniform argument stability upper bounds
for the sampling with replacement variant of Algorithm <ref>.

Let $\A$ denote the sampling with replacement NISPP method (Algorithm <ref>) where $\bB_k$ is chosen uniformly at random from subsets of $\bS$ of a given size $B_k$. Then $\A$
satisfies the following uniform argument stability bounds:

$$ \sup_{\bS\simeq \bS^{\prime}}\EE_{\A}[\delta_{\A}(\bS,\bS^{\prime})]
\leq  \sum_{k=1}^{K}\Big( \frac{2M}{n\lambda_k}+\sqrt{\frac{4\nu}{\lambda_k}} \Big) .$$

Furthermore, if $|\bB_{k}|= B$ and $\lambda_k= \lambda$ for all $k$ (i.e., constant batch size and regularization parameter throughout iterations) then w.p. at least $1-n\exp\{-KB/[4n]\}$ (over both sampling and noise addition)

$$ \sup_{\bS\simeq \bS^{\prime}}[\delta_{\A}(\bS,\bS^{\prime})] \leq \frac{4MK}{\lambda n}+K\sqrt{\frac{4\nu}{\lambda}}.$$

Let $\bS\simeq \bS^{\prime}$ and $(u_{k+1})_k$, $(u_{k+1}^{\prime})_k$ the trajectories of
the algorithm on $\bS$ and $\bS^{\prime}$, respectively.
By Lemma <ref>, letting $\delta_{k+1}=\|\tilde w_{k+1}-\tilde w_{k+1}^{\prime}\|$, we get
$\delta_{j} \leq \sum_{k=1}^{j}\Big( \frac{2\mathbf{a}_k}{\lambda_k}+\sqrt{\frac{4\nu}{\lambda_k}} \Big),$
where $\mathbf{a}_k=\gnorm{F_{\bB_{k+1}}(u_{k+1})-F_{\bB_{k+1}'}(u_{k+1}')}{}{}$ is a random
variable. By the law of total probability, $\EE[\mathbf{a}_k]
\leq \frac{|\bB_{k+1}|}{n}\frac{2M}{|\bB_{k+1}|}+\big(1-\frac{|\bB_{k+1}|}{n}\big)\cdot 0=\frac{2M}{n}.$
Hence, $\EE[\delta_{j}]\leq \sum_{k=1}^{j}\Big( \frac{2M}{n\lambda_k}+\sqrt{\frac{4\nu}{\lambda_k}} \Big)
\leq \sum_{k=1}^{K}\Big( \frac{2M}{n\lambda_k}+\sqrt{\frac{4\nu}{\lambda_k}} \Big).$
Since $ \gnorm{\proj_{\W}(\wb{w}_K)-\proj_{\W}(\wb{w}_K^{\prime})}{}{} \le \gnorm{\wb{w}_K - \wb{w}'_K}{}{} \le \max_{k\in[K]}\delta_k$, and since $\bS\simeq \bS^{\prime}$ are arbitrary,

$$ \sup_{\bS\simeq \bS^{\prime}}\EE_{\A}[\delta_{\cal A}(\bS,\bS^{\prime})]
\leq \sum_{k=1}^{K}\Big( \frac{2M}{n\lambda_k}+\sqrt{\frac{4\nu}{\lambda_k}} \Big) .$$

We proceed now to the high-probability bound. Let $\br_k\sim\mbox{Ber}(p)$, for $k\in[K]$,
with $Kp< 1$. Then, for any $0<\theta<1/2$,

\begin{eqnarray*}
\PP\left[  \sum_{k=1}^K\br_k \geq Kp+\tau \right]
&\leq& \exp\big(-\theta(\tau+Kp))\Big[ 1+p(e^{\theta}-1) \Big]^K \leq \exp\{Kp\theta^2-\theta\tau\}.
\end{eqnarray*}

Choosing $\theta=\tau/[2Kp]<1/2$, we get that the probability above is upper bounded by $\exp\{-\tau^2/[4Kp]\}$.
Finally, choosing $\tau=Kp$, we get

$$\PP\Big[  \sum_{k=1}^K\br_k \geq 2Kp \Big] \leq \exp\{-Kp/4\}. $$

Next, fix the coordinate $i$ where $S\simeq S^{\prime}$ may differ. Noticing that $\mathbf{a}_k$
is a.s. upper bounded by $(2M/B)\br_k$ with $\br_k\sim\mbox{Ber}(p)$, with $p=B/n$, we get

$$ \PP\Big[ \sum_{k=1}^K\frac{2\mathbf{a}_k}{\lambda} \geq \frac{2}{\lambda}\frac{2M}{n}\Big] \leq \exp\{-\frac{KB}{4n}\}. $$

In particular, w.p. at least $1-\exp\{-\frac{KB}{4n}\}$, we have
$ \mathbf{a}_k \leq \frac{4M}{\lambda n}+\sqrt{\frac{4\nu}{\lambda}}.$
Using a union bound over $i\in[n]$ (and noticing that averaging preserves the stability bound), we conclude that

$$ \PP\Big[ \sup_{\bS\simeq \bS^{\prime}} \delta_{\A}(\bS,\bS^{\prime}) >\frac{4MK}{\lambda n}+K\sqrt{\frac{4\nu}{\lambda}}\Big] \leq n\exp\{-KB/4n\}. $$

Hence, we conclude the proof.

An important observation, for the high-probability guarantee to be useful, is necessary that the algorithm is run for sufficiently many iterations; in particular, we require $K=\omega(n/B)$. Whether this assumption can be completely avoided is an interesting question. Nevertheless, as we will see in the section, our policy for DP-SVI and DP-SSP problem satisfies this requirement.

The following corollary is in order.

For Algorithm <ref>, we have

\[\EE[\GenP{x(\A; \bS)}{\bS} +\GenD{y(\A; \bS)}{\bS}] \le 2M \tsum_{k =1}^K \paran*{\frac{2M}{n\lambda_k} + \sqrt{\frac{4\nu}{\lambda_k}}}.\]

Proof follows directly from the fact $w = [x,y]^T$ is stable then $\gnorm{x(\bS)- x(\bS')}{}{} + \gnorm{y(\bS) - y(\bS')}{}{} \le 2\gnorm{w(\bS)-w(\bS')}{}{}$. Noting this relation along with bound on expected stability of $w$ (provided in Proposition <ref>) and Lemma <ref>, we conclude the proof.

§.§ Optimal risk for DP-SVI and DP-SSP by the NISPP method

In previous three sections, we provided bounds on optimization error, generalization error and value of $\sigma$ for obtaining $(\priva,\privb)$-differential privacy. In this section, we specify a policy for selecting $\lambda_k, B_k, \gamma_k, \sigma_k $ and $\nu$ such that requirement in previous three sections are satisfied and we can get optimal risk bounds while maintaining $(\priva,\privb)$-privacy. In particular, consider the multi-pass NISPP method where each sample batch $\bB_k$ is chosen randomly from subsets of $\bS$ with replacement. Then, we have the following theorem:

Let $\A$ be the multi-pass NISPP method (Algorithm <ref>). Set the following constant stepsize and batchsize policy for $\A$:

γ_k = 1,           λ_k = λ_0= max*M/D,Lmax*√(n), √(dlog1/δ)/θ,           B_k = B= √(n),

K = n,        ν= M^2/λ_0n^2,         σ_k+1 = 8M/B λ_0√(ln(1/η))/ .

Then, Algorithm <ref> is $(\priva, \privb)$-differential private. Moreover, output $\A(\bS)$ satisfies the following bound on $\EE_{{\cal A}}[\wVIgap{\A(\bS)}{F}]$ for SVI problem (or $\EE_{{\cal A}}[\wSPgap{\A(\bS)}{f}]$ for SSP problem)

\[ O\paran*{(M+LD)D \paran*{\frac{1}{\sqrt{n}} + \frac{\sqrt{d\ln{1/\eta}}}{\eps n}}},\]

Moreover, such solution is obtained in total of $\wt{O}(n\sqrt{n})$ stochastic operator evaluations.

Note that since $\nu$ satisfies assumption in Proposition <ref>, we have $\ell_2$-sensitivity of the update of $u_{k+1}$ is $\frac{4M}{\lambda_0B_{k+1}}$. Then, in view of Theorem <ref> along with value of $\sigma_{k+1}$, we conclude that Algorithm <ref> is $(\priva, \privb)$-differential private.

Now, for convergence, we first bound the empirical gap. Given that our bounds for (<ref>) and (<ref>) are analogous, we proceed indistinctively for both problems. By Theorem <ref>, along with the fact that sampling with replacement is an unbiased stochastic oracle for the empirical operator, we have for any $\bS$

\begin{align*}
\EE_{\A}[\mbox{EmpGap}(\A,F_\bS)] &\le \nu + \frac{\lambda_0D^2}{n} + \frac{7M^2}{\lambda_0} + \frac{160M^2d}{\eps^2B^2\lambda_0}\ln{\frac{1}{\eta}} + M\sqrt{2d}\frac{8M}{Bn \eps \lambda_0}\\
&\le O\paran*{MD \paran*{\frac{1}{\sqrt{n}} + \frac{\sqrt{d \ln{1/\eta}}}{\eps n} 	} }.
%\frac{MD}{n^{5/2}} + \frac{MD}{\sqrt{n}} + \frac{MD\sqrt{d\log{1/\delta}}}{\theta n} +
\end{align*}

Similar claims can be made for empirical gap for (<ref>) problem: $\EE_{{\cal A}}\big[ \mbox{EmpGap}(\A, f_{\bS})\big]$ where output of $\A$ is $(x(\bS), y(\bS))$.

Next, by Theorem <ref>, we have that $\A(\bS)$ (or $x(\bS)$ and $y(\bS)$ for the SSP case) are UAS with parameter

\[\delta = \frac{2M}{\lambda_0} + n\sqrt{\frac{4\nu}{\lambda_0}} = \frac{4M}{\lambda_0} \le \frac{4D}{\sqrt{n}}. \]

Hence, noting that empirical risk upper bounds weak VI or SP gap, i.e., using Proposition <ref> or Proposition <ref> (depending on whether
the problem is an SSP or SVI, respectively), we have that the risk is upper bounded by its empirical risk plus $M\delta$, where
$\delta$ is the UAS parameter of the algorithm; in particular, if $\mbox{WeakGap}(\A;\bS)$ is the (SVI or SSP, respectively) gap function for the
expectation objective, then

\begin{align*}
\EE_{\A, \bS } \big[ \wVIgap{\A}{F}\big] \le O\Big(MD\Big(\frac{1}{\sqrt{n}} + \frac{\sqrt{d\ln(1/\eta)} } {\eps n}\Big)\Big) + \frac{20MD}{\sqrt{n}} = O\Big(MD\Big(\frac{1}{\sqrt{n}} + \frac{\sqrt{d\ln(1/\eta)}}{\eps n}\Big)\Big)
\end{align*}

Similar claim can be made for $\wSPgap{\A}{f}$.

Finally, we analyze the running time performance. As in Remark <ref> , number of OE method iterations for obtaining $\nu$-approximate solution is $O\big(\frac{L+\lambda_0}{\lambda_0}\ln\big(\frac{LD^2+MD+\lambda_0D^2}{\nu}\big)\big)$. Now note that $\frac{L+\lambda_0}{\lambda_0} \le \frac{\sqrt{n}+1}{\sqrt{n}} \le 2$ since $n \ge 1$. Moreover, in view of (<ref>), we have  $\ln\big(\frac{LD^2+MD+\lambda_0D^2}{\nu}\big) \le \ln\paran[\big]{\frac{4\lambda_0D^2}{\nu}} \le \ln\paran[\big]{n^2\max\braces[\big]{1, \frac{L^2D^2}{M^2}}\max\braces[\big]{n, \frac{d\ln(1/\privb)}{\priva^2}}}$.
Each iteration of OE method costs $B$ stochastic operator evaluations and we run outer-loop of NISPP method $K$ times. Hence, total stochastic operator evaluations (after ignoring the $\ln$-term) $\wt{O}(KB) = \wt{O}(n\sqrt{n})$.
Hence, we conclude the proof.

§ LOWER BOUNDS AND OPTIMALITY OF OUR ALGORITHMS

In this Section, we show the optimality of our obtained rates from Sections <ref> and <ref>. The first observation is that, since DP-SCO corresponds to a DP-SSP problem where ${\cal Y}$ is a singleton, the complexity of DP-SSP is lower bounded
by $\Omega(MD \big( \frac{1}{\sqrt n}+\min\big\{1,\frac{\sqrt{d}}{\varepsilon n}\big\} \big)\big)$: this is a known lower bound for DP-SCO [4]. It is important to note as well that this reduction applies to the weak generalization gap, as defined in (<ref>), as in the case ${\cal Y}=\{\bar{y}\}$ is a singleton:

\begin{align*}
\wSPgap{{\cal A}}{f}
&=\mathbb{E}_{\cal A}[ \sup_{y\in {\cal Y}} \mathbb{E}_{\bS}[f(x(\bS),y)]-\inf_{x\in {\cal X}} \mathbb{E}_{\bS}[f(x, y(\bS))] ]\\
&=  \mathbb{E}_{\cal A}\mathbb{E}_{\bS}[f(x(\bS),\bar{y})]-\inf_{x\in {\cal X}} f(x, \bar{y})\\
&=  \mathbb{E}_{{\cal A},\bS}[f(x(\bS),\bar{y})-\inf_{x\in {\cal X}} f(x, \bar{y})],
\end{align*}

which is simply the expected optimality gap. Using this reduction, together with a lower bound for DP-SCO [4], we conclude that

Let $n,d\in\NN$, $L,M,D,\varepsilon>0$ and $\delta=o(1/n)$. The class of problems DP-SSP with gradient operators within the class ${\cal M}_{\W}^1(L,M)$, and domain ${\cal W}$ containing an Euclidean ball of diameter $D/2$, satisfies the lower bound

$$ \Omega\Big(MD\Big( \frac{1}{\sqrt n}+\min\Big\{1,\frac{\sqrt{d}}{\varepsilon n}\Big\} \Big)\Big).$$

Next, we study the case of DP-SVI. The situation is more subtle here. Our approach is to first prove a reduction from population weak VI gap to empirical strong VI gap for the case where operators are constant w.r.t. $w$. In fact, it seems unlikely this reduction works for more general monotone operators, however this suffices for our purposes, as we will later prove a lower bound construction used for DP-ERM [5] leads to a lower bound for strong VI gap with constant  operators.

The formal reduction to the empirical version of the problem is presented in the following lemma. Its proof follows closely the reduction from DP-SCO to DP-ERM from [4]. Below, given a dataset $\bS\in {\cal Z}^n$, let ${\cal P}_{\bS}=\frac1n\sum_{\bbeta\in\bS}\delta_{\bbeta}$ be the empirical distribution associated with $\bS$.

Let ${\cal A}$ be an $(\priva/[4\log(1/\privb)], e^{-\priva}\privb/[8\log(1/\privb)])$-DP algorithm for SVI problems. Then there exists an $(\priva,\privb)$-DP algorithm ${\cal B}$ such that for any empirical VI problem with constant operators,

\[\EmpVIgap{{\cal B}}{F_{\bS}} \leq \wVIgap{{\cal A}}{F_{{\cal P}_{\bS}}} \qquad (\forall \bS\in {\cal Z}^n).\]

Consider the algorithm ${\cal B}$ that does the following: first, it extracts a sample $\mathbf{T}\sim ({\cal P}_{\bS})^n$; next, executes ${\cal A}$ on $\mathbf{T}$; and finally, outputs ${\cal A}(\mathbf{T})$. We claim that this algorithm is $(\priva,\privb)$-DP w.r.t. $\bS$, which follows easily by bounding the number of repeated examples with high probability, together with the group privacy property applied to ${\cal A}$ (for a more detailed proof, see Appendix C in [4]). Now, given a constant operator $F_{\bbeta}(w)$, let $R(\bbeta)\in\mathbb{R}^d$ be its unique evaluation. Let now $R_{\bS}$ be the unique evaluation of $F_{\bS}$, and given a distribution ${\cal P}$ let $R_{\cal P}$ be the unique evaluation of $F_{{\cal P}}(w)=\mathbb{E}_{\bbeta\sim {\cal P}}[F_{\bbeta}(w)]$.

Noting that $\mathbb{E}_{\mathbf{T}}[R_{\mathbf{T}}]=R_{\bS},$ we have that

\begin{align*}
\EmpVIgap{{\cal B}(\bS)}{F_{\bS}}
&= \mathbb{E}_{\cal B}[\sup_{w\in {\cal W}} \langle R_{\bS},{\cal B}(\bS)-w \rangle]\\
&= \mathbb{E}_{{\cal A},\mathbf{T}}[ \langle R_{\bS},{\cal A}(\mathbf{T})\rangle- \inf_{w\in {\cal W}} \langle R_{\bS},w \rangle]\\
&=\mathbb{E}_{{\cal A}} \sup_{w\in {\cal W}} \mathbb{E}_{\mathbf{T}}[ \langle R_{\bS},{\cal A}(\mathbf{T})-w \rangle]\\
&= \wVIgap{{\cal A}}{F_{{\cal P}_{\bS}}},
\end{align*}

where third equality holds since the optimal choice of $w$ is independent of $\mathbf{T}$, and the
last equality holds by definition of the weak gap function and the fact that $\mathbf{T}\sim ({\cal P}_{\bS})^n$.

In this Section, we show the optimality of our obtained rates from Sections <ref> and <ref>.

The first straightforward observation is that, since DP-SCO reduces to DP-SSP, then the complexity of DP-SSP is lower bounded
by $\Omega(MD \big( \frac{1}{\sqrt n}+\min\big\{1,\frac{\sqrt{d}}{\varepsilon n}\big\} \big)\big)$: this is a known
lower bound for DP-SCO [4]. Now we establish a lower bound on the risk of DP-SVI, which is also
obtained from a DP-SCO problem.

Next, we prove a lower bound for the empirical VI problem over constant operators.

Let $n,d\in\NN$, $L,M,D,\varepsilon>0$ and $2^{-o(n)}\leq \delta \leq o(1/n)$. The class of DP empirical VI problems with constant operators within the class ${\cal M}_{\W}^1(L,M)$,
and domain ${\cal W}$ containing an Euclidean ball of diameter $D/2$ satisfies the lower bound

$$ \Omega\Big(MD\Big( \min\Big\{1,\frac{\sqrt{d\log(1/\privb)}}{\varepsilon n}\Big\} \Big)\Big).$$

Consider the following empirical VI problem: $F_{\bbeta}(u)= M\bbeta$, $\W={\cal B}(0,D)$ and dataset $\bS$ with points contained in $\{-1/\sqrt d,+1/\sqrt d\}^d$. Notice that, since the operator in this case is constant, the VI gap coincides with the excess risk of the associated convex optimization problem

$$ (P)~~~ \min_{u\in \W} \Big\langle\frac{M}{n}\sum_{i\in[n]} \bbeta_i, u\Big\rangle. $$

Indeed, for any $u\in {\cal W}$,

\begin{eqnarray*}
\EmpVIgap{u}{F_\bS} &=& \sup_{v\in {\cal B}(0,D)} \Big\langle \frac{M}{n} \sum_{i\in[n]} \bbeta_i, u- v\Big\rangle
= \Big\langle \frac{M}{n} \sum_{i\in[n]} \bbeta_i, u+ \frac{D\sum_i \bbeta_i}{\|\sum_i \bbeta_i\|}\Big\rangle \\
&=& \Big\|\frac{MD}{n}\sum_{i\in[n]} \bbeta_i\Big\|+MD\big\langle \frac{u}{D}, \frac{1}{n}\sum_i \bbeta_i\big\rangle.
%\Big\langle u, \frac{MD\sum_i \bbeta_i}{\|\sum_i \bbeta_i\|}\Big\rangle.
\end{eqnarray*}

This, together with the lower bounds on excess risk proved for this problem in <cit.> and <cit.>
show that any $(\priva,\privb)$-DP algorithm for this problem must incur in worst-case VI gap $\Omega(MD\min\{1,\frac{\sqrt{d\log(1/\privb)}}{\priva n}\})$,
which proves the result.

The two results above provide the claimed lower bound for the weak SVI gap of any differentially private algorithm.

Let $n,d\in\NN$, $L,M,D,\varepsilon>0$ and $2^{-o(n)}\leq \delta \leq o(1/n)$. The class of DP-SVI problems with
operators within the class ${\cal M}_{\W}^1(L,M)$,
and domain ${\cal W}$ containing an Euclidean ball of diameter $D/2$ satisfies a lower bound for the weak VI gap

$$ \tilde\Omega\Big(MD\Big( \frac{1}{\sqrt n}+\min\Big\{1,\frac{\sqrt{d}}{\varepsilon n}\Big\} \Big)\Big).$$

Before we prove the result, let us observe that the presented lower bound shows the optimality of our algorithms in the
range where $M\geq LD$. Obtaining a matching lower bound for any choice of $M,L,D$ is an interesting question, which
unfortunately our proof technique does not address: this is a limitation that the lower bound is based on constant operators, and therefore their Lipschitz constants are always zero.

Let ${\cal A}$ be any algorithm for SVI. By the classical (nonprivate) lower bounds for SVI [30, 22], we have that the minimax SVI gap attainable is lower bounded by $\Omega(MD/\sqrt n)$. On the other hand, using Lemma <ref> the accuracy of any $(\priva,\privb)$-DP algorithm for weak SVI gap is lower bounded by the strong gap achieved by $(4\priva \ln(1/\privb),e^{\priva}\tilde O(\privb))$-DP algorithms on empirical VI problems with constant operators. Finally, by Proposition <ref>, the latter class of problems enjoys a lower bound $\Omega(\min\{1,\sqrt{d\ln(1/[e^{\priva}\tilde O(\privb)])}/[\varepsilon  n\ln(1/\privb)]\})=\tilde\Omega(\min\{1,\sqrt{d}/[n\varepsilon])$, which implies a lower bound on the former class of this order. We conclude by combining both the private and nonprivate
lower bounds established above.

§ ACKNOWLEDGEMENTS

CG would like to thank Roberto Cominetti for valuable discussions on stochastic variational inequalities and nonexpansive iterations. CG's research is partially supported by INRIA through the INRIA Associate Teams project, ANID – Millenium Science Initiative Program – NCN17_059, and FONDECYT 1210362 project. We would also like to thank the authors of [40] for pointing out a gap in a proof of an earlier version of this paper.

[1]
Martin Abadi, Andy Chu, Ian Goodfellow, H Brendan McMahan, Ilya Mironov, Kunal
Talwar, and Li Zhang.
Deep learning with differential privacy.
In Proceedings of the 2016 ACM SIGSAC Conference on Computer and
Communications Security, pages 308–318. ACM, 2016.

[2]
Alfred Auslender and Marc Teboulle.
Interior projection-like methods for monotone variational
Mathematical Programming, 104(1):39–68, 2005.

[3]
Raef Bassily, Vitaly Feldman, Cristóbal Guzmán, and Kunal Talwar.
Stability of stochastic gradient descent on nonsmooth convex losses.
In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin,
editors, Advances in Neural Information Processing Systems, volume 33,
pages 4381–4391. Curran Associates, Inc., 2020.

[4]
Raef Bassily, Vitaly Feldman, Kunal Talwar, and Abhradeep Guha Thakurta.
Private stochastic convex optimization with optimal rates.
In Advances in Neural Information Processing Systems,
volume 32. Curran Associates, Inc., 2019.

[5]
Raef Bassily, Adam Smith, and Abhradeep Thakurta.
Differentially private empirical risk minimization: Efficient
algorithms and tight error bounds, 2014.

[6]
Olivier Bousquet and André Elisseeff.
Stability and generalization.
J. Mach. Learn. Res., 2:499?526, March 2002.

[7]
Olivier Bousquet, Yegor Klochkov, and Nikita Zhivotovskiy.
Sharper bounds for uniformly stable algorithms.
In Jacob Abernethy and Shivani Agarwal, editors, Proceedings of
Thirty Third Conference on Learning Theory, volume 125 of Proceedings
of Machine Learning Research, pages 610–626. PMLR, 09–12 Jul 2020.

[8]
Kamalika Chaudhuri and Claire Monteleoni.
Privacy-preserving logistic regression.
In NIPS, 2008.

[9]
Kamalika Chaudhuri, Claire Monteleoni, and Anand D Sarwate.
Differentially private empirical risk minimization.
Journal of Machine Learning Research, 12(Mar):1069–1109, 2011.

[10]
Cynthia Dwork and Aaron Roth.
The algorithmic foundations of differential privacy.
Found. Trends Theor. Comput. Sci., 9(3-4):211–407, 2014.

[11]
Cynthia Dwork, Guy N. Rothblum, and Salil P. Vadhan.
Boosting and differential privacy.
In 51th Annual IEEE Symposium on Foundations of Computer
Science, FOCS 2010, October 23-26, 2010, Las Vegas, Nevada, USA, pages
51–60. IEEE Computer Society, 2010.

[12]
Francisco Facchinei and Jong-Shi Pang.
Finite-dimensional variational inequalities and complementarity
Springer Science & Business Media, 2007.

[13]
Vitaly Feldman, Tomer Koren, and Kunal Talwar.
Private stochastic convex optimization: Optimal rates in linear time.
In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory
of Computing, STOC 2020, page 439–449, New York, NY, USA, 2020.
Association for Computing Machinery.

[14]
Vitaly Feldman and Jan Vondrak.
High probability generalization bounds for uniformly stable
algorithms with nearly optimal rate.
In Alina Beygelzimer and Daniel Hsu, editors, Proceedings of the
Thirty-Second Conference on Learning Theory, volume 99 of Proceedings
of Machine Learning Research, pages 1270–1279, Phoenix, USA, 25–28 Jun
2019. PMLR.

[15]
Moritz Hardt, Ben Recht, and Yoram Singer.
Train faster, generalize better: Stability of stochastic gradient
In Maria Florina Balcan and Kilian Q. Weinberger, editors, Proceedings of The 33rd International Conference on Machine Learning,
volume 48 of Proceedings of Machine Learning Research, pages
1225–1234, New York, New York, USA, 20–22 Jun 2016. PMLR.

[16]
Yu-Guan Hsieh, Franck Iutzeler, Jérôme Malick, and Panayotis
On the convergence of single-call stochastic extra-gradient methods.
In Hanna M. Wallach, Hugo Larochelle, Alina Beygelzimer, Florence
d'Alché-Buc, Emily B. Fox, and Roman Garnett, editors, Advances
in Neural Information Processing Systems 32: Annual Conference on Neural
Information Processing Systems 2019, NeurIPS 2019, 8-14 December 2019,
Vancouver, BC, Canada, pages 6936–6946, 2019.

[17]
Alfredo N. Iusem, Alejandro Jofré, Roberto I. Oliveira, and Philip
Variance-based extragradient methods with line search for stochastic
variational inequalities.
SIAM J. Optim., 29(1):175–206, 2019.

[18]
Prateek Jain, Pravesh Kothari, and Abhradeep Thakurta.
Differentially private online learning.
In 25th Annual Conference on Learning Theory (COLT), pages
24.1–24.34, 2012.

[19]
Prateek Jain and Abhradeep Thakurta.
(near) dimension independent risk bounds for differentially private
In ICML, 2014.

[20]
A. Jalilzadeh and U. V. Shanbhag.
A proximal-point algorithm with variable sample-sizes (ppawss) for
monotone stochastic variational inequality problems.
In 2019 Winter Simulation Conference (WSC), pages 3551–3562,

[21]
Anatoli Juditsky and Arkadi Nemirovski.
First order methods for nonsmooth convex large-scale optimization, i:
General purpose methods, 2012.
MIT Press.

[22]
Anatoli Juditsky, Arkadii S. Nemirovski, and Claire Tauvel.
Solving variational inequalities with Stochastic Mirror-Prox
Stochastic Systems, 1(1):17–58, 2011.

[23]
Daniel Kifer, Adam Smith, and Abhradeep Thakurta.
Private convex empirical risk minimization and high-dimensional
In Conference on Learning Theory, pages 25–1, 2012.

[24]
GM Korpelevich.
The extragradient method for finding saddle points and other
Matecon, 12:747–756, 1976.

[25]
Georgios Kotsalis, Guanghui Lan, and Tianjiao Li.
Simple and optimal methods for stochastic variational inequalities,
i: operator extrapolation, 2020.

[26]
Yunwen Lei, Zhenhuan Yang, Tianbao Yang, and Yiming Ying.
Stability and generalization of stochastic gradient methods for
minimax problems.
In International Conference on Machine Learning, pages
6175–6186. PMLR, 2021.

[27]
Tian Li, Maziar Sanjabi, Ahmad Beirami, and Virginia Smith.
Fair resource allocation in federated learning.
In 8th International Conference on Learning Representations,
ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020. OpenReview.net, 2020.

[28]
A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro.
Robust stochastic approximation approach to stochastic programming.
SIAM Journal on Optimization, 19(4):1574–1609, 2009.

[29]
Arkadi Nemirovski.
Prox-method with rate of convergence o(1/t) for variational
inequalities with lipschitz continuous monotone operators and smooth
convex-concave saddle point problems.
SIAM Journal on Optimization, 15(1):229–251, 2004.

[30]
A.S. Nemirovsky and D.B. Yudin.
Problem Complexity and Method Efficiency in Optimization.
J. Wiley @ Sons, New York, 1983.

[31]
Yurii Nesterov.
Dual extrapolation and its applications to solving variational
inequalities and related problems.
Mathematical Programming, 109(2):319–344, 2007.

[32]
J von Neumann.
Zur theorie der gesellschaftsspiele.
Mathematische Annalen, 100:295–320, 1928.

[33]
R. Tyrrell Rockafellar.
Monotone operators and the proximal point algorithm.
SIAM Journal on Control and Optimization, 14(5):877–898, 1976.

[34]
Shai Shalev-Shwartz, Ohad Shamir, Nathan Srebro, and Karthik Sridharan.
Learnability, stability and uniform convergence.
The Journal of Machine Learning Research, 11:2635–2670, 2010.

[35]
Maurice Sion.
On general minimax theorems.
Pacific Journal of Mathematics, 8(1):171–176, 1958.

[36]
Thomas Steinke and Jonathan R. Ullman.
Between pure and approximate differential privacy.
J. Priv. Confidentiality, 7(2), 2016.

[37]
Jonathan Ullman.
Private multiplicative weights beyond linear queries.
In Proceedings of the 34th ACM SIGMOD-SIGACT-SIGAI Symposium on
Principles of Database Systems, pages 303–312. ACM, 2015.

[38]
Roman Vershynin.
High-Dimensional Probability: An Introduction with Applications
in Data Science.
Cambridge Series in Statistical and Probabilistic Mathematics.
Cambridge University Press, 2018.

[39]
Robert Williamson and Aditya Menon.
Fairness risk measures.
In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning,
volume 97 of Proceedings of Machine Learning Research, pages
6786–6797. PMLR, 09–15 Jun 2019.

[40]
Zhenhuan Yang, Shu Hu, Yunwen Lei, Kush R. Varshney, Siwei Lyu, and Yiming
Differentially private SGDA for minimax problems.
CoRR, abs/2201.09046, 2022.

[41]
Junyu Zhang, Mingyi Hong, Mengdi Wang, and Shuzhong Zhang.
Generalization bounds for stochastic saddle point problems.
In Proceedings of The 24th International Conference on
Artificial Intelligence and Statistics, volume 130 of Proceedings of
Machine Learning Research, pages 568–576. PMLR, 13–15 Apr 2021.

§ PROOF OF PROPOSITION <REF>

Let $\bS' = (\bbeta'_1, \dots, \bbeta'_n)$ be independent of $\bS$. For $i \in [n]$, we denote $\bS^{i} = (\bbeta_1, \dots, \bbeta_{i-1}, \bbeta'_{i},  \bbeta_{i+1}, \dots, \bbeta_n)$. Then, for any $w \in \W$, we have

\begin{align}
\EE_{\bS}\inprod{F(w)}{\A(\bS) - w} &= \EE_{\bS, \bS'}\frac{1}{n}\sum_{i = 1}^n \inprod{F_{\bbeta'_{i}}(w)}{\A(\bS)-w} \nonumber\\
&=\EE_{\bS, \bS'}\frac{1}{n}\sum_{i = 1}^n \inprod{F_{\bbeta_{i}}(w)}{\A(\bS^{i})-w} \nonumber\\
&=\EE_{\bS, \bS'}\frac{1}{n}\sum_{i = 1}^n \inprod{F_{\bbeta_{i}}(w)}{\A(\bS)-w} + \inprod{F_{\bbeta_{i}}(w)}{\A(\bS^{i}) - \A(\bS)} \nonumber\\
&\le \EE_{\bS, \bS'}\frac{1}{n}\sum_{i = 1}^n \inprod{F_{\bbeta_{i}}(w)}{\A(\bS)-w} + \gnorm{F_{\bbeta_{i}}(w)}{}{}
\gnorm{\A(\bS^{i}) - \A(\bS)}{}{} \nonumber\\
&\le \EE_{\bS, \bS'}\frac{1}{n}\sum_{i = 1}^n \inprod{F_{\bbeta_{i}}(w)}{\A(\bS)-w} + M\gnorm{\A(\bS^{i}) - \A(\bS)}{}{} \nonumber \\%\quad\text{a.s.}
&=\EE_{\bS }\inprod{F_\bS(w)}{\A(\bS) - w} + \EE_{\bS, \bS'} \frac{1}{n}\sum_{i = 1}^n M\gnorm{\A(\bS^{i}) - \A(\bS)}{}{}
\label{eq:a.s_weak_bound}
\end{align}

Now, taking supremum over $w \in \W$ and taking expectation over $\A$ which is $\delta$-UAS, we have,

\begin{align*}
\EE_{{\cal A}}[\wVIgap{\A(\bS)}{F}] &\le \EE_{\A}[\sup_{w \in \W}\EE_{\bS}\inprod{F_{\bS}(w) }{\A(\bS) -w}] + \EE_{\bS, \bS'}\frac{1}{n}\sum_{i = 1}^n M \EE_{\A}\gnorm{\A(\bS^i) - \A(\bS)}{}{}\\
&\le\EE_{\A,\bS}[\sup_{w \in \W}\inprod{F_{\bS}(w) }{\A(\bS) -w}] + \EE_{\bS, \bS'}\frac{1}{n}\sum_{i = 1}^n M \EE_{\A}\gnorm{\A(\bS^i) - \A(\bS)}{}{}\\
&\le \EE_{\bS}[\EmpVIgap{\A}{F_{\bS}}] + M\delta.
\end{align*}

§ OPERATOR EXTRAPOLATION METHOD [25]

Suppose we want to solve VI problem associated with operator $F_k(\cdot) = F(\cdot)+ \lambda_k(\cdot-w_k)$ whose (unique) solution be $w_{k+1}^*$. It is clear that $F_k$ is an $(L+\lambda_k)$-Lipschitz continuous operator which is $\lambda_k$-strongly monotone as well. Denote $\kappa := \frac{\lambda_k}{L+\lambda_k} + 1$ and consider the following algorithm for solving this problem:
Operator Extrapolation (OE) method

[1]
Let $z_0 = z_1 = w_k$ be given.
$t = 1, \dots, T$
$z_{t+1} = \argmin_{w \in \W} \frac{1}{2(L+\lambda_k)} \inprod{F_k(z_t) + \frac{1}{\kappa}\big[F_k(z_t)-F_k(z_{t-1}) \big] }{w} + \frac{1}{2}\gnorm{w-z_t}{}{2}$

We have the following convergence guarantee for this algorithm:
$$\gnorm{z_{T}-w^*_{k+1}}{}{2} \le \kappa^{-T}\gnorm{z_1 - w^*_{k+1}}{}{2}$$

In particular, in order to ensue that $\gnorm{z_T-w_{k+1}^*}{}{} \le \frac{\nu}{LD^2+MD+ 2\lambda_kD^2}$, we require
$$T = 2\kappa \ln\Big(\frac{LD^2+MD+2\lambda_kD^2}{\nu}\Big) = \frac{2(L+2\lambda_k)}{L+\lambda_k}\ln\Big(\frac{LD^2+MD+2\lambda_kD^2}{\nu}\Big)$$


# Updated analysis of near-threshold heavy quarkonium production
for probe of proton’s gluonic gravitational form factors

Yuxun Guo<EMAIL_ADDRESS>Department of Physics, University of Maryland,
College Park, MD 20742, USA  Nuclear Science Division, Lawrence Berkeley
National Laboratory, Berkeley, CA 94720, USA  Xiangdong Ji<EMAIL_ADDRESS>Department of Physics, University of Maryland, College Park, MD 20742, USA
Yizhuang Liu<EMAIL_ADDRESS>Institute of Theoretical Physics,
Jagiellonian University, 30-348 Kraków, Poland  Jinghong Yang<EMAIL_ADDRESS>Department of Physics, University of Maryland, College Park, MD 20742, USA

###### Abstract

There has been growing interest in the near-threshold production of heavy
quarkonium which can access the gluonic structure in the nucleon. Previously
[1] we studied this process with quantum chromodynamics (QCD) and showed that
it can be factorized with the gluon generalized parton distributions (GPDs) in
the heavy quark limit. We further argued that the hadronic matrix element is
dominated by its leading moments corresponding to the gluonic gravitational
form factors (GFFs) in this limit. Since then, there have been many new
developments on this subject. More experimental measurements have been made
and published, and the lattice simulation of gluonic GFFs has been improved as
well. In this work, we make an important revision to a previous result and
perform an updated analysis with the new inputs. We also study the importance
of the large momentum transfer to extract these gluonic structures reliably in
this framework.

##  I Introduction

The gluonic structures have been an important topic in quantum chromodynamics
(QCD) and hadronic physics. On the one hand, gluons that mediate the strong
interaction play a prominent role at the non-perturbative scale in the hadron.
On the other hand, they are free of electroweak interaction and much harder to
probe than quarks. Consequently, there has been rising interest in the
exclusive electro/photo-productions of heavy quarkonium. Assuming suppressed
contributions from the intrinsic heavy quarks, these processes are dominated
by the exchange of gluons and thus can be used to access the gluonic
structures. Experiments with such purposes are planned in the future Electron-
Ion Collider (EIC) [2], whereas at Jefferson Laboratory (JLab) the near-
threshold production of $J/\psi$ that requires less energy has been proposed
and measured which revealed exciting results [3, 4, 5, 6].

There have been many theoretical developments in the literature to analyze
this process [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. In
the previous work of three of the authors [1], we showed that the near-
threshold photoproduction of heavy quarkonium can be factorized with gluon
generalized parton distributions (GPDs), extending the factorization proved
for the diffractive production in the collinear limit [23, 24]. Utilizing the
heavy quark limit, we also argued that the hadronic matrix element will be
dominated by the leading moments that correspond to the gluonic gravitational
form factors (GFFs) near the threshold. Although a complete determination of
the gluonic GFFs from these measurements alone is still model-dependent [5],
it provides us with an effective tool to handle such problems.

In this work, we perform an updated analysis of the near-threshold $J/\psi$
production following the previous work for two main reasons. First, we note
that there was a missing factor of 2 in the hadronic matrix element $G(\xi,t)$
in ref. [1] due to a mix of conventions. Besides, there have been many new
developments on this subject ever since our previous work. More data have been
measured and published recently by the $J/\psi$ 007 experiment at JLab Hall C
[5] and the GlueX collaboration at Hall D [6]. Moreover, the lattice QCD
simulation of gluonic GFFs has been improved [25]. The large $\xi$ expansion
essential for relating this process to the gluonic GFFs requires careful
treatments. Consequently, we will present more detailed analyses of these new
results in this work.

##  II Revised formula

We start by noting that the definition of gluon GPDs in eq. (16) of the
previous work [1] includes an extra factor of $\frac{1}{2}$ from the trace
operator Tr, which should be removed. Correspondingly, the extra factor of
$\frac{1}{2}$ in eq. (24) in ref. [1] shall be removed as well. The revision
will leave the main cross-section formula unchanged, i.e., eqs. (13) and (17)
in ref. [1]. However, one needs to substitute the gluon GPDs with the
corrected definition that reads,

$\displaystyle F_{g}(x,\xi,t)\equiv$
$\displaystyle\frac{1}{(\bar{P}^{+})^{2}}\int\frac{\text{d}\lambda}{2\pi}e^{i\lambda
x}\left<P^{\prime}\left|F^{a+i}_{\;\;\;\;\;}\left(-\frac{\lambda
n}{2}\right)F^{a+}_{\;\;i}\left(\frac{\lambda n}{2}\right)\right|P\right>\ ,$
(1)

where the index $a$ sums over all colors. The $F_{g}(x,\xi,t)$ can be
parameterized as [26, 27]

$\displaystyle
F_{g}(x,\xi,t)=\frac{1}{2\bar{P}^{+}}\bar{u}(P^{\prime})\left[H_{g}\gamma^{+}+E_{g}\frac{i\sigma^{+\alpha}\Delta_{\alpha}}{2M_{N}}\right]u(P)\
,$ (2)

where $H_{g}$ and $E_{g}$ are the well-known $H_{g}(x,\xi,t)$ and
$E_{g}(x,\xi,t)$ GPDs. The hadronic matrix element $G(\xi,t)$ remains to be
[1]

$\displaystyle G(t,\xi)=\frac{1}{2\xi}\int_{-1}^{1}\text{d}x{\cal
A}(x,\xi)F_{g}(x,\xi,t)\ ,$ (3)

where the Wilson coefficient ${\cal A}(x,\xi)$ reads

${\cal A}(x,\xi)\equiv\frac{1}{x+\xi-i0}-\frac{1}{x-\xi+i0}\ .$ (4)

On the other hand, the extra factor of $\frac{1}{2}$ does affect the relations
between the gluon GPDs and GFFs. Removing the extra factor of $\frac{1}{2}$ in
the gluon GPDs causes an extra factor of 2 to be multiplied to each gluonic
GFF when expanding the GPDs in terms of their moments. Consequently, in the
unpolarized case we have

$\displaystyle|G(t,\xi)|^{2}=\frac{4}{\xi^{4}}\Bigg{\\{}\left(1-\frac{t}{4M_{N}^{2}}\right)E_{2}^{2}$
$\displaystyle-2E_{2}(H_{2}+E_{2})+\left(1-\xi^{2}\right)(H_{2}+E_{2})^{2}\Bigg{\\}}\
,$ (5)

where an extra factor of $4$ appears accordingly, compared to the eq. (26) in
ref. [1]. The same factor of $4$ should be multiplied to the polarized
$|G(t,\xi)|^{2}$ in eqs. (58) and (59) in ref. [1] as well. Recall that the
GFFs $H_{2}\equiv H_{2}(t,\xi)$ and $E_{2}\equiv E_{2}(t,\xi)$ follow the same
definition as ref. [1] that reads

$\displaystyle\int_{0}^{1}dxH_{g}(x,\xi,t)=A_{2,0}^{g}(t)+(2\xi)^{2}C^{g}_{2}\equiv
H_{2}(t,\xi)\ ,$
$\displaystyle\int_{0}^{1}dxE_{g}(x,\xi,t)=B^{g}_{2,0}(t)-(2\xi)^{2}C^{g}_{2}\equiv
E_{2}(t,\xi)\ .$ (6)

Figure 1: A comparison of the total cross-sections predicted with gluonic GFFs
from the improved lattice simulation [25] to the two measurements by the GlueX
collaboration in 2019 [4] and 2023 [6]. The $J/\psi$ 007 experiment [5] did
not measure the total cross-section directly and will be discussed in the next
section for the differential cross-section analysis.

In FIG. 1, we compare the total cross-sections predicted with the improved
lattice simulation of the gluonic GFFs [25] to the two measurements by the
GlueX collaboration in 2019 [4] and 2023 [6] respectively without tuning any
parameters. The consistency seems better compared to the one in previous work
before the revision [1]. More detailed analyses regarding the extraction of
GFFs with differential cross-section measurements by both the GlueX
collaboration and the $J/\psi$ 007 experiment [5] will be presented in the
next section.

##  III  Analyses of the near-threshold $J/\psi$ production data

To extract the gluonic GFFs from the near-threshold $J/\psi$ photoproduction
measurements, recall that the cross-section formula in the previous work reads
[1]:

$\begin{split}\frac{d\sigma}{dt}=\frac{\alpha_{\rm
EM}e_{Q}^{2}}{4\left(W^{2}-M_{N}^{2}\right)^{2}}\frac{(16\pi\alpha_{S})^{2}}{3M_{V}^{3}}|\psi_{\rm
NR}(0)|^{2}|G(t,\xi)|^{2}\ .\end{split}$ (7)

The same formula applies after the revision except that the revised hadronic
matrix element $G(t,\xi)$ in eq. (II) has an extra factor of $4$.

We parameterize the two GFFs $A_{g}(t)$ and $C_{g}(t)$ in tripole forms:

$\displaystyle A_{g}(t)$
$\displaystyle=\frac{A_{g}(0)}{\left(1-\frac{t}{m_{A}^{2}}\right)^{3}}\ ,$ (8)
$\displaystyle C_{g}(t)$
$\displaystyle=\frac{C_{g}(0)}{\left(1-\frac{t}{m_{C}^{2}}\right)^{3}}\ ,$ (9)

ignoring the $B_{g}(t)$. The forward $A_{g}(0)$ is fixed according to the
gluon PDF from global analysis to be 0.414 [28]. Then we are left with three
parameters: $m_{A}$, $C_{g}(0)$, and $m_{C}$ to be determined from the near-
threshold production measurements, for which we consider the combination of
the recently published data from the $J/\psi$ 007 experiment [5] and the GlueX
collaboration [6].

Figure 2: Differential cross-section data points together with contours of
equal $\xi$ on the $(E_{b},-t)$ plane in the kinematically allowed region with
$M_{J/\psi}=3.097$ GeV. Each dot represents a data point from the $J/\psi$ 007
experiment at JLab Hall C (circle) [5] or GlueX collaboration at Hall D
(square) [6].

Before moving on to the data analysis, we must first comment that the relation
between the near-threshold $J/\psi$ production cross-sections and the gluonic
GFFs is justified in the heavy quark limit only, where the momentum transfer
squared $|t|$ approaches infinity and the skewness parameter $\xi$ approaches
1. Thus, the extraction of the gluonic GFFs relies on the so-called large
$\xi$ expansion that is most applicable in the $\xi\to 1$ limit. However, we
can only have measurements with finite momentum transfer squared $|t|$ and
skewness $\xi<1$ in reality — the $J/\psi$ 007 measurements only cover the
region with $\xi<0.6$ whereas the GlueX measurements extend to larger-$\xi$
region but with limited data points as shown in FIG. 2. Out of the 124
differential cross-section data points combining the $J/\psi$ 007 and GlueX
measurements, 85 have $\xi>0.4$ while only 33 have $\xi>0.5$. Also associated
with the increasing $\xi$ is the decreasing quality in the data due to the
lack of events at large $|t|$, putting additional challenge to the reliable
extraction of the gluonic GFFs.

Accordingly, we can either focus on the large-$|t|$ data where the statistical
uncertainties will be huge, or consider the medium or even lower-$|t|$ data as
well where the systematical uncertainties could be the dominant effect. These
possibilities will be studied in more details in the following subsections.

###  III.1 Analysis with large-$|t|$ data

We start with a rather extreme option where we consider the data with large
momentum transfer squared $|t|$ only, i.e., we select the data with
$\xi>0.5$111The large $|t|$ and large $\xi$ conditions can be used
interchangeably near the threshold as shown in FIG. 2 — it requires a minimum
$|t|$ to reach a certain $\xi$ and larger $\xi$ requires larger $|t|$.. As
mentioned above, there are only 33 data points in this region, which are too
few to determine all three parameters $m_{A}$, $C_{g}(0)$, and $m_{C}$.
Indeed, by fitting all three of them to the data with $\xi>0.5$ by a standard
$\chi^{2}$ analysis with the iminuit Python interface of Minuit2 package [29,
30], we obtain $m^{\xi>0.5}_{A}=2.38\pm 0.08\text{ GeV}$,
$C^{\xi>0.5}_{g}(0)=-16\pm 34$ and $m^{\xi>0.5}_{C}=0.60\pm 0.26\text{ GeV}$.
The unreasonably large best-fit value and statistical uncertainty of
$C_{g}(0)$ indicate that the $C_{g}(t)$ form factor cannot be effectively
constrained with the large-$|t|$ data. The $A_{g}(t)$ form factor, on the
other hand, is better constrained benefiting from the forward constraint from
gluon PDFs. Consequently, we have to utilize other information about the
gluonic GFFs, e.g., the ones from lattice simulations [31, 25], as the
reference values to fix the undetermined parameters and avoid potential
overfitting.

Set |  $A_{g}(0)$ |  $m_{A}(\text{GeV})$ |  $C_{g}(0)$ |  $m_{C}(\text{GeV})$ |  red. $\chi^{2}$
---|---|---|---|---|---
ref. | 0.414 | 1.64 | -0.48 | 1.07 | -
1 | fixed | 2.38(08) | -16(34) | 0.60(26) | 1.34
2 | fixed | 2.53(12) | fixed | 1.41(09) | 1.39
3 | fixed | 2.46(08) | -1.20(20) | fixed | 1.33
4 | fixed | fixed | 0.014(03) | 13(14) | 1.25
5 | fixed | 2.14(03) | fixed | fixed | 1.70
6 | fixed | fixed | 0.29(06) | fixed | 2.63
Table 1: A summary of the best-fit parameters of the $\xi>0.5$ data with
reduced $\chi^{2}$ fitting to the combined differential cross-section data
from the $J/\psi$ 007 experiment [5] and the GlueX collaboration [6]. The
reference values are from the gluon PDF [28] for the $A_{g}(0)$ and lattice
simulations [25] for the other three. Parameters listed as “fixed” are fixed
to be the reference values.

In TABLE 1, we summarize the results from various fits that fix the parameters
differently, where in the first row we list the reference values from the
global-fitted gluon PDF for the $A_{g}(0)$ [28] and lattice simulations of
GFFs [25] for the other three parameters. Among all the fits, set 2 and 3 seem
more realistic that exhibit no signs of overfitting and describe the data fine
with reduced $\chi^{2}$s around 1.3. On the other hand, both set 1 and 4 have
undetermined parameters with unreasonably large values and uncertainties
though their reduced $\chi^{2}$s look fine. Such observation indicates the
potential overfitting in these sets. As for set 5 and 6, their large reduced
$\chi^{2}$s imply that a one-parameter fit may not be able to describe the
large-$|t|$ data well.

We also comment on the size of reduced $\chi^{2}$s in the fits which seems
rather large for fits with only 33 data points. The main reason is the
anomalously rising behavior in the $t$-dependence of the measured differential
cross-sections by the GlueX collaboration, which has been discussed with more
details in ref. [6] (see for instance the FIG. 13 therein). Such behaviors
were not observed in the $J/\psi$ 007 experiment, partially due to its limit
kinematical coverage at large $|t|$ as shown in FIG. 2. Other than that, the
GlueX measurements are in good agreement with the measurements by the $J/\psi$
007 experiment. Also note that the GlueX data have about $20\%$ normalization
uncertainties (about $4\%$ for the $J/\psi$ 007 experiment) that we did not
include since the two data sets seem consistent except for this rising
$t$-dependence observed by the GlueX. This would lead to lower reduced
$\chi^{2}$ if included.

Figure 3: A comparison of the $A_{g}(t)$ (up) and $C_{g}(t)$ (down) form
factor extracted with the differential cross-section data from the $J/\psi$
007 experiment only [5] selecting $\xi>0.5$ and the lattice simulation results
[25]. Bands correspond to $1\sigma$ statistical uncertainties.

For the analysis here, when we select large-$|t|$ data, such behaviors
observed in the large-$|t|$ region will be more pronounced, which cannot be
described well by the tripole form used here. This caused the overfitting
behavior observed for the $C_{g}(0)$ here, and it will continue to make the
reduced $\chi^{2}$s larger even when we extend to include the medium or
lower-$|t|$ data as we will show in the following subsections. To check this
argument, we perform the same fits but to the $J/\psi$ 007 data only, and
obtain much lower reduced$\chi^{2}$s for the fits here and in the following
subsections. For instance, we obtain $m^{\xi>0.5}_{A}=2.10\pm 0.29\text{
GeV}$, $C^{\xi>0.5}_{g}(0)=-0.5\pm 2.6$ and $m^{\xi>0.5}_{C}=1\pm 1.6\text{
GeV}$ with a reduced $\chi^{2}$ of 0.46 when fitting to only the $J/\psi$ 007
data with $\xi>0.5$. Although the relative uncertainties are larger due to the
reduced sample size, the overfitting behavior seems to disappear. In FIG. 3,
we also present the corresponding extracted GFFs and compare them with the
lattice results. We note that the rising $t$-dependence of the measured
differential cross-section by the GlueX collaboration might be due to other
sources beyond the scope of this work.

With these results, we remark that with the limited data at large $|t|$, not
much information can be extracted about the gluonic GFFs, especially the
$C_{g}(t)$ form factor. Excluding the fits with undetermined parameters or
large reduced $\chi^{2}$s, the remaining ones cannot constrain the $C_{g}(t)$
form factor well unless when one of the parameters is fixed. While the
$A_{g}(t)$ form factor is better constrained owing to the forward constraint
from the gluon PDF, the extracted tripole mass $m_{A}$ still depends on the
values of the fixed parameters, the $C_{g}(0)$ particularly. Since the lattice
simulation of the $C_{g}(t)$ form factor has the largest uncertainties at
$t=0$ [25], the reliability of the extracted $m_{A}$ would be affected
accordingly. Therefore, it is crucial to obtain more data with higher quality
in the analysis to better constrain the gluonic GFFs.

###  III.2 Analysis including medium-$|t|$ data

Given the limited constraining power with only the large-$|t|$ data due to the
insufficient data, one compromising choice is to include more data with medium
$|t|$. As we decrease the cut in $|t|$ or the skewness $\xi$, we effectively
reduce the statistical uncertainties while enhancing the systematical
uncertainties. Then, the general philosophy is to find the cut that balances
the two uncertainties to maximally utilize the data. Although the proper
estimation of the systematical uncertainties from the higher order effects and
large $\xi$ expansion will be extremely involving and beyond the scope of this
work, we could still include the medium-$|t|$ data in the fit and study their
effects for discussion.

Set |  $A_{g}(0)$ |  $m_{A}(\text{GeV})$ |  $C_{g}(0)$ |  $m_{C}(\text{GeV})$ |  red. $\chi^{2}$
---|---|---|---|---|---
ref. | 0.414 | 1.64 | -0.48 | 1.07 | -
1 | fixed | 2.07(05) | -1.21(37) | 0.91(10) | 1.42
2 | fixed | 2.25(05) | fixed | 1.38(04) | 1.64
3 | fixed | 2.14(03) | -0.81(03) | fixed | 1.43
4 | fixed | 1.88(01) | fixed | fixed | 2.38
5 | fixed | fixed | -0.29(02) | fixed | 3.74
Table 2: A summary of the best-fit parameters of the $\xi>0.4$ data with
reduced $\chi^{2}$ fitting to the combined differential cross-section data
from the $J/\psi$ 007 experiment [5] and the GlueX collaboration [6]. The
reference values are from the gluon PDF [28] for the $A_{g}(0)$ and lattice
simulations [25] for the other three. Parameters listed as “fixed” are fixed
to be the reference values.

With that in mind, we repeat the above analysis with the 85 data when
selecting $\xi>0.4$. The extended set of data does constrain the three
parameters better. Consequently, we obtain $m^{\xi>0.4}_{A}=2.07\pm 0.05\text{
GeV}$, $C^{\xi>0.4}_{g}(0)=-1.2\pm 0.4$, and $m^{\xi>0.4}_{C}=0.91\pm
0.10\text{ GeV}$ from a three-parameter fit. This, together with the other
fits that fix the parameters differently, is summarized in TABLE 2 as the set
1. Unlike the previous case where many of the parameters cannot be determined,
no signs of overfitting are observed here by virtue of the extra inputs. Among
all the fits, set 1 has the lowest reduced $\chi^{2}$ as expected. Meanwhile,
since the $m_{C}$ obtained in set 1 is close to the reference value, the two-
parameter fit in set 3 that fixes the $m_{C}$ looks similar. The other two-
parameter fit in set 2 shows slightly worse results, whereas the one-parameter
fits in set 4 and 5 do not seem to work. Thus, we consider the three-parameter
fit in set 1 for further analysis.

Figure 4: A comparison of the $A_{g}(t)$ (up) and $C_{g}(t)$ (down) form
factor extracted with the combined differential cross-section data from the
$J/\psi$ 007 experiment [5] and the GlueX collaboration [6] selecting
$\xi>0.4$ and the lattice simulation results [25]. Bands correspond to
$1\sigma$ statistical uncertainties.

In FIG. 4, we compare the GFFs extracted in set 1 with the ones from lattice
simulation [25]. Still, a negative $C_{g}(0)$ is favored when including the
medium-$|t|$ data which is consistent with the lattice simulation. The
uncertainty associated with the extracted $C_{g}(0)$ has been improved though
still sizable, whereas the $A_{g}(t)$ form factor remains to be better
constrained due to the forward constraint from the gluon PDF. Thus, we note
that when including the medium-$|t|$ data the observation is consistent with
the large-$|t|$ extraction. However, the gluonic GFFs are better constrained
with the extra data, allowing us to roughly extract the gluon GFFs. Keeping in
mind that there are still model-dependence and systematical uncertainties to
be clarified in the extraction, we shall note that, as mentioned above, more
data with higher quality at large/medium $|t|$ are crucial for better
determination of the gluonic GFFs in this framework.

Figure 5: A comparison of the total cross-sections calculated based on the fit
to the combined differential cross-section data from the $J/\psi$ 007
experiment [5] and the GlueX collaboration [6] selecting $\xi>0.4$ and the two
GlueX measurements [4, 6]. The band corresponds to $1\sigma$ statistical
uncertainties.

In FIG. 5, we also compare the total cross-section calculated based on the fit
in set 1 to the two measurements by the GlueX collaboration in 2019 [4] and
2023 [6]. Although these total cross-section data are not considered in the
fit, they still show good agreement.

###  III.3 Analysis with full data

In the last subsection, we also present the fit to all the 124 differential
cross-section data. Since about 1/3 of these data have $\xi<0.4$ and they will
be weighted more than the others in the fit resulting from their lower
relative uncertainties, the systematical uncertainties in this case can get
out of control. Therefore, we consider it as an exercise to illustrate the
effect of the lower-$|t|$ data in the analysis. It has to be kept in mind that
these results should NOT be taken too seriously due to the potentially large
systematical uncertainties associated.

Figure 6: A comparison of the $A_{g}(t)$ (up) and $C_{g}(t)$ (down) form
factor extracted with all the combined differential cross-section data from
the $J/\psi$ 007 experiment [5] and the GlueX collaboration [6] and the
lattice simulation results [25]. Bands correspond to $1\sigma$ statistical
uncertainties.

Being cautious of that, we fit the three parameters $m_{A}$, $C_{g}(0)$, and
$m_{C}$ to all the data, and obtain $m_{A}=1.88\pm 0.03\text{ GeV}$,
$C_{g}(0)=-1.49\pm 0.27$, and $m_{C}=0.8\pm 0.05\text{ GeV}$ with the reduced
$\chi^{2}$ as 1.89. An explicit comparison of the extracted GFFs with the
lattice simulation [25] is shown in FIG. 6. Quite surprisingly, the extracted
GFFs agree well with the lattice simulation, which is in accord with the
previous finding in FIG. 1 that the revised cross-section formula produces
nice consistency between the lattice simulations and the near-threshold
$J/\psi$ production cross-sections. However, since the GFFs from lattice
simulation are mostly in the small-$|t|$ region where the systematical
uncertainty could be large, the reliability of these statements might be
challenged by that accordingly.

##  IV Extracted gluonic mass and scalar radii

Figure 7: Comparisons of the extracted gluonic mass (up) and scalar (down)
radii. The four values from left to right are from lattice simulation [25],
holographic QCD extraction [20] with the $J/\psi$ 007 measurements [5], the
extraction here with $\xi>0.4$ data combining the $J/\psi$ 007 experiment and
the GlueX measurements [6], and the extraction with all data combining the
$J/\psi$ 007 experiment and the GlueX measurements. The gray bands indicate
potentially large systematical uncertainties besides the statistical ones.

To make the comparison more intuitive, in this section we also consider the
extracted gluonic mass and scalar radii of the proton. They are defined as
[32],

$\displaystyle\left\langle r^{2}_{m}\right\rangle_{g}$
$\displaystyle=\frac{1}{A_{g}(0)}\left[6\left.\frac{dA\left(t\right)}{dt}\right|_{t=0}-6\leavevmode\nobreak\
\frac{C(0)}{M_{N}^{2}}\right]\ ,$ (10) $\displaystyle\left\langle
r^{2}_{s}\right\rangle_{g}$
$\displaystyle=\frac{1}{A_{g}(0)}\left[6\left.\frac{dA\left(t\right)}{dt}\right|_{t=0}-18\leavevmode\nobreak\
\frac{C(0)}{M_{N}^{2}}\right]\ ,$

and renormalized with the gluonic momentum fraction $A_{g}(0)$. In FIG. 7, we
compare the gluonic mass and scalar radii extracted here with the lattice
simulation results [25] as well as the holographic QCD extraction [20] with
the $J/\psi$ 007 measurements [5]. Consistent with the arguments before, as we
increase the cut in $\xi$, there will be fewer measurements with larger
uncertainties, leading to more uncertain mass and scalar radii
correspondingly. Actually, it is apparent with eq. (10) that the uncertainties
in the extracted radii are mainly from the uncertainties in $C_{g}(0)$, as the
$A_{g}(0)$ is fixed by the gluon PDF. On the other hand, the fits with lower
cut in $\xi$ will be associated with larger systematical uncertainties, e.g.,
in the fit with all data.

Therefore, we consider the extraction with $\xi>0.4$ data as the reference
values for this work, which gives $\left<r_{m}\right>_{g}=0.97\pm
0.12\rm{\leavevmode\nobreak\ fm}$ and $\left<r_{s}\right>_{g}=1.58\pm
0.22\rm{\leavevmode\nobreak\ fm}$. Note that the sizable statistical
uncertainties here may still be accompanied by comparable systematical
uncertainties. Accordingly, more measurements with higher quality, especially
at large $|t|$ are of vital importance to further improve the extraction of
the gluonic GFFs and the mass/scalar radii from the near-threshold $J/\psi$
production measurements.

We also notice that the scalar and mass radii extracted with the holographic
QCD approach appear to agree well with the lattice simulation with a smaller
statistical uncertainty compared to ours. There are two reasons that account
for that. First, the holographic approach works better in the Regge limit
$t\to 0$ where more data with higher quality exist. Besides, there is an
undetermined normalization constant in the holographic QCD approach that is
manually fixed [20], with which the statistical/systematical uncertainties
associated are not accounted. Therefore, one could consider the difference
between the gluonic GFFs and the scalar/mass radii extracted with the
holographic QCD approach and the ones here as an estimation of the
systematical uncertainties or model dependence.

##  V Summary and outlook

To summarize, in this work we revise the previous cross-section formula [1]
with the factor of 2 mismatch, and perform an updated analysis with the latest
lattice simulation of gluonic GFFs [25] and the recently published data from
$J/\psi$ 007 experiment [5] and GlueX collaboration [6]. We show that with the
revised formula the agreement between the gluonic GFFs from lattice simulation
and extraction with GPD factorization gets improved. On the other hand, we
also argue that this framework requires large momentum transfer squared $|t|$
and skewness $\xi$, and thus we perform a series of analyses to properly
address their effects. We show that the gluonic GFFs can be roughly
constrained with a cut of $\xi>0.4$. However, it is crucial to have more high-
quality data at large $|t|$ to improve the extraction.

The critical future developments include studying the higher order corrections
such as the next-to-leading order effects in the strong coupling $\alpha_{S}$
as well as the finite quarkonium mass correction in $M_{p}/M_{V}$. In
addition, the systematical uncertainties from the large $\xi$ expansion also
require proper treatment.

Acknowledgment: We thank F. Yuan for useful discussions and correspondences.
We particularly thank Z.-E. Meziani for sharing the data of the recent
$J/\psi$-007 experiment. This research is partly supported by the U.S.
Department of Energy, Office of Nuclear Physics, under contract number DE-
SC0020682.

## References

  * Guo  _et al._ [2021] Y. Guo, X. Ji, and Y. Liu, Phys. Rev. D 103, 096010 (2021), arXiv:2103.11506 [hep-ph] .
  * Accardi  _et al._ [2016] A. Accardi _et al._ , Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Joosten and Meziani [2018] S. Joosten and Z. E. Meziani, PoS QCDEV2017, 017 (2018), arXiv:1802.02616 [hep-ex] .
  * Ali  _et al._ [2019] A. Ali _et al._ (GlueX), Phys. Rev. Lett. 123, 072001 (2019), arXiv:1905.10811 [nucl-ex] .
  * Duran  _et al._ [2023] B. Duran _et al._ , Nature 615, 813 (2023), arXiv:2207.05212 [nucl-ex] .
  * Adhikari  _et al._ [2023] S. Adhikari _et al._ , (2023), arXiv:2304.03845 [nucl-ex] .
  * Voloshin [1979] M. B. Voloshin, Nucl. Phys. B 154, 365 (1979).
  * Gottfried [1978] K. Gottfried, Phys. Rev. Lett. 40, 598 (1978).
  * Appelquist and Fischler [1978] T. Appelquist and W. Fischler, Phys. Lett. B 77, 405 (1978).
  * Bhanot and Peskin [1979] G. Bhanot and M. E. Peskin, Nucl. Phys. B 156, 391 (1979).
  * Kharzeev [1996] D. Kharzeev, Proc. Int. Sch. Phys. Fermi 130, 105 (1996), arXiv:nucl-th/9601029 .
  * Kharzeev  _et al._ [1999] D. Kharzeev, H. Satz, A. Syamtomov, and G. Zinovjev, Eur. Phys. J. C 9, 459 (1999), arXiv:hep-ph/9901375 .
  * Kharzeev [2021] D. E. Kharzeev, Phys. Rev. D 104, 054015 (2021), arXiv:2102.00110 [hep-ph] .
  * Gryniuk and Vanderhaeghen [2016] O. Gryniuk and M. Vanderhaeghen, Phys. Rev. D 94, 074001 (2016), arXiv:1608.08205 [hep-ph] .
  * Gryniuk  _et al._ [2020] O. Gryniuk, S. Joosten, Z.-E. Meziani, and M. Vanderhaeghen, Phys. Rev. D 102, 014016 (2020), arXiv:2005.09293 [hep-ph] .
  * Hatta and Yang [2018] Y. Hatta and D.-L. Yang, Phys. Rev. D 98, 074003 (2018), arXiv:1808.02163 [hep-ph] .
  * Hatta  _et al._ [2019] Y. Hatta, A. Rajan, and D.-L. Yang, Phys. Rev. D 100, 014032 (2019), arXiv:1906.00894 [hep-ph] .
  * Mamo and Zahed [2020] K. A. Mamo and I. Zahed, Phys. Rev. D 101, 086003 (2020), arXiv:1910.04707 [hep-ph] .
  * Mamo and Zahed [2021] K. A. Mamo and I. Zahed, Phys. Rev. D 103, 094010 (2021), arXiv:2103.03186 [hep-ph] .
  * Mamo and Zahed [2022] K. A. Mamo and I. Zahed, Phys. Rev. D 106, 086004 (2022), arXiv:2204.08857 [hep-ph] .
  * Sun  _et al._ [2022] P. Sun, X.-B. Tong, and F. Yuan, Phys. Rev. D 105, 054032 (2022), arXiv:2111.07034 [hep-ph] .
  * Sun  _et al._ [2021] P. Sun, X.-B. Tong, and F. Yuan, Phys. Lett. B 822, 136655 (2021), arXiv:2103.12047 [hep-ph] .
  * Collins  _et al._ [1997] J. C. Collins, L. Frankfurt, and M. Strikman, Phys. Rev. D 56, 2982 (1997), arXiv:hep-ph/9611433 .
  * Ivanov  _et al._ [2004] D. Y. Ivanov, A. Schafer, L. Szymanowski, and G. Krasnikov, Eur. Phys. J. C 34, 297 (2004), [Erratum: Eur.Phys.J.C 75, 75 (2015)], arXiv:hep-ph/0401131 .
  * Pefkou  _et al._ [2022] D. A. Pefkou, D. C. Hackett, and P. E. Shanahan, Phys. Rev. D 105, 054509 (2022), arXiv:2107.10368 [hep-lat] .
  * Ji [1997] X.-D. Ji, Phys. Rev. Lett. 78, 610 (1997), arXiv:hep-ph/9603249 .
  * Diehl [2003] M. Diehl, Phys. Rept. 388, 41 (2003), arXiv:hep-ph/0307382 .
  * Hou  _et al._ [2021] T.-J. Hou _et al._ , Phys. Rev. D 103, 014013 (2021), arXiv:1912.10053 [hep-ph] .
  * Dembinski and et al. [2020] H. Dembinski and P. O. et al. 10.5281/zenodo.3949207 (2020).
  * James and Roos [1975] F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).
  * Shanahan and Detmold [2019] P. E. Shanahan and W. Detmold, Phys. Rev. D 99, 014511 (2019), arXiv:1810.04626 [hep-lat] .
  * Ji [2021] X. Ji, Front. Phys. (Beijing) 16, 64601 (2021), arXiv:2102.07830 [hep-ph] .


# Quasinormal resonances of rapidly-spinning Kerr black holes and the
universal relaxation bound

Shahar Hod  The Ruppin Academic Center, Emeq Hefer 40250, Israel  The Hadassah
Institute, Jerusalem 91010, Israel

###### Abstract

The universal relaxation bound suggests that the relaxation times of perturbed
thermodynamical systems is bounded from below by the simple time-times-
temperature (TTT) quantum relation $\tau\times T\geq{{\hbar}\over{\pi}}$. It
is known that some perturbation modes of near-extremal Kerr black holes in the
regime $MT_{\text{BH}}/\hbar\ll m^{-2}$ are characterized by normalized
relaxation times $\pi\tau\times T_{\text{BH}}/\hbar$ which, in the approach to
the limit $MT_{\text{BH}}/\hbar\to 0$, make infinitely many oscillations with
a tiny constant amplitude around $1$ and therefore cannot be used directly to
verify the validity of the TTT bound in the entire parameter space of the
black-hole spacetime (Here $\\{T_{\text{BH}},M\\}$ are respectively the
Bekenstein-Hawking temperature and the mass of the black hole, and $m$ is the
azimuthal harmonic index of the linearized perturbation mode). In the present
compact paper we explicitly prove that all rapidly-spinning Kerr black holes
respect the TTT relaxation bound. In particular, using analytical techniques,
it is proved that all black-hole perturbation modes in the complementary
regime $m^{-1}\ll MT_{\text{BH}}/\hbar\ll 1$ are characterized by relaxation
times with the simple dimensionless property $\pi\tau\times
T_{\text{BH}}/\hbar\geq 1$.

##  I Introduction

Perturbed black-hole spacetimes are usually characterized by a relaxation
phase which is dominated by exponentially decaying quasinormal oscillations
whose discrete complex frequencies $\\{\omega_{n}\\}^{n=\infty}_{n=0}$ encode
valuable information about the fundamental physical parameters (in particular,
the mass and angular momentum) of the perturbed black hole Nollert1 ; Press ;
Cruz ; Vish ; Davis ; QNMs .

The damped quasinormal oscillations of matter and radiation fields in
perturbed black-hole spacetimes, which describe waves with outgoing boundary
conditions at spatial infinity and ingoing boundary conditions at the
absorbing black-hole horizon Detwe , reflect the fact that most perturbations
fields (fields that conform to existing no-hair theorems Whee ; Car ; Bek1 ;
Chas ; BekVec ; Hart ; Hod11 ) are eventually scattered away to infinity or
absorbed by the central black hole Notetails ; Tails ; Notescc ; Hodrc ; HerR
.

The timescale $\tau_{\text{relax}}$ associated with the relaxation dynamics of
a newly born black hole may be determined by the imaginary part of the
fundamental (least damped) quasinormal resonant frequency:

$\tau_{\text{relax}}\equiv 1/\Im\omega_{0}\ .$ (1)

Taking cognizance of the fact that black holes have thermodynamic and quantum
properties Bektem ; Hawtem , a physically important question naturally arises:
How short can the relaxation time (1) of a newly born black hole be?

A mathematically compact answer to this physically intriguing question, which
is based on standard ideas from classical thermodynamics and quantum
information theory, has been suggested in HodTTT ; Notetl :

$\tau_{\text{relax}}\times T\geq{{\hbar}\over{\pi k_{\text{B}}}}\ ,$ (2)

where $T$ is the characteristic temperature of the thermodynamic physical
system.

Intriguingly, remembering that the semi-classical Bekenstein-Hawking
temperature Bektem ; Hawtem ; Noteunit

$T_{\text{BH}}={{\kappa}\over{2\pi k_{\text{B}}}}\cdot\hbar\ $ (3)

of black holes depends linearly on the quantum Planck constant (here $\kappa$
is the characteristic surface gravity of the black-hole horizon), one realizes
that the time-times-temperature (TTT) quantum bound (2) provides a classical
lower bound Notehbarcal  on the characteristic relaxation time of a newly born
black hole. In particular, taking cognizance of Eqs. (1), (2), and (3), one
realizes that all dynamically formed black holes should be characterized by
(at least) one relaxation mode with the simple classical property

$\Im\omega_{0}\leq{{\pi T_{\text{BH}}}\over{\hbar}}={1\over 2}\kappa\ .$ (4)

##  II Relaxation dynamics of rapidly spinning Kerr black holes

Interestingly, the upper bound (4) implies that newly born near-extremal black
holes (with the property $T_{\text{BH}}\to 0$) should be characterized by
extremely long relaxation times HodTTT ; Hodan ; Als . In particular, using
the physically important and mathematically elegant Detweiler equation Detw ,
which characterizes the complex resonant spectra of rapidly-spinning (near-
extremal) Kerr black holes, one obtains the simple functional relation Hoddet
; Gruzdet ; Notehbar

$\Im\omega_{0}=\pi
T_{\text{BH}}\cdot\big{\\{}1+C\cdot\sin[2\delta\ln(MT_{\text{BH}})]\big{\\}}\cdot[1+O(MT_{\text{BH}})]\
$ (5)

for the co-rotating perturbation modes of near-extremal black holes, where
$C(l,m)$ is a constant NoteCC , $M$ is the black-hole mass, and the physical
parameter $\delta(l,m)\in\mathbb{R}$ is closely related to the characteristic
eigenvalue of the angular Teukolsky equation Teuk ; Notedelta .

Before proceeding, it is important to stress the fact that the validity of the
Detweiler resonance equation Detw , and thus also the validity of the
analytically derived expression (5) Hoddet ; Gruzdet , are restricted to the
dimensionless regime Teuk ; Starrv ; Hodmmm

$m^{2}\cdot MT_{\text{BH}}\ll 1\ ,$ (6)

where $m$ is the azimuthal harmonic index of the black-hole perturbation mode.

From Eq. (5) one immediately finds that perturbation modes of rapidly-spinning
Kerr black holes in the regime (6) are characterized by normalized relaxation
times $\pi\tau\times T_{\text{BH}}$ that oscillate infinitely many times
around $1$ as the extremal limit $MT_{\text{BH}}\to 0$ is approached. One
therefore deduces that, in the near-extremal $MT_{\text{BH}}\ll 1$ regime,
there are finite intervals of the black-hole temperature for which a
perturbation mode in the dimensionless regime (6) cannot be used to prove the
general validity of the TTT relaxation bound (2) [or equivalently, the upper
bound (4)] in black-hole physics Notenng .

The main goal of the present compact paper is to prove the physically
important fact that all rapidly-spinning Kerr black holes respect the TTT
relaxation bound (2). To this end, we shall study below the linearized
relaxation dynamics of newly born near-extremal Kerr black holes. In
particular, using analytical techniques, we shall explicitly prove that
composed black-hole-field perturbation modes in the eikonal large-$m$ regime

$m^{-1}\ll MT_{\text{BH}}\ll 1\ $ (7)

are characterized by relaxation times that respect the TTT relaxation bound.

##  III Resonant relaxation spectra of rapidly spinning Kerr black holes in
the eikonal large-frequency regime

Using a physically intuitive and mathematically elegant analysis, Mashhoon
Mash  has presented a simple analytical technique for calculating the
quasinormal resonance spectra of black holes in the large-frequency
(geometric-optics) regime. In particular, the physical idea presented in Mash
; Goeb  is to relate, in the eikonal large-frequency regime, the real part of
the black-hole quasinormal frequencies to the angular velocity which
characterizes the motion of massless particles along the equatorial null
circular geodesic of the black-hole spacetime and to relate the imaginary part
of the black-hole quasinormal frequencies to the instability timescale which
characterizes the gradual leakage of the perturbed massless particles (null
rays) from the unstable null circular geodesic of the black-hole spacetime.

In particular, using a perturbation scheme for the instability which
characterizes the dynamics of equatorial null circular geodesics in the Kerr
black-hole spacetime, Mashhoon has provided the remarkably compact formula
Mash

$\omega_{n}=m\omega_{+}-i(n+{1\over 2})\beta\omega_{+}\ \ \ \ ;\ \ \ \
n=0,1,2,...\ ,$ (8)

for the discrete quasinormal resonant spectra of spinning Kerr black holes. As
emphasized in Mash , this resonance formula is valid for perturbation modes in
the eikonal large-frequency regime

$l=m\gg 1\ ,$ (9)

where $\\{l,m\\}$ are the (spheroidal and azimuthal) angular harmonic indexes
which characterize the linearized field mode.

The various terms in the resonance formula (8) have the following physical
interpretations Mash :
(1) The term

$r_{\text{ph}}(a/M)\equiv 2M\cdot\Big{\\{}1+\cos\Big{[}{2\over
3}\arccos\Big{(}-{{a}\over{M}}\Big{)}\Big{]}\Big{\\}}\ $ (10)

is the radius of the co-rotating equatorial null circular geodesic, where
$\\{M.a\\}$ are respectively the mass and angular momentum per unit mass
Noteaaa  of the central spinning Kerr black hole.
(2) The term

$\omega_{+}(a/M)\equiv{{M^{1/2}\over{r^{3/2}_{\text{\text{ph}}}+aM^{1/2}}}}\ $
(11)

is the Kepler frequency which characterizes the co-rotating equatorial null
circular geodesic of the black-hole spacetime.
(3) The dimensionless functional expression of $\beta=\beta(M,a)$ is given by
Mash

$\beta(a/M)={{(12M)^{1/2}(r_{\text{ph}}-r_{+})(r_{\text{ph}}-r_{-})}\over{r^{3/2}_{\text{ph}}(r_{\text{ph}}-M)}}\
,$ (12)

where

$r_{\pm}=M\pm(M^{2}-a^{2})^{1/2}\ $ (13)

are the (outer and inner) horizon radii of the spinning Kerr black hole.

We shall henceforth focus our attention on rapidly-spinning (near-extremal)
Kerr black holes, which are characterized by the simple dimensionless relation

${{a}\over{M}}=1-\epsilon\ \ \ \text{with}\ \ \ 0\leq\epsilon\ll 1\ .$ (14)

From Eqs. (13) and (14) one finds the relation

${{r_{\pm}}\over{M}}=1\pm\sqrt{2}\sqrt{\epsilon}\mp{{\epsilon^{3/2}}\over{2\sqrt{2}}}+O(\epsilon^{5/2})\
.$ (15)

Substituting (14) into Eq. (10), one obtains the near-extremal
(small-$\epsilon$) dimensionless functional expression

${{r_{\text{ph}}}\over{M}}=1+{{2\sqrt{2}}\over{\sqrt{3}}}\cdot\sqrt{\epsilon}+{4\over
9}\cdot\epsilon+O(\epsilon^{3/2})\ $ (16)

for the radius of the co-rotating equatorial null circular geodesic.
Substituting Eqs. (14), (15), and (16) into Eqs. (11) and (12), one finds

$M\omega_{+}={1\over
2}-{{\sqrt{3}}\over{2\sqrt{2}}}\cdot\sqrt{\epsilon}+{7\over
12}\cdot\epsilon+O(\epsilon^{3/2})\ $ (17)

and

$\beta=\sqrt{2}\cdot\sqrt{\epsilon}-{{4}\over{3\sqrt{3}}}\cdot\epsilon+O(\epsilon^{3/2})\
.$ (18)

Substituting Eqs. (17) and (18) into the resonance formula (8), one obtains
the dimensionless functional expression

$M\omega_{n}=m\cdot\Big{(}{1\over
2}-{{\sqrt{3}}\over{2\sqrt{2}}}\cdot\sqrt{\epsilon}+{7\over
12}\cdot\epsilon\Big{)}-i\big{(}n+{1\over
2}\big{)}\cdot\Big{(}{{\sqrt{\epsilon}}\over{\sqrt{2}}}-{{13}\over{6\sqrt{3}}}\cdot\epsilon\Big{)}+O(\epsilon^{3/2})\
.$ (19)

for the complex quasinormal resonant frequencies of rapidly-spinning (near-
extremal) Kerr black holes in the eikonal large-frequency regime $l=m>>1$ [see
Eq. (9)].

It is physically interesting and mathematically convenient to express the
analytically derived Kerr resonance spectrum (19) in terms of the black-hole
Bekenstein-Hawking temperature Bektem ; Hawtem

$T_{\text{BH}}={{r_{+}-r_{-}}\over{4\pi(r^{2}_{+}+a^{2})}}={{1}\over{\pi
M}}\cdot\Big{[}{{\sqrt{\epsilon}}\over{2\sqrt{2}}}-{{\epsilon}\over{2}}+O(\epsilon^{3/2})\Big{]}\
$ (20)

and the black-hole angular velocity

$\Omega_{\text{H}}={a\over{r^{2}_{+}+a^{2}}}={{1}\over{M}}\cdot\Big{[}{1\over
2}-{{\sqrt{\epsilon}}\over{\sqrt{2}}}+{{\epsilon}\over{2}}+O(\epsilon^{3/2})\Big{]}\
.$ (21)

Substituting (20) and (21) into Eq. (19), one finds the black-hole resonance
formula

$\omega_{n}=\Big{\\{}m\Omega_{\text{H}}\cdot\big{[}1+\pi(4-2\sqrt{3})MT_{\text{BH}}\big{]}-i\big{(}n+{1\over
2}\big{)}\cdot 2\pi
T_{\text{BH}}\cdot\big{[}1-{{2\pi}\over{9}}(13\sqrt{3}-18)MT_{\text{BH}}\big{]}\Big{\\}}\cdot[1+O(MT_{\text{BH}})]\
.$ (22)

##  IV Summary

Motivated by the suggested time-times-temperature (TTT) universal relaxation
bound (2), we have studied, using analytical techniques, the quasinormal
resonance spectra of near-extremal (rapidly-spinning) Kerr black holes.

It was first noted that perturbation modes of near-extremal Kerr black holes
in the regime $m^{2}\cdot MT_{\text{BH}}\ll 1$ are characterized by normalized
relaxation times $\pi\tau\times T_{\text{BH}}$ which, in the approach to the
extremal limit $MT_{\text{BH}}\to 0$, oscillate infinitely many times around
$1$ and therefore cannot be used to verify the general validity of the
suggested TTT relaxation bound (2) [or equivalently, the upper bound (4)] in
black-hole physics.

We have therefore analyzed the black-hole relaxation spectra in the
complementary regime $m^{-1}\ll MT_{\text{BH}}\ll 1$. In particular, we have
explicitly proved that the equatorial black-hole-field perturbation modes in
the eikonal large-frequency regime (7) are characterized by the simple near-
extremal relation

$\Im\omega_{0}=\pi
T_{\text{BH}}\cdot\big{[}1-{{2\pi}\over{9}}(13\sqrt{3}-18)MT_{\text{BH}}\big{]}<\pi
T_{\text{BH}}\ \ \ \ \text{for}\ \ \ \ m^{-1}\ll MT_{\text{BH}}\ll 1\ .$ (23)

[It is worth noting that the functional relation (23) characterizes the
fundamental ($n=0$) resonant mode of the analytically derived black-hole
relaxation spectrum (22)]. We therefore conclude that the eikonal relaxation
modes (23), which characterize the composed black-hole-field system in the
dimensionless large-frequency regime $m^{-1}\ll MT_{\text{BH}}\ll 1$,
guarantee that newly born near-extremal Kerr black holes respect the suggested
TTT universal relaxation bound (2) Notenln .

ACKNOWLEDGMENTS

This research is supported by the Carmel Science Foundation. I would like to
thank Yael Oren, Arbel M. Ongo, Ayelet B. Lata, and Alona B. Tea for helpful
discussions.

## References

  * (1) H. P. Nollert, Class. Quantum Grav. 16, R159 (1999).
  * (2) W. H. Press, Astrophys. J. 170, L105 (1971).
  * (3) V. de la Cruz, J. E. Chase and W. Israel, Phys. Rev. Lett. 24, 423 (1970).
  * (4) C.V. Vishveshwara, Nature 227, 936 (1970).
  * (5) M. Davis, R. Ruffini, W. H. Press and R. H. Price, Phys. Rev. Lett. 27, 1466 (1971).
  * (6) E. W. Leaver, Proc. R. Soc. A 402, 285 (1985); B. Mashhoon, Phys. Rev. D 31, 290 (1985); H. P. Nollert, Phys. Rev. D 47, 5253 (1993); S. Hod, Phys. Rev. Lett. 81, 4293 (1998) [arXiv:gr-qc/9812002]; S. Hod, Phys. Rev. D 67, 081501 (2003) [arXiv:gr-qc/0301122]; S. Hod and U. Keshet, Class. Quant. Grav. 22, L71 (2005) [arXiv:gr-qc/0505112]; S. Hod and U. Keshet, Phys.Rev. D 73, 024003 (2006) [arXiv:hep-th/0506214]; U. Keshet and S. Hod, Phys. Rev. D 76, R061501 (2007) [arXiv:0705.1179]; Y. Décanini, A. Folacci, and B. Raffaelli, Phys. Rev. D 84, 084035 (2011); S. Hod, Phys. Rev. D 88, 084018 (2013) [arXiv:1311.3007].
  * (7) S. L. Detweiler, in Sources of Gravitational Radiation, edited by L. Smarr (Cambridge University Press, Cambridge, England, 1979).
  * (8) R. Ruffini and J. A. Wheeler, Phys. Today 24, 30 (1971).
  * (9) B. Carter, in Black Holes, Proceedings of 1972 Session of Ecole d’ete de Physique Theorique, edited by C. De Witt and B. S. De Witt (Gordon and Breach, New York, 1973).
  * (10) A. E. Mayo and J. D. Bekenstein, Phys. Rev. D 54, 5059 (1996).
  * (11) J. E. Chase, Commun. Math. Phys. 19, 276 (1970); J. D. Bekenstein, Phys. Rev. Lett. 28, 452 (1972); C. Teitelboim, Lett. Nuovo Cimento 3, 326 (1972); I. Pena and D. Sudarsky, Class. Quant. Grav. 14, 3131 (1997).
  * (12) J. D. Bekenstein, Phys. Rev. D 5, 1239 (1972); 5, 2403 (1972); M. Heusler, J. Math. Phys. 33, 3497 (1992); D. Sudarsky, Class. Quant. Grav. 12, 579 (1995).
  * (13) J. Hartle, Phys. Rev. D 3, 2938 (1971); C. Teitelboim, Lett. Nuovo Cimento 3, 397 (1972).
  * (14) S. Hod, Phys. Rev. D 84, 124030 (2011) [arXiv:1112.3286].
  * (15) It is interesting to note that some non-decaying matter fields can be supported in asymptotically flat black-hole spacetimes without being absorbed into the central supporting black holes Hodrc ; HerR .
  * (16) S. Hod, Phys. Rev. D 86, 104026 (2012) [arXiv:1211.3202]; S. Hod, The Euro. Phys. Journal C 73, 2378 (2013) [arXiv:1311.5298]; S. Hod, Phys. Rev. D 90, 024051 (2014) [arXiv:1406.1179]; S. Hod, Phys. Lett. B 739, 196 (2014) [arXiv:1411.2609]; S. Hod and O. Hod, Phys. Rev. D 81, 061502 Rapid communication (2010) [arXiv:0910.0734]; S. Hod, Phys. Lett. B 708, 320 (2012) [arXiv:1205.1872]; S. Hod, Phys. Lett. B 751, 177 (2015); S. Hod, Class. and Quant. Grav. 33, 114001 (2016) [arXiv:1705.08905]; S. Hod, J. High Energy Phys. 01, 30 (2017) [arXiv:1612.00014].
  * (17) C. A. R. Herdeiro and E. Radu, Phys. Rev. Lett. 112, 221101 (2014); C. A. R. Herdeiro and E. Radu, Phys. Rev. D 89, 124018 (2014); C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D 23, 1442014 (2014); Y. Brihaye, C. Herdeiro, and E. Radu, Phys. Lett. B 739, 1 (2014); C. L. Benone, L. C. B. Crispino, C. Herdeiro, and E. Radu, Phys. Rev. D 90, 104024 (2014); J. C. Degollado and C. A. R. Herdeiro, Phys. Rev. D 90, 065019 (2014); C. Herdeiro, E. Radu, and H. Rúnarsson, Phys. Lett. B 739, 302 (2014); C. Herdeiro and E. Radu, Class. Quantum Grav. 32 144001 (2015); C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D 24, 1542014 (2015); C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D 24, 1544022 (2015); J. C. Degollado and C. A. R. Herdeiro, Gen. Rel. Grav. 45, 2483 (2013); P. V. P. Cunha, C. A. R. Herdeiro, E. Radu, and H. F. Rúnarsson, Phys. Rev. Lett. 115, 211102 (2015); B. Kleihaus, J. Kunz, and S. Yazadjiev, Phys. Lett. B 744, 406 (2015); C. A. R. Herdeiro, E. Radu, and H. F. Rúnarsson, Phys. Rev. D 92, 084059 (2015); C. Herdeiro, J. Kunz, E. Radu, and B. Subagyo, Phys. Lett. B 748, 30 (2015); C. A. R. Herdeiro, E. Radu, and H. F. Rúnarsson, Class. Quant. Grav. 33, 154001 (2016); C. A. R. Herdeiro, E. Radu, and H. F. Rúnarsson, arXiv:1604.06202; Y. Brihaye, C. Herdeiro, and E. Radu, Phys. Lett. B 760, 279 (2016); Y. Ni, M. Zhou, A. C. Avendano, C. Bambi, C. A. R. Herdeiro, and E. Radu, JCAP 1607, 049 (2016); M. Wang, arXiv:1606.00811; J. F. M. Delgado, C. A. R. Herdeiro, and E. Radu, Phys. Rev. D 94, 024006 (2016); J. F. M. Delgado, C. A. R. Herdeiro, E. Radu, H. Rúnarsson, Phys. Lett. B 761, 234 (2016).
  * (18) It is worth mentioning that the late-time relaxation dynamics of perturbed black holes in asymptotically flat spacetimes is also characterized by inverse power-law decaying tails which reflect the late-time backscattering of dynamical perturbation fields by the small curvature potential at asymptotically far regions Tails .
  * (19) R. H. Price, Phys. Rev. D 5, 2419 (1972); C. Gundlach, R. H. Price, and J. Pullin, Phys. Rev. D 49, 883 (1994); J. Bicák, Gen. Relativ. Gravitation 3, 331 (1972); E. S. C. Ching, P. T. Leung, W. M. Suen, and K. Young, Phys. Rev. Lett. 74, 2414 (1995); E. S. C. Ching, P. T. Leung, W. M. Suen, and K. Young, Phys. Rev. D 52, 2118 (1995); S. Hod and T. Piran, Phys. Rev. D 58, 044018 (1998) [arXiv:gr-qc/9801059]; S. Hod, Phys. Rev. D 58, 104022 (1998) [arXiv:gr-qc/9811032]; S. Hod, Phys. Rev. D 61, 024033 (2000) [arXiv:gr-qc/9902072]; S. Hod, Phys. Rev. D 61, 064018 (2000) [arXiv:gr-qc/9902073]; L. Barack, Phys. Rev. D 61, 024026 (2000); S. Hod, Phys. Rev. Lett. 84, 10 (2000) [arXiv:gr-qc/9907096]; S. Hod, Phys. Rev. D 60, 104053 (1999) [arXiv:gr-qc/9907044]; S. Hod, Class. Quant. Grav. 26, 028001 (2009) [arXiv:0902.0237]; S. Hod, Class. Quant. Grav. 18, 1311 (2001) [arXiv:gr-qc/0008001]; S. Hod, Phys. Rev. D 66, 024001 (2002) [arXiv:gr-qc/0201017]; R. J. Gleiser, R. H. Price, and J. Pullin, Class. Quantum Grav. 25, 072001 (2008); M. Tiglio, L. E. Kidder, and S. A. Teukolsky, Class. Quant. Grav. 25, 105022 (2008); R. Moderski and M. Rogatko, Phys. Rev. D 77, 124007 (2008); X. He and J. Jing, Nucl. Phys.B 755, 313 (2006); H. Koyama and A. Tomimatsu, Phys. Rev. D 65, 084031 (2002); B. Wang, C. Molina, and E. Abdalla, Phys. Rev. D 63, 084001 (2001); A. Zenginoglu and M. Tiglio, Phys. Rev.D 80, 024044 (2009); S. Hod, Class. and Quant. Grav. 30, 237002 (2013) [arXiv:1402.4819]; S. Hod, J. High Energy Phys. 1309, 056 (2013) [arXiv:1310.4247].
  * (20) J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973).
  * (21) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975).
  * (22) S. Hod, Phys. Rev. D 75, 064013 (2007) [arXiv:0611004].
  * (23) It is interesting to note that the quantum relaxation bound (2) provides a quantitative formulation for the fundamental third law of thermodynamics.
  * (24) We use natural units in which $G=c=k_{\text{B}}=1$.
  * (25) Substituting the semi-classical black-hole relation (3) into the quantum bound (2), one finds that the lower bound on the relaxation timescale of a perturbed black hole is $\hbar$-independent [see Eq. (4) below].
  * (26) S. Hod, Class. Quant. Grav. 24, 4235 (2007) [arXiv:0705.2306].
  * (27) A. Pesci, Class. Quant. Grav. 24, 6219 (2007); S. Hod, Phys. Lett. B 666, 483 (2008) [arXiv:0810.5419]; S. Hod, Phys. Rev. D 80, 064004 (2009) [arXiv:0909.0314]; S. Hod, Phys. Rev. D. 84, 044046 (2011) [arXiv:1109.4080]; S. Hod, Phys. Lett. B 715, 348 (2012) [arXiv:1207.5282]; S. Hod, The Euro. Phys. Jour. C (Letter) 75, 272 (2015) [arXiv:1410.2252].
  * (28) S. Detweiler, Astrophys. J. 239, 292 (1980).
  * (29) A. Gruzinov, arXiv:0705.1725.
  * (30) S. Hod, Phys. Rev. D 78, 084035 (2008) [arXiv:0811.3806].
  * (31) For brevity, we shall henceforth use units in which $\hbar=1$.
  * (32) The amplitude $C=C(l,m)$ depends on the angular harmonic indexes of the composed black-hole-field perturbation mode Gruzdet ; Hoddet . It is worth noting that one finds $|C|\ll 1$ for most perturbation modes Gruzdet ; Hoddet .
  * (33) S. A. Teukolsky, Phys. Rev. Lett. 29, 1114 (1972); S. A. Teukolsky, Astrophys. J. 185, 635 (1973); S. A. Teukolsky and W. H. Press, Astrophys. J. 193, 443 (1974); W. H. Press and S. A. Teukolsky, Astrophys. J. 185, 649 (1973).
  * (34) It is worth noting that the simple relation (5) characterizes all near-extremal perturbation modes of rapidly-spinning Kerr black holes with the property $\delta^{2}>0$ Gruzdet ; Hoddet . For near-extremal perturbation modes with the property $\delta^{2}<0$ one finds the relation $\Im\omega_{0}=\pi T_{\text{BH}}\cdot(1+2|\delta|)\cdot[1+O(MT_{\text{BH}})]>\pi T_{\text{BH}}$ Gruzdet ; Hoddet .
  * (35) A. A. Starobinsky, Zh. Eksp. Teor. Fiz. 64, 48 (1973) [Sov. Phys. JETP 37, 28 (1973)]; A. A. Starobinsky and S. M. Churilov, Zh. Eksp. Teor. Fiz. 65, 3 (1973) [Sov. Phys. JETP 38, 1 (1973)].
  * (36) S. Hod, Class. and Quant. Grav. 32, 134002 (2015) [arXiv:1607.00003]; S. Hod, Phys. Rev. D 94, 044036 (2016) [arXiv:1609.07146].
  * (37) It is interesting to note that numerical studies indicate that, for a given black-hole temperature, there is at least one black-hole-field perturbation mode in the regime (6) that conforms to the TTT relaxation bound (4) HodTTT ; Gruzdet ; Hoddet . However, to the best of our knowledge there is no general argument which guarantees that all black holes must be characterized by this property and it is therefore physically important to present a more generic mathematical proof for the validity of the suggested TTT relaxation bound (2) [or equivalently, the upper bound (4)] in black-hole physics.
  * (38) B. Mashhoon, Phys. Rev. D 31, 290 (1985).
  * (39) C. J. Goebel, Astrophys. J. 172, L95 (1972).
  * (40) We shall assume $a\geq 0$ without loss of generality.
  * (41) It is worth noting again that one learns from the functional relation (5) that some linearized perturbation modes in the dimensionless regime (6) may not conform to the TTT relaxation bound (2) for particular black-hole temperatures. The TTT relaxation bound asserts that a dynamically formed black hole should be characterized by (at least) one relaxation mode with the property (4). In the present compact paper we have explicitly proved that the resonant perturbation modes (23) of near-extremal Kerr black holes in the dimensionless regime (7) conform to the suggested TTT relaxation bound (2) [or equivalently, to the upper bound (4)].