for any $\widetilde{w} \in \widetilde{\WW}$. Moreover, if $\left(\frac{r}{2 \mu \sqrt{d}}\right) \mathcal{A} = \widetilde{\mathcal{A}}: \widetilde{X}^n \to \widetilde{\WW}$ is any algorithm and $\widetilde{X} \sim \widetilde{\mathcal{D}}^n$, then by <ref> and <ref>, we have
\[
\expec \widetilde{F}(\widetilde{\mathcal{A}}(\widetilde{X})) - \widetilde{F}^* = \frac{\mu}{2}\expec \|\widetilde{\mathcal{A}}(\widetilde{X}) - \expec \widetilde{\mathcal{D}} \|^2 = \frac{\mu}{2}\left(\frac{r}{\mu \sqrt{5d}}\right)^2 \expec \|\mathcal{A}(X) -  \expec \mathcal{D}\|^2 \geq \frac{r^2}{1280 \mu n}.
\]

Next, we prove the second (private) term in the lower bound. Let $f$ be as defined above. For our unscaled hard distribution, we follow [10, 45] and define a family of distributions $\{Q_\nu\}_{\nu \in \mathcal{V}}$ on $\mathbb{R}^d$, where $\mathcal{V} \subset \{\pm 1\}^d$ will be defined later. For any given $\nu \in \mathcal{V}$, we define the distribution $Q_\nu$ as follows: $X_\nu \sim Q_\nu$ iff
\begin{equation}
\label{eq:harddist}
X_\nu =
\begin{cases}
0 &\mbox{with probability $1 - p$} \\
p^{-1/k} \nu &\mbox{with probability $p$}
\end{cases}
\end{equation}

where $p := \min\left(1, \frac{\sqrt{d}}{n \sqrt{\rho}}\right)$. Now, we select a set $\mathcal{V} \subset \{\pm 1\}^d$ such that $|\mathcal{V}| \geq 2^{d/20}$ and $d_{\text{Ham}}(\nu, \nu') \geq \frac{d}{8}$ for all $\nu, \nu' \in \mathcal{V}, \nu \neq \nu'$: such $\mathcal{V}$ exists by standard Gilbert-Varshamov bound (see e.g. <cit.>). For any $\nu \in \mathcal{V}$, if $x \sim Q_\nu$ and $w \in \WW := B_2^d(0,\sqrt{d} p^{-1/k})$, then
\[
\expec \| \nabla f(w,x) \|^k = \expec \| w - x \|^k = \|w\|^k(1-p) + \|w - p^{-1/k} \nu\|^k p \leq 1 + p[p^{-1}(4d)^{k/2}] \leq 2(4d)^{k/2},
\]

since $\|w - p^{-1/k} \nu\|^2 \leq 4dp^{-2/k}$ for any $w \in B_2^d(0,\sqrt{d} p^{-1/k}), \nu \in \{\pm 1\}^d$.
Note also that our choice of $\WW$ and $p \leq 1$ ensures that $\expec[Q_\nu] \in \WW$.
Moreover, as in the proof of <cit.>, zCDP Fano's inequality (see <cit.>) implies that for any $\rho$-zCDP algorithm $\mathcal{A}$,
\begin{equation}
\label{eq:tt}
\sup_{\nu \in \mathcal{V}} \expec_{X \sim Q_\nu^n, \mathcal{A}} \|\mathcal{A}(X) - \expec Q_\nu\|^2 = \Omega\left(d \min\left\{1, \left(\frac{\sqrt{d}}{n\sqrt{\rho}}\right)^{\frac{2k-2}{k}}\right\} \right).
\end{equation}

\[
\expec_{X \sim Q_\nu^n, \mathcal{A}} F(\mathcal{A}(X)) - F^* = \Omega\left(d \min\left\{1, \left(\frac{\sqrt{d}}{n\sqrt{\rho}}\right)^{\frac{2k-2}{k}}\right\} \right)
\]

for some $\nu \in \mathcal{V}$, by <ref>. Now we scale our hard instance: $f \to \widetilde{f} = \mu f$, $\WW \to \widetilde{\WW} := \frac{r}{4 \mu \sqrt{d}} \WW$, $\XX \to \widetilde{\XX} := \frac{r}{4 \mu \sqrt{d}}\XX$ and $\mathcal{D} \to \widetilde{\mathcal{D}} = \frac{r}{4 \mu \sqrt{d}} \mathcal{D}$. Then $\wt{f}(\cdot, \widetilde{x})$ is $\mu$-strongly convex and $\mu$-smooth, and
\[
\expec \| \nabla \widetilde{f}(\widetilde{w}, \widetilde{x}) \|^k = \expec \| \widetilde{w} - \widetilde{x} \|^k = \mu^k \left(\frac{r}{4 \mu \sqrt{d}}\right)^k \expec \| \nabla f(w,x) \|^k \leq r^k,
\]

for any $\widetilde{w} \in \widetilde{\WW}$. Moreover, if $\left(\frac{r}{4 \mu \sqrt{d}}\right) \mathcal{A} = \widetilde{\mathcal{A}}: \widetilde{X}^n \to \widetilde{\WW}$ is any algorithm and $\widetilde{X} \sim \widetilde{\mathcal{D}}^n$, then
\[
\expec \widetilde{F}(\widetilde{\mathcal{A}}(\widetilde{X})) - \widetilde{F}^* = \frac{\mu}{2}\expec \|\widetilde{\mathcal{A}}(\widetilde{X}) - \expec \widetilde{\mathcal{D}} \|^2 = \frac{\mu}{2}\left(\frac{r}{4 \mu \sqrt{d}}\right)^2 \expec \|\mathcal{A}(X) -  \expec \mathcal{D}\|^2 \geq \frac{r^2}{32 \mu d}\Omega\left(d \min\left\{1, \left(\frac{\sqrt{d}}{n\sqrt{\rho}}\right)^{\frac{2k-2}{k}}\right\} \right),
\]

by <ref>. This completes the proof.

§ DETAILS AND PROOFS FOR <REF>

§.§ Noisy Clipped Accelerated SGD for Smooth Convex Losses (<ref>)

We first present <ref>, which is a generalized version of <ref> that allows for any  and will be useful for our analysis.
Framework for DP Accelerated Stochastic Approximation (AC-SA)

[1]
Data $X \in \XX^n$, iteration number $T \leq n$,

stepsize parameters $\{\eta_t \}_{t \in [T]}, \{\alpha_t \}_{t \in [T]}$ with $\alpha_1 = 1, \alpha_t \in (0,1)$ $\forall t \geq 2$,

DP $\texttt{MeanOracle}$.

Initialize $w_0^{ag} = w_0 \in \WW$ and $t = 1$.
$t \in [T]$
$w_t^{md} :=
%\frac{(1- \alpha_t)(\mu + \eta_t)}{\eta_t + (1 - \alpha_t^2)\mu}w_{t-1}^{ag} + \frac{\alpha_t[(1-\alpha_t)\mu + \eta_t]}{\eta_t + (1-\alpha_t^2)\mu}w_{t-1}
(1- \alpha_t)w_{t-1}^{ag} + \alpha_t w_{t-1}$.
Draw new batch $\mathcal{B}_t$ (without replacement) of $n/T$ samples from $X$.
$\tilt(w_t^{md}) := \texttt{MeanOracle}(\{\nabla f(w_t^{md}, x)\}_{x \in \mathcal{B}_t}; \frac{n}{T}; \frac{\epsilon^2}{2})$
$w_{t} :=
%\argmin_{w \in \mathcal{W}}\left\{\alpha_t \left[\langle \tilt(w_t^{md}), w\rangle + \frac{\mu}{2}\|w_t^{md} - w\|^2 \right] + \left[(1-\alpha_t) \frac{\mu}{2} + \frac{\eta_t}{2}\right]\|w_{t-1} - w\|^2\right\}.
\argmin_{w \in \mathcal{W}}\left\{\alpha_t\langle \tilt(w_t^{md}), w\rangle + \frac{\eta_t}{2}\|w_{t-1} - w\|^2\right\}.
$w_{t}^{ag} := \alpha_t w_t + (1-\alpha_t)w_{t-1}^{ag}.$

Output: $w_T^{ag}$.

Proposition <ref> provides excess risk guarantees for <ref> in terms of the bias and variance of the .

Consider <ref> run with a  satisfying $\tilt(w_t^{md}) = \nabla F(w_t^{md}) + b_t + N_t$, where $\|b_t\| \leq B$ (with probability $1$), $\expec N_t = 0$, $\expec\|N_t\|^2 \leq \Sigma^2$ for all $t \in [T-1]$, and $\{N_t\}_{t=1}^T$ are independent. Assume that $F: \WW \to \mathbb{R}$ is convex and  $\beta$-smooth, $F(w_0) - F^* \leq \Delta$, and $\|w_0 - \ws\| \leq D$. Suppose parameters are chosen in <ref> so that for all $t \in [T]$, $\eta_t > \beta \alpha_t^2$ and $\eta_t/\Gamma_t = \eta_1/\Gamma_1$, where
\[
\Gamma_t :=
\begin{cases}
1, &\mbox{$t = 1$}\\
(1 - \alpha_t) \Gamma_t, &\mbox{$t \geq 2$}.
\end{cases}
\]

\[
\expec F(w_T^{ag}) - F^* \leq  \frac{\Gamma_T \eta_1 D^2}{2} + \Gamma_T \sum_{t=1}^T\left[\frac{2 \alpha_t^2(\Sigma^2 + B^2)}{\Gamma_t(\eta_t - \beta \alpha_t^2)} + \frac{\alpha_t}{\Gamma_t} BD\right].
\]

In particular, choosing $\alpha_t = \frac{2}{t + 1}$ and $\eta_t = \frac{4\eta}{t(t+1)}, ~\forall t \geq 1$, where $\eta \geq 2\beta$ implies
\[
\expec F(w_T^{ag}) - F^* \leq \frac{4 \eta D^2}{T(T+1)} + \frac{4(\Sigma^2 + B^2)(T+2)}{3\eta} + BD.
\]

Further, setting $\eta = \max\left\{ 2 \beta, \frac{T^{3/2}\sqrt{\Sigma^2 + B^2}}{D}\right\}$ implies
\begin{equation}
%\label{eq:ACSA rate T}
\expec F(w_T^{ag}) - F^* \lesssim \frac{\beta D^2}{T^2} + \frac{D(\Sigma + B)}{\sqrt{T}} + BD.
\end{equation}

We begin by extending <cit.> to biased/noisy stochastic gradients. Fix any $w_{t-1}, w_{t-1}^{ag} \in \WW$. By <cit.>, we have
\begin{equation}
\label{eq: gl lem3}
F(w_t^{ag}) \leq (1-\alpha_t) F(w_{t-1}^{ag}) + \alpha[F(z) + \langle \nabla F(z), w_t - z \rangle] + \frac{\beta}{2}\|w_t^{ag} - z\|^2,
\end{equation}

for any $z \in \WW$. Denote
\[
\Upsilon_t(w) := \alpha_t \langle N_t + b_t, w - w_{t-1} \rangle + \frac{\alpha_t^2 \|N_t + b_t\|^2}{\eta_t - \beta \alpha_t^2}
\]

and $d_t := w_t^{ag} - w_t^{md} = \alpha_t(w_t - w_{t-1})$. Then using <ref> with $z = w_t^{md}$, we have
\begin{align}
\label{eq: 3.20}
F(w_t^{ag}) &\leq (1 - \alpha_t) F(w_{t-1}^{ag}) + \alpha_t[F(w_t^{md}) + \langle \nabla F(w_t^{md}), w_t - w_t^{md} \rangle] + \frac{\beta}{2}\|d_t\|^2 \nonumber \\
&= (1 - \alpha_t) F(w_{t-1}^{ag}) + \alpha_t[F(w_t^{md}) + \langle \nabla F(w_t^{md}), w_t - w_t^{md}\rangle] + \frac{\eta_t}{2}\|w_{t-1} - w_t\|^2 - \frac{\eta_t - \beta \alpha_t^2}{2 \alpha_t^2}\|d_t\|^2,
\end{align}

by the expression for $d_t$. Now we apply <cit.> with $p(u) = \alpha_t [\langle \tilt(w_t^{md}, u \rangle], ~\mu_1 = 0, ~\mu_2 = \eta_t, ~ \tilde{x} = w_t^{md}$, and $\tilde{y} = w_{t-1}$ to obtain (conditional on all randomness) for any $w \in \WW$:

\begin{align*}
&\alpha_t[F(w_t^{md}) + \langle \tilt(w_t^{md}), w_t - w_t^{md} \rangle] + \frac{\eta_t}{2}\|w_{t-1} - w_t\|^2 \\
&\leq \alpha_t[F(w_t^{md}) + \langle \nabla F(w_t^{md}), w - w_t^{md}\rangle] \\
&\;\;\;+ \alpha_t \langle N_t + b_t, w - w_t^{md} \rangle + \frac{\eta_t}{2}\|w_{t-1} - w\|^2 - \frac{\eta_t}{2}\|w_t - w\|^2.
\end{align*}

Next, we combine the above inequality with <ref> to get
\begin{align}
\label{eq: 75}
F(w_t^{ag}) &\leq (1-\alpha_t)F(w_{t-1}^{ag}) + \alpha_t[F(w_t^{md}) + \langle \nabla F(w_t^{md}), w - w_t^{md}\rangle] + \frac{\eta_t}{2}\left[\|w_{t-1} - w\|^2 - \|w_t - w\|^2 \right] \nonumber \\
&\;\;\;+ \underbrace{-\frac{\eta_t - \beta \alpha_t^2}{2\alpha_t^2}\|d_t\|^2 + \alpha_t\langle N_t + b_t, w - w_t\rangle }_{U_t},
\end{align}

for all $w \in \WW$. By Cauchy-Schwartz, we can bound
\begin{align}
\label{eq: U_t}
U_t &\leq - \frac{\eta_t - \beta \alpha_t^2}{2 \alpha_t^2}\|d_t\|^2 + \|N_t + b_t\| \|d_t\| + \alpha_t \langle N_t + b_t, w - w_{t-1}\rangle \nonumber \\
&\leq \Upsilon_t(w),
\end{align}

where the last inequality follows from maximizing the concave quadratic function $q(\|d_t\|) := -\left[\frac{\eta_t - \beta \alpha_t^2}{2\alpha_t^2}\right]\|d_t\|^2 + \|N_t + b_t\| \|d_t\|$ with respect to $\|d_t\|$. Plugging the bound <ref> back into <ref> shows that
\begin{align}
\label{eq: 77}
F(w_t^{ag}) &\leq (1-\alpha_t)F(w_{t-1}^{ag}) + \alpha_t[F(w_t^{md}) + \langle \nabla F(w_t^{md}), w - w_t^{md}\rangle] + \frac{\eta_t}{2}\left[\|w_{t-1} - w\|^2 - \|w_t - w\|^2 \right] \nonumber \\
\end{align}

Then it can be shown (see <cit.>) that the assumptions on $\eta_t$ and $\alpha_t$ imply that
\begin{align}
\label{eq: gl prop5}
F(w_T^{ag}) - \Gamma_T \sum_{t=1}^T\left[\frac{\alpha_t}{\Gamma_t}\left(F(w_t^{md}) + \langle \nabla F(w_t^{md}), w - w_t^{md} \rangle \right) \right]
&\leq \Gamma_T \sum_{t=1}^T \frac{\eta_t}{2\Gamma_t}[\|w_{t-1} - w\|^2 - \|w_t - w\|^2] \\
&\;\;\; + \Gamma_T \sum_{t=1}^T \frac{\Upsilon_t(w)}{\Gamma_t},
\end{align}

for any $w \in \WW$ and any $T \geq 1$. Now,
\[
\sum_{t=1}^T \frac{\alpha_t}{\Gamma_t} = \frac{1}{\Gamma_T}
\]

by definition. Hence by convexity of $F$,
\[
\sum_{t=1}^T\left[\frac{\alpha_t}{\Gamma_t}\left(F(w_t^{md}) + \langle \nabla F(w_t^{md}), w - w_t^{md} \rangle \right) \right] \leq F(w), ~\forall w \in \WW.
\]

Also, since $\Gamma_t/\eta_t = \Gamma_1/\eta_1$ for all $t \geq 1$, we have
\[
\Gamma_T \sum_{t=1}^T \frac{\eta_t}{2\Gamma_t}[\|w_{t-1} - w\|^2 - \|w_t - w\|^2] = \Gamma_t \frac{\eta_1}{2\Gamma_1}[\|w_{0} - w\|^2 - \|w_T - w\|^2] \leq \Gamma_T \eta_1 \frac{1}{2}\|w_0 - w\|^2,
\]

since $\Gamma_1 = 1$. Substituting the above bounds into <ref>, we get
\begin{align}
F(w_T^{ag}) - F(w) &\leq \Gamma_T \eta_1  \frac{1}{2}\|w_0 - w\|^2 + \Gamma_T \sum_{t=1}^T \frac{\Upsilon_t(w)}{\Gamma_t}, ~\forall w \in \WW.
\end{align}

Now, setting $w = \ws$ and taking expectation yields
\begin{align}
\expec[F(w_T^{ag}) - F^*] &\leq \frac{\Gamma_T \eta_1  D^2}{2} + \Gamma_T\sum_{t=1}^T \frac{\expec \Upsilon_t(\ws)}{\Gamma_t} \\
&\leq \frac{\Gamma_T \eta_1  D^2}{2} + \Gamma_T\sum_{t=1}^T \left[\frac{1}{\Gamma_t}\left(\alpha_t \expec\langle b_t, \ws - w_{t-1} \rangle + \frac{2\alpha_t^2(\Sigma^2 + B^2)}{\eta_t - \beta \alpha_t^2}\right) \right] \\
&\leq \frac{\Gamma_T \eta_1  D^2}{2} + \Gamma_T\sum_{t=1}^T \left[\frac{1}{\Gamma_t}\left(\alpha_t BD + \frac{2\alpha_t^2(\Sigma^2 + B^2)}{\eta_t - \beta \alpha_t^2}\right) \right],
\end{align}

where we used conditional independence of $N_t$ and $\ws - w_{t-1}$ given $w_{t-1}$, Young's inequality, Cauchy-Schwartz, and the definitions of $B^2$ and $\Sigma^2$. This establishes the first claim of the proposition. The second and third claims are simple corollaries, which can be verified as in <cit.> and the ensuing discussion.

Let $\epsilon > 0$ and assume $F$ is convex and $\beta$-smooth. Then, there are algorithmic parameters such that <ref>

is $\frac{\epsilon^2}{2}$-zCDP. Further,

if $n \geq T := \left \lceil \min\left\{
% \left(\frac{\beta D}{r}\right)^{(k-1)/(5k-1)}
\left(\frac{\beta D}{r}\right)^{2k/(5k-1)}
\left(\frac{\epsilon n}{\sqrt{d}}\right)^{(2k-2)/(5k-1)}, \sqrt{\frac{\beta D}{r}}n^{1/4}\right\} \right \rceil$, then,

\[
% \eplac \lesssim r D\left[\frac{1}{\sqrt{n}} +
% \min\left\{
% \left(\frac{\sqrt{d}}{\epsilon n}\right)^{4(k-1)/(5k-1)} \left(\frac{\beta D}{r}\right)^{(k-1)/(5k-1)}
% , \left(\frac{\sqrt{d}}{\epsilon n} \right)^{(k-1)/k}
% \left(\frac{\beta D \sqrt{n}}{r}\right)^{(k-1)/4k}
% \right\}
% \right].
\eplac \lesssim
r D\left[
\frac{1}{\sqrt{n}} +
% \min\left\{
% \left(\left(\frac{\beta D}{r}\right)^{1/4} \frac{\sqrt{d}}{\epsilon n} \right)^{\frac{4(k-1)}{5k-1}}
% , \left(\left(\frac{\beta D \sqrt{n}}{r}\right)^{1/4} \frac{\sqrt{d}}{\epsilon n}  \right)^{\frac{k-1}{k}}
% \right\}
\max\left\{
\left(\left(\frac{\beta D}{r}\right)^{1/4} \frac{\sqrt{d}}{\epsilon n} \right)^{\frac{4(k-1)}{5k-1}}
, \left(\frac{\sqrt{d}}{\epsilon n}\right)^{\frac{k-1}{k}}
% \left(\left(\frac{\beta D \sqrt{n}}{r}\right)^{1/4} \frac{\sqrt{d}}{\epsilon n}  \right)^{\frac{k-1}{k}}
\right\}
\right].
\]

Privacy: Choose $\sigma^2 = \frac{4C^2 T^2}{\epsilon^2 n^2}$. First, the collection of all $\tilt(w_t^{md}), ~t \in [T]$ is $\frac{\epsilon^2}{2}$-zCDP: since the batches of data drawn in each iteration are disjoint, it suffices (by parallel composition [59]) to show that $\tilt(w_t^{md})$ is $\frac{\epsilon^2}{2}$-zCDP for all $t$. Now, the $\ell_2$ sensitivity of each clipped gradient update is bounded by $\Delta = \sup_{w, X \sim X'} \|\frac{T}{n} \sum_{x \in \mathcal{B}_t} \Pi_{C}(\nabla f(w, x)) - \sum_{x' \in \mathcal{B}'_t} \Pi_{C}(\nabla f(w, x'))\| =  \sup_{w, x, x'} \|\frac{T}{n} \Pi_{C}(\nabla f(w, x)) - \Pi_{C}(\nabla f(w, x'))\| \leq \frac{2CT}{n}$. Thus, $\tilt(w_t^{md})$ is $\frac{\epsilon^2}{2}$-zCDP by Proposition <ref>. Second, the iterates $w_t^{ag}$ are deterministic functions of $\tilt(w_t^{md})$, so the post-processing property of differential privacy [29, 16] ensures that <ref> is $\frac{\epsilon^2}{2}$-zCDP.

Excess risk:
Consider round $t \in [T]$ of <ref>, where <ref> is run on input data $\{\nabla f(w_t, x_i^t)\}_{i=1}^{n/T}$. Denote the bias of <ref> by $b_t:= \expec \tilt(w_t) - \nabla F(w_t)$, where $\tilt(w_t) = \widetilde{\nu}$ in the notation of <ref>. Also let $\hilt(w_t) := \hat{\mu}$ (in the notation of Lemma <ref>) and denote the noise by $N_t = \tilt(w_t) - \nabla F(w_t) - b_t = \tilt(w_t) - \expec \tilt(w_t)$. Then we have $B := \sup_{t \in [T]}\|b_t\| \leq \frac{r^k}{(k-1) C^{k-1}}$ and $\Sigma^2 := \sup_{t \in [T]} \expec[\|N_t\|^2] \leq d\sigma^2 + \frac{r^2 T}{n} \lesssim \frac{d C^2 T^2}{\epsilon^2 n^2} + \frac{r^2 T}{n}$, by Lemma <ref>. Plugging these estimates for $B$ and $\Sigma^2$ into Proposition <ref> and setting $C = r(\frac{\epsilon n}{\sqrt{dT}})^{1/k}$, we get
\begin{align}
\label{eq: thingg}
\expec F(w_T^{ag}) - F^* &\lesssim \frac{\beta D^2}{T^2} + \frac{D(\Sigma + B)}{\sqrt{T}} + BD \nonumber \\
&\lesssim \frac{\beta D^2}{T^2} + \frac{CD \sqrt{dT}}{\epsilon n} + \frac{r D}{\sqrt{n}} + \frac{r^k D}{C^{k-1}} \nonumber \\
&\lesssim \frac{\beta D^2}{T^2} + r D\left[\frac{1}{\sqrt{n}} + \left(\frac{\sqrt{dT}}{\epsilon n}\right)^{(k-1)/k}\right].
\end{align}

Now, our choice of $T$ implies that $\frac{\beta D^2}{T^2} \leq r D\left[\frac{1}{\sqrt{n}} + \left(\frac{\sqrt{dT}}{\epsilon n}\right)^{(k-1)/k}\right]$ and we get the result upon plugging in $T$.

§.§ Noisy Clipped SGD for Strongly Convex Losses (<ref>)

We begin by presenting <ref>, which is a generalized version of <ref> that allows for any  and will be useful in our analysis.

Noisy Clipped SGD Framework for Heavy-Tailed SCO

[1]
Data $X \in \XX^n$, $T \leq n$,

stepsizes $\{\eta_t\}_{t=0}^{T}$, averaging weights $\{\zeta_t\}_{t=0}^{T}$.

Initialize $w_0 \in \WW$.
$t \in \{0, 1, \cdots, T\}$
Draw new batch $\mathcal{B}_t$ (without replacement) of $n/T$ samples from $X$.
$\tilt(w_t) := \texttt{MeanOracle}(\{\nabla f(w_t, x)\}_{x \in \mathcal{B}_t}; \frac{n}{T}; \frac{\epsilon^2}{2})$
$w_{t+1} = \Pi_{\WW}\left[w_t - \eta_t \tilt(w_t) \right]
%\prox_{\eta_t f^1}\left(w_t - \eta_t \tilt^0(w_t) \right)

Output: $\widehat{w}_T := \frac{1}{Z_T} \sum_{t=0}^{T} \zeta_t w_{t+1}$, where $Z_T = \sum_{t=0}^T \zeta_t$.

In Proposition <ref>, we provide the convergence guarantees for <ref> in terms of the bias and variance of the .

Let $F: \WW \to \mathbb{R}$ be $\mu$-strongly convex and $\beta$-smooth with condition number $\kappa := \frac{\beta}{\mu}$. Let $w_{t+1} := \Pi_{\WW}[w_t - \eta_t \tilt(w_t)]$, where $\tilt(w_t) = \nabla F(w_t) + b_t + N_t$, such that the bias and noise (which can depend on $w_t$ and the samples drawn) satisfy $\|b_t\| \leq B$ (with probability $1$), $\expec N_t = 0$, $\expec\|N_t\|^2 \leq \Sigma^2$ for all $t \in [T-1]$, and that $\{N_t\}_{t=1}^T$ are independent. Then, there exist stepsizes $\{\eta_t\}_{t=1}^T$ and weights $\{\zeta_t\}_{t=0}^T$ such that the average iterate $\widehat{w}_T := \frac{1}{\sum_{t=0}^T \zeta_t} \sum_{t=0}^T \zeta_t w_{t+1}$ satisfies

\[
\EPL \leq 32\beta D^2 \exp\left(-\frac{T}{4 \kappa}\right) + \frac{72\Sigma^2}{\mu T} + \frac{2 B^2}{\mu}.
\]

Define $g(w_t) = -\frac{1}{\eta_t}(w_{t+1} - w_t)$. Then

\begin{align}
\label{eq:ving}
\expec\|w_{t+1} - \ws\|^2 &= \expec\|w_t - \eta_t g(w_t) - \ws\|^2 \nonumber \\
&= \expec\|w_t - \ws\|^2 - 2\eta_t \expec\langle g(w_t), w_t - \ws \rangle + \eta_t^2 \expec\|g(w_t)\|^2.
\end{align}

Now, conditional on all randomness, we use smoothness and strong convexity to write:

\begin{align*}
F(w_{t+1}) - F(\ws) &= F(w_{t+1}) - F(w_t) + F(w_t) - F(\ws) \\
&\leq \langle F(w_t), w_{t+1} - w_t \rangle + \frac{\beta}{2}\|w_{t+1} - w_t\|^2 + \langle \nabla F(w_t), w_t - \ws \rangle - \frac{\mu}{2}\|w_t - \ws\|^2 \\
&= \langle \tilt(w_t), w_{t+1} - \ws \rangle + \langle \nabla F(w_t) - \tilt(w_t), w_{t+1} - \ws \rangle + \frac{\beta \eta_t^2}{2}\|g(w_t)\|^2 - \frac{\mu}{2}\|w_t - \ws\|^2  \\
&\leq \langle g(w_t), w_{t+1} - \ws \rangle + \langle \nabla F(w_t) - \tilt(w_t), w_{t+1} - \ws \rangle + \frac{\beta \eta_t^2}{2}\|g(w_t)\|^2 - \frac{\mu}{2}\|w_t - \ws\|^2 \\
&= \langle g(w_t), w_{t+1} - w_t \rangle + \langle g(w_t), w_{t} - \ws \rangle - \langle b_t + N_t, w_{t+1} - \ws \rangle + \frac{\beta \eta_t^2}{2}\|g(w_t)\|^2 - \frac{\mu}{2}\|w_t - \ws\|^2 \\
&= \langle g(w_t), w_{t} - \ws \rangle - \langle b_t + N_t, w_{t+1} - \ws \rangle + \left(\frac{\beta \eta_t^2}{2} - \eta_t\right)\|g(w_t)\|^2 - \frac{\mu}{2}\|w_t - \ws\|^2,
\end{align*}

where we used the fact that $\langle \Pi_{\WW}(y) - x, \Pi_{\WW}(y) - y \rangle \leq 0$ for all $x \in \WW, y\in \mathbb{R}^d$ (c.f. <cit.>) to obtain the last inequality. Thus,
\[
-2 \eta_t \expec \langle g(w_t), w_{t} - \ws \rangle \leq -2\eta_t \expec[F(w_{t+1}) - F^*] + 2\eta_t \expec\left[- \langle b_t + N_t, w_{t+1} - \ws \rangle + \left(\frac{\beta \eta_t^2}{2} - \eta_t\right)\|g(w_t)\|^2 - \frac{\mu}{2}\|w_t - \ws\|^2 \right].
\]

Combining the above inequality with <ref>, we get
\begin{equation}
\label{eq:bing}
\expec\|w_{t+1} - \ws\|^2 \leq (1 - \mu \eta_t)\expec\|w_t - \ws\|^2 - 2\eta_t\expec[F(w_{t+1}) - F^*] - 2\eta_t \expec \langle b_t + N_t, w_{t+1} - \ws \rangle + 2\eta_t\left(\frac{\eta_t^2 \beta}{2} - \eta_t\right)\expec\|g(w_t)\|^2.
\end{equation}

Next, consider

\begin{align*}
|\expec \langle b_t + N_t, w_{t+1} - \ws \rangle| &\leq |\expec \langle b_t + N_t, w_{t+1} - w_t \rangle| + |\expec \langle b_t + N_t, w_t - \ws \rangle|\\
&= |\expec \langle b_t + N_t, w_{t+1} - w_t \rangle| + |\expec \langle b_t, w_t - \ws \rangle| \\
&\leq |\expec \langle b_t + N_t, w_{t+1} - w_t \rangle| + \frac{B^2}{\mu} + \frac{\mu}{4} \expec\|w_t - \ws\|^2
\end{align*}

by independence of $N_t$ (which has zero mean) and $w_t - \ws$, and Young's inequality. Next, note that $v := w_t -\eta_t(\nabla F(w_t) + b_t)$ is independent of $N_t$, so $\expec \langle N_t, \Pi_{\WW}(v)\rangle = 0$. Thus,

\begin{align*}
|\expec \langle N_t, w_{t+1} - w_t \rangle| &= |\expec \langle N_t, w_{t+1} \rangle | \\
&= |\expec \langle N_t, \Pi_{\WW}\left[w_t - \eta_t\left(\nabla F(w_t) + b_t + N_t\right)\right] \rangle | \\
&= |\expec \langle N_t, \Pi_{\WW}\left[v - \eta_t N_t\right] \rangle |\\
&= |\expec \langle N_t, \Pi_{\WW}\left[v] - \Pi_{\WW}[v - \eta_t N_t\right] \rangle | \\
&\leq \expec\left[\|N_t\| \|\Pi_{\WW}\left[v] - \Pi_{\WW}[v - \eta_t N_t\right]\|\right]\\
&\leq \expec\left[\|N_t\| \|\eta_t N_t\|\right] \\
&\leq \eta_t \Sigma^2,
\end{align*}

by Cauchy-Schwartz and non-expansiveness of projection.
\begin{align*}
|\expec \langle b_t, w_{t+1} - w_t \rangle| &= |\expec \langle b_t, -\eta_t g(w_t) \rangle|\\
&\leq \frac{B^2}{\mu} + \frac{\eta_t^2 \mu}{4}\expec\|g(w_t)\|^2,
\end{align*}

by Young's inequality. Therefore,
\begin{align*}
-2\eta_t \expec \langle b_t + N_t, w_{t+1} - \ws \rangle \leq 2\eta_t\left[\frac{2B^2}{\mu} + \frac{\eta_t^2 \mu}{4}\expec\|g(w_t)\|^2 + \eta_t \Sigma^2 + \frac{\mu}{4}\expec\|w_t - \ws\|^2\right].
\end{align*}

Plugging this bound back into <ref> and choosing $\eta_t \leq \frac{1}{\beta} \leq \frac{1}{\mu}$
\begin{align*}
\expec\|w_{t+1} - \ws\|^2 &\leq \left(1 - \frac{\mu \eta_t}{2}\right)\expec\|w_t - \ws\|^2 - 2\eta_t\expec[F(w_{t+1}) - F^*] + \frac{4\eta_t B^2}{\mu} + 2\eta_t^2 \Sigma^2
+ 2\eta_t\left(\frac{\eta_t^2 \beta}{2} - \eta_t + \frac{\eta_t^2 \mu}{4}\right)\expec\|g(w_t)\|^2 \\
&\leq \left(1 - \frac{\mu \eta_t}{2}\right)\expec\|w_t - \ws\|^2 - 2\eta_t\expec[F(w_{t+1}) - F^*] + \frac{4\eta_t B^2}{\mu} + 2\eta_t^2 \Sigma^2.
\end{align*}

Next, we apply Lemma <ref> (see below) with $r_t := \expec \|w_{t} - \ws\|^2, ~s_t := \expec F(w_{t+1}) - F^* - \frac{2B^2}{\mu}$, $a := \frac{\mu}{2}$, $b := 2$, $c = 2 \Sigma^2$, and $g = \beta$. We may assume $s_t \geq 0$ for all $t$: if this inequality breaks for some $t$, then simply return $w_{t+1}$ instead of $\widehat{w}_T$ to obtain $\expec F(w_t) - F^* < \frac{2B^2}{\mu}$. Thus,
\[
\frac{1}{\Gamma_T}\sum_{t=0}^T \gamma_t \expec[F(w_{t+1}) - F^*] \leq \frac{1}{2}\left[32 \beta D^2 \exp\left(\frac{-\mu T}{4 \beta}\right) + \frac{144 \Sigma^2}{\mu T} + \frac{2 B^2}{\mu}\right]
\]

Finally, Jensen's inequality yields the proposition.

Let $b > 0$, let $a, c \geq 0,$ and $\{\eta_t\}_{t \geq 0}$ be non-negative step-sizes such that $\eta_t \leq \frac{1}{g}$ for all $t \geq 0$ for some parameter $g \geq a$. Let $\{r_t\}_{t \geq 0}$ and $\{s_t\}_{t \geq 0}$ be two non-negative sequences of real numbers which satisfy
\[
r_{t+1} \leq (1 - a \eta_t)r_t - b \eta_t s_t + c \eta_t^2
\]

for all $t \geq 0.$
Then there exist particular choices of step-sizes $\eta_t \leq \frac{1}{g}$ and averaging weights $\zeta_t \geq 0$ such that
\[
\frac{b}{\Gamma_T}\sum_{t=0}^T s_t \zeta_t + a r_{T+1} \leq 32 g r_0 \exp\left(\frac{-aT}{2g}\right) + \frac{36c}{aT},
\]

where $\Gamma_T := \sum_{t=0}^{T} \gamma_t.$

We are now prepared to prove <ref>.

Let $\epsilon > 0$, and assume $F$ is $\mu$-strongly convex and $\beta$-smooth with $\kappa = \frac{\beta}{\mu} \leq n/\ln(n)$.

Then, there are parameters such that <ref> instantiated with  <ref> is $\frac{\epsilon^2}{2}$-zCDP, and
\begin{equation}
%\label{eq: smooth sc upper l2 clip}
%\EPL = \mathcal{O}\left(\frac{r^2}{\mu}\left(\frac{1}{n} + \left(\frac{\sqrt{d \kappa \ln(n)}}{\epsilon n}\right)^{\frac{2k-2}{k}} \right)\right).
\EPL \leq \frac{Ar^2}{\mu}\left(\frac{1}{n} + \left(\frac{\sqrt{d \kappa \ln(n)}}{\epsilon n}\right)^{\frac{2k-2}{k}} \right).
\end{equation}

Privacy: Choose $\sigma^2 = \frac{4C^2 T^2}{\epsilon^2 n^2}$. Since the batches of data drawn in each iteration are disjoint, it suffices (by parallel composition [59]) to show that $\tilt(w_t)$ is $\frac{\epsilon^2}{2}$-zCDP for all $t$. Now, the $\ell_2$ sensitivity of each clipped gradient update is bounded by $\Delta = \sup_{w, X \sim X'} \|\frac{T}{n} \sum_{x \in \mathcal{B}_t} \Pi_{C}(\nabla f(w, x)) - \sum_{x' \in \mathcal{B}'_t} \Pi_{C}(\nabla f(w, x'))\| =  \sup_{w, x, x'} \|\frac{T}{n} \Pi_{C}(\nabla f(w, x)) - \Pi_{C}(\nabla f(w, x'))\| \leq \frac{2CT}{n}$. Hence Proposition <ref> implies that the algorithm is $\frac{\epsilon^2}{2}$-zCDP.

Excess risk: For any iteration $t \in [T]$, denote the bias of <ref> by $b_t:= \expec \tilt(w_t) - \nabla F(w_t)$, where $\tilt(w_t) = \widetilde{\nu}$ in the notation of <ref>. Also let $\hilt(w_t) := \hat{\nu}$ (in the notation of <ref>) and denote the noise by $N_t = \tilt(w_t) - \nabla F(w_t) - b_t = \tilt(w_t) - \expec \tilt(w_t)$. Then we have $B := \sup_{t \in [T]}\|b_t\| \leq \frac{r^k}{(k-1) C^{k-1}}$ and $\Sigma^2 := \sup_{t \in [T]} \expec[\|N_t\|^2] \leq d\sigma^2 + \frac{r^2 T}{n} \lesssim \frac{d C^2 T^2}{\epsilon^2 n^2} + \frac{r^2 T}{n}$, by <ref>. Plugging these bias and variance estimates into <ref>, we get
\[
\EPL \lesssim \beta D^2 \exp\left(- \frac{T}{4 \kappa}\right) + \frac{1}{\mu T}\left(\frac{d C^2 T^2}{\epsilon^2 n^2} + \frac{r^2 T}{n}\right) + \frac{r^{2k}}{C^{2k-2} \mu}.
\]

Choosing $C = r\left(\frac{\epsilon^2 n^2}{dT}\right)^{1/2k}$ implies
\[
\EPL \lesssim \beta D^2 \exp\left(- \frac{T}{4 \kappa}\right) + \frac{r^2}{\mu}\left(\frac{1}{n} + \left(\frac{dT}{\epsilon^2 n^2}\right)^{(k-1)/k}\right).
\]

Finally, choosing $T = \left\lceil 4\kappa\ln\left(\frac{\mu \beta D^2}{r^2}\left(n + \left(\frac{\epsilon^2 n^2}{d}\right)^{(k-1)/k}\right)\right) \right\rceil \lesssim \kappa \ln(n)$ yields the result.

§ DETAILS AND PROOFS FOR <REF>

Consider <ref>, a generalization of <ref> which allows for arbitrary .

Noisy Proximal SGD Framework for Heavy-Tailed SO

[1]
Data $X \in \XX^n$, $T \leq n$,

$\texttt{MeanOracle}$ (and truncation/minibatch parameters),

privacy parameter $\rho = \epsilon^2/2$, stepsizes $\{\eta_t\}_{t=0}^{T-1}$.

Initialize $w_0 \in \WW$.
$t \in \{0, 1, \cdots, T-1\}$
Draw new batch $\mathcal{B}_t$ (without replacement) of $n/T$ samples from $X$.
$\tilt^0(w_t) := \texttt{MeanOracle}(\{\nabla f^0(w_t, x)\}_{x \in \mathcal{B}_t}; \frac{n}{T}; \frac{\epsilon^2}{2})$
$w_{t+1} =
%\Pi_{\WW}\left[w_t - \eta \tilt(w_t) \right]
\prox_{\eta_t f^1}\left(w_t - \eta_t \tilt^0(w_t) \right)

Output: $w_T$.

Proposition <ref> provides a convergence guarantee for <ref> in terms of the bias and variance of the .

Consider <ref> with biased, noisy stochastic gradients:

$\tilt^0(w_t) = \nabla F^0(w_t) + b_t + N_t$, and stepsize $\eta = \frac{1}{2\beta}$. Assume that the bias and noise (which can depend on $w_t$ and the samples drawn) satisfy $\|b_t\| \leq B$ (with probability $1$), $\expec N_t = 0$, $\expec\|N_t\|^2 \leq \Sigma^2$ for all $t \in [T-1]$, and that $\{N_t\}_{t=1}^T$ are independent. Assume further that $F$ is $\mu$-Proximal-PL, $F^0$ is $\beta$-smooth, and $F(w_0) - F^* \leq \Delta$. Then,
\[
\expec F(w_T) - F^* \leq \left(1 - \frac{\mu}{2\beta}\right)^T \Delta + \frac{4(B^2 + \Sigma^2)}{\mu}.
\]

Our proof extends the ideas in [52] to generic biased and noisy gradients without using Lipschitz continuity of $f$.

By $\beta$-smoothness, for any $r \in [T-1]$, we have
\begin{align}
\expec F(w_{r+1}) &= \expec[F^0(w_{r+1}) + f^1(w_r) + f^1(w_{r+1}) - f^1(w_r)] \nonumber \\
&\leq \expec\left\{F(w_r) + \left[\langle \tilr(w_r), w_{r+1} - w_r \rangle + \frac{\beta}{2}\|w_{r+1} - w_r\|^2 + f^1(w_{r+1}) - f^1(w_r)\right] \right\}\nonumber \\
&\;\;\; + \expec \langle \nabla F^0(w_r) - \tilr(w_r), w_{r+1} - w_r \rangle \nonumber \\
&= \expec F(w_r) + \expec \left[\langle \nabla F^0(w_r), w_{r+1} - w_r \rangle + \frac{\beta}{2}\|w_{r+1} - w_r\|^2 + f^1(w_{r+1}) - f^1(w_r) + \langle b_r + N_r, w_{r+1} - w_r \rangle \right] \nonumber \\
&\;\;\; - \expec \langle b_r + N_r, w_{r+1} - w_r \rangle \nonumber \\
&\leq \expec F(w_r) + \expec \left[\langle \nabla F^0(w_r), w_{r+1} - w_r \rangle + \beta\|w_{r+1} - w_r\|^2 + f^1(w_{r+1}) - f^1(w_r) + \langle b_r + N_r, w_{r+1} - w_r \rangle \right] \nonumber \\
&\;\;\; + \frac{B^2 + \Sigma^2}{\beta},
\label{eq27}
\end{align}

where we used Young's inequality to bound
\begin{align}
-\expec \langle b_r + N_r, w_{r+1} - w_r \rangle &\leq \frac{B^2 + \Sigma^2}{\beta} + \frac{\beta}{2}\|w_{r+1} - w_r\|^2.
\end{align}

Next, we will bound $\expec \left[\langle \nabla F^0(w_r), w_{r+1} - w_r \rangle + \beta\|w_{r+1} - w_r\|^2 + f^1(w_{r+1}) - f^1(w_r) + \langle b_r + N_r, w_{r+1} - w_r \rangle \right]$.

$H^{\text{priv}}_r(y):= \langle \nabla F^0(w_r), y - w_r \rangle + \beta\|y - w_r\|^2 + f^1(y) - f^1(w_r) + \langle b_r + N_r, y - w_r \rangle$ and $H_r(y):=  \langle \nabla F^0(w_r), y - w_r \rangle + \beta\|y - w_r\|^2 + f^1(y) - f^1(w_r)$.
Note that $H_r$ and $H^{\text{priv}}_r$ are $2\beta$-strongly convex. Denote the minimizers of these two functions by $y_*$ and $y_*^{\text{priv}}$ respectively. Now, conditional on $w_r$ and $N_r + b_r$, we claim that
\begin{equation}
\label{claimyyy}
H_r(y_*^{\text{priv}}) - H_r(y_*) \leq \frac{\|N_r + b_r\|^2}{2\beta}.
\end{equation}

To prove <ref>, we will need the following lemma:

Let $H(y), h(y)$ be convex functions on some convex closed set $\mathcal{Y} \subseteq \mathbb{R}^d$ and suppose that $H$ is $2\beta$-strongly convex. Assume further that $h$ is $L_{h}$-Lipschitz. Define $y_{1} = \arg\min_{y \in \mathcal{Y}} H(y)$ and $y_{2} = \arg\min_{y \in \mathcal{Y}} [H(y) + h(y)]$. Then $\|y_{1} - y_{2}\|_2 \leq \frac{L_{h}}{2\beta}.$

We apply <ref> with $H(y):= H_r(y), ~h(y):= \langle N_r + b_r, y\rangle$, $L_h = \|N_r + b_r\|$, $y_1 = y_*$, and $y_2 = y_*^{\text{priv}}$ to get
\[
\|y_* - y_*^{\text{priv}}\| \leq \frac{\|N_r + b_r\|}{2\beta}.
\]
On the other hand,
\[
H_r^{\text{priv}}(y_*^{\text{priv}}) = H_r(y_*^{\text{priv}}) + \langle N_r + b_r, y_*^{\text{priv}} \rangle \leq H_r^{\text{priv}}(y_*) = H_r(y_*) + \langle N_r + b_r, y_* \rangle.
\]

Combining these two inequalities yields
\begin{align}
H_r(y_*^{\text{priv}}) - H_r(y_*) &\leq \langle N_r + b_r,   y_* - y_*^{\text{priv}}  \rangle \nonumber \\
&\leq \|N_r + b_r\| \|y_* - y_*^{\text{priv}}\| \nonumber \\
&\leq \frac{\|N_r + b_r\|^2}{2\beta},
\end{align}

as claimed. Also, note that $w_{r+1} = y_*^{\text{priv}}$. Hence
\begin{equation}
\expec \left[\langle \nabla F^0(w_r), w_{r+1} - w_r \rangle + \beta\|w_{r+1} - w_r\|^2 + f^1(w_{r+1}) - f^1(w_r) + \langle b_r + N_r, w_{r+1} - w_r \rangle \right] =\expec\left[ \min_{y \in \mathbb{R}^d} H_r^{\text{priv}}(y) \right]
\end{equation}

\begin{align}
\expec\left[ \min_{y \in \mathbb{R}^d} H_r^{\text{priv}}(y) \right]&\leq \expec \left[\min_y \left\{ \langle \nabla F^0(w_r), y - w_r \rangle + \beta \|y - w_r\|^2 + f^1(y) - f^1(w_r)\right\} \right] + \frac{\Sigma^2 + B^2}{\beta} \\
&\leq -\frac{\mu}{2 \beta}\expec\left[F(w_r) - F^* \right] + \frac{\Sigma^2 + B^2}{\beta},
\end{align}

where we used the assumptions that $F$ is $\mu$-PPL and $F^0$ is $2\beta$-smooth in the last inequality.

Plugging the above bounds back into <ref>, we obtain
\begin{equation}
\expec F(w_{r+1}) \leq \expec F(w_r) - \frac{\mu}{2\beta}[F(w_r) - F^*] + \frac{2(\Sigma^2 + B^2)}{\beta},
\end{equation}

\begin{equation}
\label{eq: thingz}
\expec[F(w_{r+1}) - F^*] \leq \expec[F(w_r) - F^*](1 - \frac{\mu}{2\beta}) + \frac{2(\Sigma^2 + B^2)}{\beta}.
\end{equation}

Using <ref> recursively and summing the geometric series, we get
\begin{equation}
\expec[F(w_T) - F^*] \leq \Delta \left(1 - \frac{\mu}{2\beta}\right)^T + \frac{4(\Sigma^2 + B^2)}{\mu}.
\end{equation}

$F = F^0 + f^1$ be $\mu$-PPL for $\beta$-smooth $F^0$,
with $
%\kappa =
\frac{\beta}{\mu} \leq n/\ln(n)$.

Then, there are parameters such that <ref> is $\frac{\epsilon^2}{2}$-zCDP, and (for an absolute constant $A$):

\[
\small
% \EPLL = \mathcal{O}\left(\frac{r^2}{\mu}\left(\left(\frac{\sqrt{d}}{\epsilon n} \kappa \ln(n) \right)^{\frac{2k-2}{k}} + \frac{\kappa \ln(n)}{n}\right) \right).
\EPLL \leq \frac{Ar^2}{\mu}\left(\left(\frac{\sqrt{d}}{\epsilon n} \left(\frac{\beta}{\mu}\right) \ln(n) \right)^{\frac{2k-2}{k}} + \frac{(\beta/\mu) \ln(n)}{n}\right).
\]

Let us denote $\kappa = \beta/\mu$. We choose $\sigma^2 = \frac{4C^2 T^2}{\epsilon^2 n^2}$.

Privacy: By parallel composition (since each sample is used only once) and the post-processing property of DP (since the iterates are deterministic functions of the output of ), it suffices to show that $\tilt(w_t)$ is $\frac{\epsilon^2}{2}$-zCDP for all $t \geq 0$. By our choice of $\sigma^2$ and Proposition <ref>, $\tilt(w_t)$ is $\frac{\epsilon^2}{2}$-zCDP, since it's sensitivity is bounded by $\sup_{X \sim X', w} \frac{T}{n}\left\|\sum_{x \in \mathcal{B}_t} \Pi_C[\nabla f^0(w, x)] -  \sum_{x' \in \mathcal{B'}_t} \Pi_C[\nabla f^0(w, x')]\right\| \leq \frac{T}{n} \sup_{x, x', w} \| \Pi_C[\nabla f^0(w, x)] -  \Pi_C[\nabla f^0(w, x')]\| \leq \frac{2C T}{n}$.

Excess risk: For any iteration $t \in [T]$, denote the bias of <ref> by $b_t:= \expec \tilt(w_t) - \nabla F(w_t)$, where $\tilt(w_t) = \widetilde{\nu}$ in the notation of <ref>. Also let $\hilt(w_t) := \hat{\nu}$ (in the notation of <ref>) and denote the noise by $N_t = \tilt(w_t) - \nabla F(w_t) - b_t = \tilt(w_t) - \expec \tilt(w_t)$. Then we have $B := \sup_{t \in [T]}\|b_t\| \leq \frac{r^k}{(k-1) C^{k-1}}$ and $\Sigma^2 := \sup_{t \in [T]} \expec[\|N_t\|^2] \leq d\sigma^2 + \frac{r^2 T}{n} \leq \frac{4d C^2 T^2}{\epsilon^2 n^2} + \frac{r^2 T}{n}$, by <ref>. Plugging these bounds on $B^2$ and $\Sigma^2$ into <ref>, and choosing $T = 2\left\lceil \kappa \ln\left(\frac{\Delta \mu}{B^2 + \Sigma^2} \right)\right\rceil \lesssim \kappa \ln(n)$ where $\Delta \geq F(w_0) - F^*$,
we have:
\[
\EPLL \leq \frac{5(B^2 + \Sigma^2)}{\mu} \leq \frac{5}{\mu}\left(\frac{2r^2 T}{n} + \frac{2 r^{2k}}{(k-1)^2 C^{2k-2}} + \frac{2d C^2 T^2}{\epsilon^2 n^2}\right),
\]

for any $C > 0$.
Choosing $C = r\left(\frac{\epsilon^2 n^2}{d T^2}\right)^{1/2k}$ makes the last two terms in the above display equal, and we get
\[
\EPLL =
%\mathcal{O}\left(\frac{r^2}{\mu}\left(\left(\frac{\sqrt{d}}{\epsilon n} \kappa \ln(n) \right)^{(2k-2)/k} + \frac{\kappa \ln(n)}{n}\right) \right),
\EPLL \leq \frac{Ar^2}{\mu}\left(\left(\frac{\sqrt{d}}{\epsilon n} \kappa \ln(n) \right)^{\frac{2k-2}{k}} + \frac{\kappa \ln(n)}{n}\right)
\]

for an absolute constant $A$, as claimed.

§ SHUFFLE PRIVATE MEAN ESTIMATORS

The definitions of DP and zCDP given earlier do not dictate how the algorithm $\Al$ operates. In particular, they allow $\Al$ to send sensitive data to a third party curator/analyst, who can then add noise to the data. However, in certain practical applications (e.g. federated learning [46]), there is no third party that individuals trust to handle their sensitive data. On the other hand, it is often more realistic to have a secure shuffler (a.k.a. mixnet):

in each iteration of the algorithm, the shuffler receives encrypted, noisy reports (e.g. encrypted, noisy, private stochastic gradients) from each user and applies a uniformly random permutation to the $n$ reports, thereby anonymizing them (and amplifying privacy) [12, 23, 32, 37]. An algorithm is shuffle private if all of these “shuffled” reports are DP:
(Shuffle Differential Privacy [12, 23])

A randomized algorithm is $(\epsilon, \delta)$-shuffle DP (SDP) if the collection of reports output by the shuffler satisfies Definition <ref>.

In this next two subsections, we present two SDP algorithms for DP heavy-tailed mean estimation. The first is an SDP version of <ref> and the second is an SDP version of the coordinate-wise protocol of [50, 45]. Both of our algorithms offer the same utility guarantees as their zCDP counterparts (up to logarithms). In particular, this implies that the upper bounds obtained in the main body of this paper can also be attained via SDP protocols that do not require individuals to trust any third party with their sensitive data (assuming the existence of a secure shuffler).

§.§ $\ell_2$ Clip Shuffle Private Mean Estimator

For heavy-tailed SO problems satisfying <ref>, we propose using the SDP mean estimation protocol described in <ref>.  <ref> relies on the shuffle private vector summation protocol of [20], which is given in <ref>. The useful properties of <ref> are contained in <ref>.
$\ell_2$ Clip $\texttt{ShuffleMeanOracle1}(\{x_i\}_{i=1}^s; C; (\epsilon, \delta))$

[1]
$X = \{x_i\}_{i=1}^s$, $x_i = (x_{i,1}, \cdots, x_{i,d}) \in \mathbb{R}^d$, $C>0$, $(\epsilon, \delta) \in (\mathbb{R}_+ \times (0,1/2))$.
$i \in [s]$
$z_i := \Pi_{C}(x_i)$.

$\widetilde{\nu} := \mathcal{P}_{\text{vec}}(\{z_i\}_{i=1}^s; s; C; (\epsilon, \delta))$.
Output: $\widetilde{\nu}$.

$\mathcal{P}_{\text{vec}}$, a shuffle private protocol for vector summation

[1]
Input: database of $d$-dimensional vectors $\mathbf{X} = (\mathbf{x}_1, \cdots, \mathbf{x}_s$ with maximum norm bounded by $C > 0$; privacy parameters $(\epsilon, \delta)$.

procedure: Local Randomizer $\mathcal{R}_{\text{vec}}(\mathbf{x}_i)$
\beginALC@g
$j \in [d]$
Shift component to enforce non-negativity: $\mathbf{w}_{i,j} \gets \mathbf{x}_{i,j} + C$
$\mathbf{m}_j \gets \mathcal{R}_{1D}(\mathbf{w}_{i,j})$
Output labeled messages $\{(j, \mathbf{m}_j)\}_{j \in [d]}$
\endALC@g
end procedure
procedure: Analyzer $\mathcal{A}_{\text{vec}}(\mathbf{y})$
\beginALC@g
$j \in [d]$
Run analyzer on coordinate $j$'s messages $z_j \gets \mathcal{A}_{\text{1D}}(\mathbf{y}_j)$
Re-center: $o_j \gets z_j - L$
Output the vector of estimates $\mathbf{o} = (o_1, \cdots o_d)$
\endALC@g
end procedure

Let $\epsilon \leq 15, \delta \in (0, 1/2)$, $d, s \in \mathbb{N}$ and $C > 0$. There are choices of parameters $b, g, p$ for $\mathcal{P}_{1D}$ such that for an input data set $\mathbf{X} = (\bx_1, \cdots, \bx_s)$ of vectors with maximum norm $\|\bx_i\| \leq C$, the following holds:

1)  <ref> is $(\epsilon, \delta)$-SDP.

2) $\mathcal{P}_{\text{vec}}(\mathbf{X})$ is an unbiased estimate of $\sum_{i=1}^s \bx_i$ with bounded variance
\[
\expec\left[\left\|\mathcal{P}_{\text{vec}}(\mathbf{X}) - \sum_{i=1}^s \bx_i\right\|^2 \right] = \mathcal{O}\left(\frac{d C^2}{\epsilon^2} \ln^2(d/\delta) \right).
\]

By the post-processing property of DP, we immediately obtain:

Let $\{x_i\}_{i=1}^s \sim \DD^s$ have mean $\expec x_i = \nu$ and $\expec\|x_i\|^k \leq r^k$ for some $k \geq 2$. Denote the noiseless average of clipped samples in <ref> by $\widehat{\nu} := \frac{1}{n}\sum_{i=1}^s \Pi_C(x_i)$. Then, there exist algorithmic parameters such that <ref> is $(\epsilon, \delta)$-SDP and such that the following bias and variance bounds hold:

\begin{equation}
\label{eq: sdp bd14 bias}
\|\expec \widetilde{\nu} - \nu \| = \|\expec \widehat{\nu} - \nu \| \leq \expec \|\widehat{\nu} - \nu \| \leq \frac{r^k}{(k-1)C^{k-1}},
\end{equation}

\begin{equation}
\label{eq: sdp bd14 variance}
\expec\|\widetilde{\nu} - \expec \widetilde{\nu}\|^2 = \expec \| \widetilde{\nu} - \expec \widehat{\nu}\|^2 = \mathcal{O}\left(\frac{d C^2 \ln^2(d/\delta)}{\epsilon^2 s^2} + \frac{r^2}{s}\right).
\end{equation}

Privacy: The privacy claim is immediate from <ref> and the post-processing property of DP <cit.>.

Bias: The bias bound follows as in Lemma <ref>, since $\mathcal{P}_{\text{vec}}$ is an unbiased estimator (by <ref>).

Variance: We have

\begin{align*}
\expec\|\widetilde{\nu} - \expec\widetilde{\nu}\|^2 &= \expec \| \widetilde{\nu} - \expec \widehat{\nu}\|^2 \\
&= \expec \| \widetilde{\nu} - \widehat{\nu}\|^2 + \expec \| \widehat{\nu} - \expec \widehat{\nu}\|^2\\
&\lesssim \frac{d C^2 \ln^2(d/\delta)}{\epsilon^2 s^2} + \frac{1}{s}\expec\|\Pi_C(x_1) - \expec\Pi_C(x_1)\|^2 \\
&\leq \frac{d C^2 \ln^2(d/\delta)}{\epsilon^2 s^2} +  \frac{1}{s}\expec\|x_1 - \expec x_1\|^2\\
&\leq \frac{d C^2 \ln^2(d/\delta)}{\epsilon^2 s^2} + \frac{r^2}{s},
\end{align*}

where we used that the samples $\{x_i\}_{i=1}^s$ are i.i.d., the variance bound in <ref>, and <cit.>, which states that $\expec\|\Pi_C(X) - \expec \Pi_C(X)\|^2 \leq \expec \|X - \expec X\|^2$ for any random vector $X$.

Comparing <ref> to Lemma <ref>, we see that the bias and variance of the two s are the same up to logarithmic factors. Therefore, replacing <ref> by <ref> in our stochastic optimization algorithms yields SDP algorithms with excess risk that matches the bounds provided in this paper (via <ref>) up to logarithmic factors.

§.§ Coordinate-wise Shuffle Private Mean Estimation Oracle

For SO problems satisfying <ref>, we propose <ref> as a shuffle private mean estimation oracle. <ref> is a shuffle private variation of <ref>, which was employed by [50, 45].
Coordinate-wise Private $\texttt{MeanOracle2}(\{x_i\}_{i=1}^s; s; \tau; \frac{\epsilon^2}{2}; m)$ [50, 45]

[1]
$X = \{x_i\}_{i=1}^s$, $x_i = (x_{i,1}, \cdots, x_{i,d}) \in \mathbb{R}^d$, $\epsilon > 0, \tau > 0$, $m \in [s]$ such that $m$ divides $s$.
$j \in [d]$
Partition $j$-th coordinates of data into $m$ disjoint groups of size $s/m$.
$i \in [m]$
Clip data in $i$-th group: $Z_j^i := \left\{\Pi_{[-\tau, \tau]}(x_{(i-1)\frac{s}{m}+1, j}), \cdots, \Pi_{[-\tau, \tau]}(x_{i\frac{s}{m}, j})\right\}$.
Compute average of $Z_j^i$: $\widehat{\nu}_j^i := \frac{m}{s} \sum_{z \in Z_j^i} z$.
Compute median of group means: $\widehat{\nu}_j := \text{median}(\widehat{\nu}_j^1, \cdots, \widehat{\nu}_j^m)$.
Draw $u \sim \mathcal{N}(0, \sigma^2 \mathbf{I}_d)$, with $\sigma^2 = \frac{4\tau^2 m^2 d}{s^2 \epsilon^2}$.

Output: $\widetilde{\nu} = (\widehat{\nu}_1, \cdots, \widehat{\nu}_d) + u$.

The bias/variance and privacy properties of <ref> are summarized in <ref>.

Let $\{x_i\}_{i=1}^s \sim \DD^s$ have mean $\expec x_i = \nu$, $\|\nu\| \leq L$, and $\expec|x_{i,j} - \nu_j|^k \leq \gamma$ for some $k \geq 2$. Denote by $\widehat{\nu}$ the output of the non-private algorithm that is identical to <ref> except without the added Gaussian noise. Then, <ref> with $\sigma^2 = \frac{72 \tau^2 m^2 d}{\epsilon^2 s^2}$
is $\frac{\epsilon^2}{2}$-zCDP. Further, if $\tau \geq 2L$, then there is $m = \widetilde{\mathcal{O}}(1) \in [s]$ such that:

\begin{equation}
\|\expec \widetilde{\nu} - \nu \| = \|\expec \widehat{\nu} - \mu \| \leq \expec \|\widehat{\nu} - \nu \| = \widetilde{\mathcal{O}}\left(\sqrt{d}\left(\frac{\gamma^{1/k}}{\sqrt{s}} + \left(\frac{2}{\tau}\right)^{k-1} \gamma\right) \right) =: B,
\end{equation}

\begin{equation}
\expec\|N\|^2 = \expec\|\widetilde{\nu} - \expec \widetilde{\nu}\|^2 = \expec \| \widetilde{\nu} - \expec \widehat{\nu}\|^2 = \widetilde{\mathcal{O}}\left(B^2 + %\frac{\tau^2 d^2}{\epsilon^2 s^2}
\right).
\end{equation}

Coordinate-wise Shuffle Private $\texttt{ShuffleMeanOracle2}(\{x_i\}_{i=1}^s; s; \tau; \epsilon, \delta; m)$

[1]
$X = \{x_i\}_{i=1}^s$, $x_i = (x_{i,1}, \cdots, x_{i,d}) \in \mathbb{R}^d$, $\epsilon > 0, \delta \in (0, 1)$, $m \in [s]$.
$j \in [d]$
$\epsilon_j := \frac{\epsilon}{4\sqrt{2d \ln(1/\delta_j)}}, \delta_j := \frac{\delta}{2d}$.
Partition $j$-th coordinates of data into $m$ disjoint groups of size $s/m$.
$i \in [m]$
Clip $j$-th coordinate of data in $i$-th group: $Z_j^i := \left\{\Pi_{[-\tau, \tau]}(x_{(i-1)\frac{s}{m}+1, j}), \cdots, \Pi_{[-\tau, \tau]}(x_{i\frac{s}{m}, j})\right\}$.
Shift to enforce non-negativity: $Z_j^i \gets Z_j^i + (\tau, \cdots, \tau)$.
Compute noisy average of $s/m$ scalars in $Z_j^i$: $\widetilde{\nu}_j^i := \frac{m}{s} \mathcal{P}_{1D}(Z_j^i; \frac{s}{m}; 2\tau; \epsilon_j, \delta_j)$.
Re-center: $\widetilde{\nu}_j^i \gets \widetilde{\nu}_j^i - \tau$.
Compute median of noisy means: $\widetilde{\nu}_j := \text{median}(\widetilde{\nu}_j^1, \cdots, \widetilde{\nu}_j^m)$.

Output: $\widetilde{\nu} = (\widetilde{\nu}_1, \cdots, \widetilde{\nu}_d)$.

The $\mathcal{P}_{1D}$ subroutine used in <ref> is an SDP protocol for summing scalars that we borrow from [20]. It is outlined in <ref>. <ref> decomposes into a local randomizer $\mathcal{R}$ that individuals execute on their own data, and an analyzer component $\mathcal{A}$ that the shuffler executes. $\mathcal{S}(\mathbf{y})$ denotes the shuffled vector $\mathbf{y}$: i.e. the vector whose components are random permutations of the components of $\mathbf{y}$. We describe the privacy guarantee, bias, and variance of <ref> in <ref>.

$\mathcal{P}_{\text{1D}}$, a shuffle private protocol for summing scalars [20]

[1]
Scalars $Z = (z_1, \cdots z_s) \in [0,\tau]^s$; design parameters $g, b \in \mathbb{N}; p \in (0, \frac{1}{2})$.
procedure: Local Randomizer $\mathcal{R}_{1D}(z_i)$
\beginALC@g
$i \in [s]$
$\widebar{z}_i \gets \lfloor z_i g/\tau \rfloor$.
Sample rounding value $\eta_1 \sim \textbf{Ber}(z_i g/\tau - \widebar{z}_i)$.
Set $\hat{z}_i \gets \widebar{z}_i + \eta_1$.
Sample privacy noise value $\eta_2 \sim \textbf{Bin}(b,p)$.
Report $y_i = (y_{i,1}, \cdots, y_{i, g+b}) \in \{0,1\}^{g + b}$ containing $\hat{z}_i + \eta_2$ copies of $1$ and $g + b - (\hat{z}_i + \eta_2)$ copies of $0$.

\endALC@g
end procedure

procedure: Shuffler $\mathcal{S}(\mathbf{y})$
\beginALC@g
Shuffler receives $\mathbf{y}:= (y_1, \cdots, y_s)$, draws a uniformly random permutation $\pi$ of $[g+b] \times [s]$, and sends $\mathcal{S}(\mathbf{y}):=  (y_{\pi(1, 1)}, \cdots, y_{\pi(s, g+b)})$ to analyzer.
\endALC@g
end procedure
procedure: Analyzer $\mathcal{A}_{\text{1D}}(\mathcal{S}(\mathbf{y}))$
\beginALC@g
Output: $\frac{\tau}{g}[(\sum_{i=1}^{s}\sum_{l=1}^{b+g} y_{\pi(i,l)}) - pbs]$.
\endALC@g

Let $\{x_i\}_{i=1}^s \sim \DD^s$ have mean $\expec x_i = \nu$, $\|\nu\| \leq L$, and $\expec\|x_{i,j} - \nu_j\|^k \leq \gamma$ for some $k \geq 2$. Let $\epsilon \leq 8\ln(2d/\delta), ~\delta \in (0, 1/2)$, and choose $\tau \geq 2L$. Then, there exist choices of parameters ($g, b, p, m$) such that <ref> is $(\epsilon, \delta)$-SDP, and has bias and variance bounded as follows:

\begin{equation}
\|\expec \widetilde{\nu} - \nu \|
% = \|\expec \widehat{\nu} - \nu \| \leq \expec \|\widehat{\nu} - \nu \|
= \widetilde{\mathcal{O}}\left(\sqrt{d}\left(\frac{\gamma^{1/k}}{\sqrt{s}} + \left(\frac{2}{\tau}\right)^{k-1} \gamma\right) \right)
\end{equation}

\begin{equation}
\expec\|N\|^2 = \expec\|\widetilde{\nu} - \expec \widetilde{\nu}\|^2
%= \expec \| \widetilde{\nu} - \expec \widehat{\nu}\|^2
\widetilde{\mathcal{O}}\left(\frac{\tau^2 d \ln^2(d/\delta)}{s^2 \epsilon^2} + d\left(\frac{1}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)\right).
\end{equation}

Comparing <ref> with <ref>, we see that the bias and variance of the two s <ref> and <ref> are the same up to logarithmic factors. Therefore, using <ref> as  for stochastic optimization results in the same excess risk bounds as one would get by using <ref>, up to logarithmic factors.

The proof of <ref> will require <ref>, which is due to [20]. First, we need the following notation:

For random variables $X$ and $Y$, define
\[
D^{\delta}_{\infty}(X||Y) = \sup_{S \subseteq \text{supp}(X): \mathbb{P}(X \in S) \geq \delta} \ln\left[\frac{\mathbb{P}(X \in S) - \delta}{\mathbb{P}(Y \in S)}\right].
\]

An important fact is that a randomized algorithm $\Al$ is $(\epsilon, \delta)$-DP if and only if $D^{\delta}_{\infty}(\Al(X)||\Al(X')) \leq \epsilon$ for all adjacent data sets $X \sim X'$ [29].

Let $s \in \mathbb{N}, \epsilon \leq 15$.

Let $g \geq \tau \sqrt{s}$, $b > \frac{180g^2 \ln(2/\delta)}{\epsilon^2 s}$ and
$p = \frac{90g^2 \ln(2/\delta)}{b \epsilon^2 s}$. Then for any adjacent scalar databases $Z, Z' \in [0, \tau]^r$ differing on user $u$ ($z_u \neq z'_u$), we have:

a) $D^{\delta}_{\infty}(\mathcal{S} \circ \mathcal{R}_{1D}^r(Z)
%; r; \tau; g; b; p)
||\mathcal{S} \circ \mathcal{R}_{1D}^r(Z')) \leq \epsilon\left(\frac{2}{g} + \frac{|z_u - z'_u|}{\tau}\right)$.

b) Unbiasedness: $\expec[\mathcal{P}_{1D}(Z)] = \sum_{i=1}^s z_i$.

c) Variance bound: $\expec[(\mathcal{P}_{1D}(Z) - \sum_{i=1}^s z_i)^2] = \mathcal{O}\left(\frac{\tau^2}{\epsilon^2}\ln(1/\delta)\right)$.

To prove the utility guarantees in <ref>, we begin by providing (<ref>) high probability bounds on the bias and noise induced by the non-private version of <ref>, in which $\mathcal{P}_{1D}$ in line 8 is replaced by the sum of $z \in Z_j^i$. <ref> is a refinement of <cit.>, with correct scaling for arbitrary $\gamma > 0$ and exact constants:

Let $\zeta \in (0,1)$ and $X \sim \mathcal{D}$ be a random $d$-dimensional vector with mean $\nu$ such that $\|\nu\| \leq L$ and $\expec|\langle X - \nu, e_j \rangle|^k \leq \gamma$ for some $k \geq 2$ for all $j \in [d]$. Consider the non-private version of <ref> run with $m = \left\lceil 20 \log(4d/\zeta) \right\rceil$ and $\tau \geq 2L$, where $\mathcal{P}_{1D}$ in line 8 is replaced by the sum of $z \in Z_j^i$. Denote the output of this algorithm by $\hat{\nu}$. Then with probability at least $1 - \zeta$, we have
\[
\|\hat{\nu} - \nu\| \leq 10\sqrt{d}\left(\sqrt{\frac{m}{s}}\gamma^{1/k} + \gamma\left(\frac{2}{\tau}\right)^{k-1}\right).
\]

Denote $X = (x_1, \cdots, x_d)$, $\nu = (\nu_1, \cdots, \nu_d)$,

and $z_j := \Pi_{[-\tau, \tau]}(x_j)$ for $j \in [d]$. By <ref> (stated and proved below),
\begin{equation}
\label{eq:one}
|\expec z_j - \nu_j | \leq 10 \gamma \left(\frac{2}{\tau} \right)^{k-1}.
\end{equation}
Now an application of <cit.> with $\rho(t):=
\begin{cases}
\frac{t^2}{2} &\mbox{if $t \in [-\tau, \tau]$} \\
\tau t &\mbox{if $t > \tau$}\\
-\tau t &\mbox{if $t < -\tau$}
\end{cases}$ shows that $\expec|z_j - \expec z_j |^k \leq \gamma$. Hence by <ref> (stated and proved below),
\begin{align}
\mathbb{P}\left(\Bigg|\hat{\nu}_j^i - \expec z_j\Bigg| \leq 10\sqrt{\frac{m}{s}}\gamma^{1/k}\right) &= \mathbb{P}\left(\Bigg|\frac{m}{s}\sum_{z \in Z_j^i} z - \expec z_j \Bigg| \leq 10\sqrt{\frac{m}{s}}\gamma^{1/k}\right) \nonumber \\
&\geq 0.99, ~\forall i \in [m], ~j \in [d].
\end{align}

Next, using the “median trick” (via Chernoff/Hoeffding bound, see e.g. [6]), we get
\begin{equation}
\label{eq:two}
\mathbb{P}\left(\Bigg|\hat{\nu}_j - \expec z_j\Bigg| \leq 10\sqrt{\frac{m}{s}}\gamma^{1/k}\right) \leq 2e^{-m/20} \leq \frac{\zeta}{2d},
\end{equation}

where the last inequality follows from our choice of $m$. Now, by union bound, we have for any $a > 0$ that
\begin{align}
\prob(\|\hat{\nu} - \nu\| \geq \sqrt{d} a) &\leq \sum_{j=1}^d \prob(|\hat{\nu}_j - \nu_j| \geq a) \\
&\leq \sum_{j=1}^d [\prob(|\hat{\nu}_j - \expec z_j| \geq a/2) + \prob(|\expec z_j - \nu_j| \geq a/2)].
\end{align}

Plugging in $a = 2\left[10 \gamma^{1/k}\sqrt{\frac{m}{s}} + 10 \gamma\left(\frac{2}{\tau}\right)^{k-1}\right]$ implies

\begin{equation}
\prob\left(\|\hat{\nu} - \nu\| \geq \sqrt{d}\left(10 \gamma^{1/k}\sqrt{\frac{m}{s}} + 10 \gamma\left(\frac{2}{\tau}\right)^{k-1} \right)\right) \leq d\left(\frac{\zeta}{2d} + 0\right) = \frac{\zeta}{2},
\end{equation}

by <ref> and <ref>.

Below we give the lemmas that we used in the proof of <ref>. The following is a refinement of <cit.> with the proper scaling in $\gamma$:

Let $x \sim \mathcal{D}$ be a random variable with mean $\nu$ and $\expec|x - \nu|^k \leq \gamma$ for some $k \geq 2$. Let $z = \Pi_{[-\tau, \tau]}(x)$ for $\tau \geq 2|\nu|$. Then,
\[
| \nu - \expec z | \leq 10 \gamma\left(\frac{2}{\tau}\right)^{k-1}.
\]

We begin by recalling the following form of Chebyshev's inequality:
\begin{equation}
\label{eq:cheb}
\prob(|x - \nu| \geq c) \leq \frac{\expec|x - \nu|^k}{c^k} \leq \frac{\gamma}{c^k}
\end{equation}

for any $c > 0$, via Markov's inequality. By symmetry, we may assume without loss of generality that $\nu \leq 0$. Note that
\begin{align}
\nu - \expec z &= \expec\left[x\mathbbm{1}_{x < -\tau} + x\mathbbm{1}_{x > \tau} + x\mathbbm{1}_{x \in [-\tau, \tau]} - (-\tau\mathbbm{1}_{x < -\tau} + \tau\mathbbm{1}_{x > \tau} + x\mathbbm{1}_{x \in [-\tau, \tau]})\right] \\
&= \expec\left[(x + \tau)\mathbbm{1}_{x < -\tau} + (x - \tau)\mathbbm{1}_{x > \tau}\right].
\end{align}

\begin{equation}
|\nu - \expec z| \leq \underbrace{|\expec(x + \tau) \mathbbm{1}_{x < - \tau}|}_{\textcircled{a}} + \underbrace{|\expec(x - \tau) \mathbbm{1}_{x > \tau}|}_{\textcircled{b}},
\end{equation}

by the triangle inequality. Now,
\begin{align}
\label{eq:q}
\textcircled{a} &= |\expec[x - \nu - (-\tau - \nu)\mathbbm{1}_{x < -\tau}| \nonumber \\
&\leq \expec[|x - \nu|\mathbbm{1}_{x < -\tau}] + |-\tau - \nu|\expec \mathbbm{1}_{x < -\tau} \nonumber \\
&\leq \left(\expec|x - \nu|^k \right)^{1/k}(\prob(x < -\tau))^{(k-1)/k} + |-\tau - \nu|\prob(x < -\tau),
\end{align}

by Holder's inequality. Also, since $-\frac{\tau}{2} \leq \nu \leq 0$, we have
\begin{equation}
\label{eq:pb}
\prob(x < - \tau) \leq \prob\left(x < \nu - \frac{\tau}{2}\right) \leq \prob\left(|x - \nu| > \frac{\tau}{2}\right) \leq \frac{\gamma 2^k}{\tau^k},
\end{equation}

via <ref>. Plugging <ref> into <ref> and using the bounded moment assumption, we get
\begin{align}
\textcircled{a} &\leq \gamma^{1/k} \left(\frac{\gamma 2^k}{\tau^k}\right)^{(k-1)/k} + (\tau + |\nu|)\frac{\gamma 2^k}{\tau^k}\\
&\leq \gamma\left(\frac{2}{\tau}\right)^{k-1} + 4\gamma \left(\frac{2}{\tau}\right)^{k-1} \\
&= 5\gamma \left(\frac{2}{\tau}\right)^{k-1}.
\end{align}

Likewise, a symmetric argument shows that $\textcircled{b} \leq 5\gamma \left(\frac{2}{\tau}\right)^{k-1}$. Hence the lemma follows.

Below is a is a more precise and general version of <cit.> (scaling with $\gamma$):

Let $\mathcal{D}$ be a distribution over $\mathbb{R}$ with mean $\nu$ and $\expec|\mathcal{D} - \nu|^k \leq \gamma$ for some $k \geq 2$. Let $x_1, \cdots, x_n$ be i.i.d. samples from $\mathcal{D}$. Then, with probability at least $0.99$,
\[
\Bigg|\frac{1}{n} \sum_{i=1}^n x_i - \nu\Bigg| \leq \frac{10 \gamma^{1/k}}{\sqrt{n}}.
\]

First, by Jensen's inequality, we have
\begin{equation}
\label{eq:c}
\expec[(x - \nu)^2] \leq \expec\left[|x - \nu|^k\right]^{2/k} \leq \gamma^{2/k}.
\end{equation}

\begin{align}
\expec\left[\left(\frac{1}{n}\sum_{i=1}^n x_i - \nu\right)^2\right] &= \frac{1}{n^2} \expec\left[\sum_{i=1}^n (x_i - \nu)^2\right] \\
&\leq \frac{\gamma^{2/k}}{n},
\end{align}

where we used the assumption that $\{x_i\}_{i=1}^n$ are i.i.d. and  <ref>. Thus, by Chebyshev's inequality,

\[
\prob\left(\Bigg|\frac{1}{n}\sum_{i=1}^n x_i - \nu \Bigg| \geq \frac{10 \gamma^{1/k}}{\sqrt{n}}\right) \leq \frac{1}{100}.
\]

Now, we provide the proof of <ref>:

Privacy: Let $X$ and $X'$ be adjacent data sets in $\XX^s$. Assume without loss of generality that $x_1 \neq x'_1$ and $x_l = x'_l$ for $l > 1$. By the post-processing property of DP and the fact that each sample is only processed once (due to disjoint batches) during the algorithm, it suffices to fix $i \in [m]$ and show that the composition of all $d$ invocations of $\mathcal{P}_{1D}$ (line 8) for $j \in [d]$ is $(\epsilon, \delta)$-DP. Assume without loss of generality that $Z_j := Z_j^1$ and $Z_j' := (Z_j^1)'$ contain the truncated $z_{1,j} := \Pi_{[0, 2\tau]}(x_1, j)$ and $z_{1,j}' := \Pi_{[0, 2\tau]}(x'_1, j)$, respectively. Then, by the first part of <ref>, there are choices of $b$ and $p$ such that
\begin{align*}
D^{\delta_j}_{\infty}(\mathcal{S} \circ \mathcal{R}_{1D}^r(Z_j)
%; r; \tau; g; b; p)
||\mathcal{S} \circ \mathcal{R}_{1D}^r(Z_j')) \leq \epsilon_j\left(\frac{2}{g} + \frac{|z_{1,j} - z'_{1,j}|}{2\tau}\right) \leq 2 \epsilon_j,
\end{align*}

for all $j$ by the triangle inequality, provided $g \geq \max\{2, 2\tau \sqrt{s}\}$. Hence each invocation of $\mathcal{P}_{1D}$ is $(2\epsilon_j, \delta_j)$-DP. Thus, by the advanced composition theorem <cit.>, <ref> is $(\epsilon', \delta)$-SDP, where
\begin{align}
\epsilon' &= \sum_{j=1}^d (2\epsilon_j) \left(e^{2\epsilon_j} - 1\right) + 2\sqrt{2\sum_{j=1}^d \epsilon_j^2 \ln(1/\delta_j)} \\
&\leq 8d \epsilon_1^2 + \frac{\epsilon}{2} \\
&\leq \frac{\epsilon^2}{4 \ln(1/\delta_1)} + \frac{\epsilon}{2} \\
&\leq \epsilon,
\end{align}

by our choices of $\epsilon_j, \delta_j$ and the assumption that $\epsilon \leq 8 \ln(1/\delta_1) = 8 \ln(2d/\delta)$.

Bias: Let $\hat{\nu} = (\hat{\nu}_1, \cdots, \hat{\nu}_d)$ denote the output of the non-private version of <ref> where $\mathcal{P}_{1D}$ in line 8 is replaced by the (noiseless) sum of $z \in Z_j^i$. Then,  <ref> tells us that
\[
\|\hat{\nu} - \nu\|^2 \leq 200 d\left(\frac{m}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)
\]

with probability at least $1 - \zeta$, if $m = \lceil 20 \log(4d/\zeta) \rceil$. Thus,
\begin{align*}
b^2 = \|\expec \widetilde{\nu} - \nu\|^2 &=  \|\expec \hat{\nu} - \nu\|^2 \leq \expec \|\hat{\nu} - \nu\|^2 \\
&\leq 200 d\left(\frac{m}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)(1 - \zeta) + 2\sup\left(\|\hat{\nu}\|^2 + \|\nu\|^2 \right) \zeta \\
&\leq 200 d\left(\frac{m}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)  + 4\tau^2 d\zeta,
\end{align*}

since $\tau \geq L \geq \|\nu\|$ by assumption. Then choosing
\[
\zeta = \frac{1}{\tau^2}\left[\left(\frac{1}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)\right]
\]

and noting that $\zeta$ is polynomial in all parameters (so that $m = \widetilde{\mathcal{O}}(1)$)
\begin{equation}
\label{eq:bias bound}
b^2 \leq \expec \|\hat{\nu} - \nu\|^2 = \widetilde{\mathcal{O}}\left(d\left(\frac{1}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)\right).
\end{equation}

Variance: We have

\begin{align}
\expec \|N\|^2  &\leq 2 \left[\underbrace{\expec\|\widetilde{\nu} - \hat{\nu}\|^2}_{\textcircled{a}} + \underbrace{\expec\|\hat{\nu} - \expec\hat{\nu} \|^2}_{\textcircled{b}} \right].
\end{align}

We bound a as follows. For any $j \in [d], i \in [m]$, denote  $\hat{\nu}_j^i:= \frac{m}{s} \sum_{z \in Z_j^i} z$ (c.f. line 8 of <ref>), the mean of the (noiseless) $s/m$ clipped $j$-th coordinates of data in group $i$. Denote $\widetilde{\nu}^i_j$ as in <ref>, which is the same as $\hat{\nu}_j^i$ except that the summation over $Z_j^i$ is replaced by $\mathcal{P}_{1D}(Z_j^i)$. Also, denote $\hat{\nu}_j = \text{median}(\hat{\nu}_j^1, \cdots, \hat{\nu}_j^m)$ and $\widetilde{\nu}_j =  \text{median}(\widetilde{\nu}_j^1, \cdots, \widetilde{\nu}_j^m)$. Then for any $j \in [d]$, we have
\begin{align*}
|\hat{\nu}_j - \widetilde{\nu}_j|^2 &= |\text{median}(\hat{\nu}_j^1, \cdots, \hat{\nu}_j^m) - \text{median}(\widetilde{\nu}_j^1, \cdots, \widetilde{\nu}_j^m)|^2 \\
&\leq \left(\sum_{i=1}^m |\hat{\nu}_j^i- \widetilde{\nu}_j^i| \right)^2 \\
&= \left(\sum_{i=1}^m \left|\frac{m}{s} \left[\sum_{z \in Z_j^i} z - \mathcal{P}_{1D}(Z_j^i) \right] \right| \right)^2 \\
&\leq \left(\frac{m}{s} \sum_{i=1}^m \left|\mathcal{P}_{1D}(Z_j^i) - \sum_{z \in Z_j^i} z \right| \right)^2 \\
&\leq \frac{m^2}{s^2} m \sum_{i=1}^m \left|\mathcal{P}_{1D}(Z_j^i) - \sum_{z \in Z_j^i} z \right|^2.
\end{align*}

Now using <ref> (part c), we get
\begin{align*}
\expec[|\hat{\nu}_j - \widetilde{\nu}_j|^2] &\leq \frac{m^4}{s^2} \expec\left[\left|\mathcal{P}_{1D}(Z_j^i) - \sum_{z \in Z_j^i} z \right|^2 \right] \\
&\lesssim \frac{m^4}{s^2} \frac{\tau^2}{\epsilon_j^2} \ln(1/\delta_j)  \\
&\lesssim \frac{m^4 \tau^2 d \ln^2(d/\delta)}{s^2 \epsilon^2}.
\end{align*}

Thus, summing over $j \in [d]$, we get
\begin{align}
\label{eq: a bound}
\expec[\|\hat{\nu} - \widetilde{\nu}\|^2] \lesssim \frac{m^4 \tau^2 d \ln^2(d/\delta)}{s^2 \epsilon^2}.
\end{align}

Next, we bound b:
\begin{align*}
\textcircled{b} &= \expec \| \hat{\nu} - \expec \hat{\nu}\|^2\\
&\leq 2 \left[\expec\|\hat{\nu} - \nu\|^2 +  \|\nu - \expec\hat{\nu}\|^2\right] \\
&\leq 2\left[\expec\|\hat{\nu} - \nu\|^2 + b^2\right] \\
& = \widetilde{\mathcal{O}}\left(d\left(\frac{1}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)\right),
\end{align*}

by the bias bound <ref>. Therefore,
\[
\expec\|N\|^2 =  \widetilde{\mathcal{O}}\left(\frac{\tau^2 d \ln^2(d/\delta)}{s^2 \epsilon^2} + d\left(\frac{1}{s} \gamma^{2/k} + \gamma^2 \left(\frac{2}{\tau}\right)^{2k-2} \right)\right).
\]


# Bell correlations of a thermal fully-connected spin chain in a vicinity of a
quantum critical point

D A Hamza and J Chwedeńczuk  Faculty of Physics, University of Warsaw, ul.
Pasteura 5, PL-02-093 Warsaw, Poland<EMAIL_ADDRESS>
###### Abstract

Bell correlations are among the most exotic phenomena through which quantum
mechanics manifests itself. Their presence signals that the system can violate
the postulates of local realism, once believed to be the nonnegotiable
property of the physical world. The importance of Bell correlations from this
fundamental point of view is even straightened by their applications—ranging
from quantum cryptography through quantum metrology to quantum computing.
Hence it is of growing interest to characterize the “Bell content” of complex,
scalable many-body systems.

Here we perform the detailed analysis of the character and strength of many-
body Bell correlations in interacting multi-qubit systems with particle-
exchange symmetry. Such configuration can be mapped onto an effective
Schrödinger-like equation, which allows for precise analytical predictions. We
show that in the vicinity of the quantum critical point, these correlations
quickly become so strong that only a fraction of qubits remains uncorrelated.
We also identify the threshold temperature, which, once overpassed, empowers
thermal fluctuations that destroy Bell correlations in the system.

We hope that the approach presented here, due to its universality, could be
useful for the upcoming research on genuinely nonclassical Bell-correlated
complex systems.

##  1 Introduction

The Bell’s rebuke of the Einstein-Podolsky-Rosen (EPR) call for complementing
quantum mechanics with a local and realistic theory [1] was conceived with the
least complex but still highly non-trivial many-body configuration: two
spin-$1/2$ particles [2]. Though its deceptive simplicity, for many years this
system has been driving the progress of both the theory and the experiment. In
1969, the Clauser-Horne-Shimony-Holt inequality [3] that refined the original
Bell’s formulation opened the way for first experiments showing the violation
of Bell inequalities [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. Finally,
loophole-free tests of Bell nonlocality were reported, discarding local
realism in quantum mechanics for good (apart from other possibilities like the
global realism [16], or some hard to exclude, though highly unlikely and
unreasonable scenarios, like a global conspiracy of producers of microchips or
detection equipment, intended to fool the experimentalists). The significance
of these efforts was acknowledged by the Nobel committee and the prize was
awarded in 2022.

With the much descried arrival of the noisy intermediate-scale quantum (NISQ)
devices, the research interest has shifted towards more complex systems. It
has become of growing relevance to understand how the fundamental quantum
relations, such as the entanglement [17, 18], the EPR steering [1, 19, 20] and
the Bell nonlocality can be created and detected, in particular in multi-qubit
systems.

The problem of detecting many-body Bell correlations, which are the main focus
of this work, has been formulated for $N$-qubit systems using $N$-body
correlation functions in various ways, most notably by means of Mermin-Bell
inequalities [21], and finally expressed in the all-encompassing form [22].
This general formulation has different variants (i.e., forms of correlation
functions, adapted for a concrete physical case) and one particularly useful
is the one formulated in Refs. [23, 24, 25]. It enabled to establish the link
between the Bell nonlocality and quantum-enhanced metrology [26] or determine
the strength of many-body Bell correlations in a commonly used method of
creating atomic spin-squeezed states [27, 28], i.e., the one-axis twisting
method [29, 30, 31, 32, 33, 34, 35].

Here we scrutinize the emergence of Bell correlations in the versatile and
scalable multi-qubit system. It is a collection of $N$ qubits where each pair
interacts with the same strength. Such setup can be realized using different
platforms well-suited for quantum technologies, such as the fully-connected
Ising spin-chain [36, 37, 38, 39, 40, 41, 42] or alternatively by trapping the
Bose-Einstein condensate in a double-well potential. While throughout this
work we will refer to the latter case, it is equivalent, from the perspective
of problems considered here, to the former. In this scenario, the separated
single-body states localized around the two minima of the trap play the role
of the pair of single-qubit levels. In the tight-binding regime, this system
is described by the Bose-Hubbard (BH) Hamiltonian, where the on-site two-body
interactions compete with the coherent Josephson tunneling across the barrier
of the double-well potential [43, 44, 45, 46, 47, 48, 49]. On the attractive
side of the interaction, there exists a quantum phase transition (QPT) [50,
51, 52, 53, 54, 55, 56, 57, 58, 59, 60], upon passing of which the properties
of the BH Hamiltonian suddenly change—from a gaussian state to the macroscopic
superposition [52, 53, 54, 55, 56, 57]. We demonstrate that in the vicinity of
this point, the Bell correlations are genuinely many-body, i.e., a macroscopic
fraction of qubits is Bell-correlated. To show this, we derive a simple yet
informative formula, that allows to predict the strength of Bell correlations
related to this emergent macroscopic superposition. Next, we consider non-zero
temperatures and analytically demonstrate that the relevant temperature scale,
above which the thermal fluctuations “kill” the Bell nonlocality, is related
to the energy gap between the ground- and the first-excited state. This allows
to derive a clear criterion for the critical—from the point of view of many-
body Bell correlations—temperature.

We hope that the presented analysis will contribute to the emerging field of
the NISQ devices, allowing for precise characterization of multi-qubit states
and for the accurate planning of future experiments.

##  2 Many-Body Quantum Systems

Consider a collection of $N$ qubits (spin-$1/2$ particles) fully connected via
the distance-independent two-body interaction of strength $U$. The spins are
subject to an external uniform magnetic field aligned by the $x$-axis, which
amplitude is proportional to $\Omega$. One possibility of modelling such
system is via the Ising model, the one on which we focus in this work. The
Hamiltonian reads

$\displaystyle\hat{H}=-\Omega\sum_{i=1}^{N}\hat{\sigma}_{x}^{(i)}+\frac{U}{2}\sum^{N}_{i\neq
j=1}\hat{\sigma}_{z}^{(i)}\hat{\sigma}_{z}^{(j)}.$ (1)

Here $\hat{\sigma}_{\xi}^{(k)}$ is the $\xi$-component of the triad of Pauli
matrices for the $k$-th spin (i.e., $\xi=x,y,z$). Note that this Hamiltonian
is invariant under the exchange of any pair of spins, hence all its eigen-
states posses this symmetry. Thus to analyze the many-body properties of the
ground- and the thermal-states, it is convenient to map it onto the bosonic BH
Hamiltonian in the tight-binding limit, namely

$\displaystyle\hat{H}_{\rm bh}=-\Omega\hat{J_{x}}+U\hat{J_{z}}^{2},$ (2)

where the collective spin operators are

$\displaystyle\hat{J}_{\xi}=\frac{1}{2}\sum_{k=1}^{N}\hat{\sigma}_{\xi}^{(k)}.$
(3)

The spectrum of this operator is spanned by the symmetrized states
$\left|N-n,n\right\rangle$ of $n$ spins in the $-1$ eigen-state of (for
instance) the $z$-component Pauli operators, and the remaining $N-n$, in the
$+1$ eigen-state. Hence the dimensionality of the Hilbert space is:
$\dim\mathcal{H}^{({\rm sym})}_{N}=N+1$ rather than:
$\dim\mathcal{H}_{N}=2^{N}$ for a generic collection of $N$ spins. Therefore
any state vector $\left|\psi\right\rangle$ in this symmetric subspace can be
expressed as

$\displaystyle\left|\psi\right\rangle=\sum_{n=0}^{N}\psi_{n}\left|N-n,n\right\rangle,\
\ \ \mathrm{with}\ \ \ \sum_{n=0}^{N}\left|\psi_{n}\right|^{2}=1.$ (4)

To extract the crucial properties of the eigen-states of the Hamiltonian from
Eq. (2) we follow the steps of Ref. [61]. First, we project the stationary
Schrödinger equaiton onto the $n$’th element of the basis, i.e, $\langle
N-n,n|\hat{H}\left|\psi\right\rangle=E\psi_{n}$, giving

$\displaystyle-\frac{\Omega}{2}\left[\psi_{n+1}\sqrt{(N-n)(n+1)}+\psi_{n-1}\sqrt{(N-n+1)n}\right]+\dfrac{U}{4}\left(N-2n\right)^{2}\psi_{n}=E\psi_{n}$
(5)

Next, we introduce the normalized population imbalance,

$\displaystyle z_{n}=\frac{(N-n)-n}{N}=1-\frac{2n}{N}$ (6)

which varies from $-1$ to $1$ with the increment equal to $\Delta z=2/N$.

Figure 1: The effective potential $V_{\rm eff}(z)$ (solid black lines) for
$N=100$ and $\gamma=-1.1$ (left) and $\gamma=-1.5$ (right). The dashed gray
lines show the ground state for each case.

The Eq. (5), expressed in terms of this new variable becomes

$\displaystyle-\frac{\Omega
N}{2}\left[\psi_{n+1}f_{+}(z_{n})+\psi_{n-1}f_{-}(z_{n})\right]+\frac{UN^{2}}{4}\psi_{n}z^{2}_{n}=E\psi_{n},$
(7) $\displaystyle\mathrm{where}\ \ \ f_{\pm}(z_{n})=\sqrt{\frac{1\pm
z_{n}}{2}\left(\frac{1\mp z_{n}}{2}+\frac{1}{N}\right)}.$

When $N$ is large, the $1/N$ term in $f_{\pm}$ can be neglected. Furthermore,
as the increment of the $z_{n}$ diminishes, the discrete variable can be
approximated with a continuous $z$. In particular, this means that the finite
difference tends to

$\displaystyle\frac{\psi_{n+1}+\psi_{n-1}-2\psi_{n}}{(\Delta
z)^{2}}\approx\frac{d^{2}}{dz^{2}}\psi(z)$ (8)

In such large-$N$ regime we obtain a stationary 1D Schrödinger-like equation
for a fictious unit-mass particle subject to an external potential

$\displaystyle V_{\rm eff}(z)=-\sqrt{1-z^{2}}+z^{2}\gamma/2,$ (9)

with $\gamma=UN/\Omega$, i.e.,

$\displaystyle\left(-\frac{2}{N^{2}}\sqrt{1-z^{2}}\frac{d^{2}}{dz^{2}}+V_{\rm
eff}(z)\right)\psi(z)=\tilde{E}\psi(z).$ (10)

The normalized energy is $\tilde{E}=\frac{2E}{\Omega N}$. Note also that here
$2/N$ is the dimensionless equivalent of the reduced Planck constant $\hbar$.
For the mathematical details of this derivation with the discussion of all its
limitations, we refer the readers to Ref. [61].

The ground state of Eq. (10) depends crucially on the value of $\gamma$, hence
on the ratio (and the sign) of the interaction-to-tunneling energies. For
large and positive $\gamma$’s the $V_{\rm eff}(z)\eqsim\frac{1}{2}\gamma
z^{2}$, hence the problem simplifies to that of the 1D harmonic potential. As
$\gamma$ grows, the width of the ground state Gaussian shrinks. This implies
the diminishing fluctuations of the population imbalance, that, since the
intra-mode coherence is mostly maintained, signals the spin-squeezing of the
state [27, 28].

Figure 2: The comparison of the potentials $V_{\rm eff}(z)$ (normalized to
$V_{\rm min}=-1$) from Eq. (9) and from Eq. (11) for $N=500$ and
$\gamma=-1.4$. The horizontal dashed blue line shows the energy scales of the
few lowest-lying states of either potential.

On the negative-$\gamma$ side, the shape of the ground state is more
diversified, depending on the value of the parameter. For
$\gamma\gtrsim\gamma_{0}\equiv-1$ the wave-function is still a Gaussian, but
now its width is growing as $\gamma$ approaches $\gamma_{0}$—this in turn
signals the phase-squeezing of the sample [62]. However, when the $\gamma_{0}$
is crossed, $V_{\rm eff}(z)$ changes abruptly indicating the passage through
the quantum phase transition (QPT). This is visualized in Fig. 1 where we show
how the potential and the ground state change across this highlighted point.
Upon crossing $\gamma_{0}$, the potential breaks and develops two-minima at
positions $\pm z_{0}$ with $z_{0}=(1-1/\gamma^{2})^{1/2}$. A quantum state
that respects the left/right symmetry of the problem is a macroscopic
superposition of two separated wave-packets. The fluctuations of the
population imbalance rapidly grow and as $\gamma\rightarrow-\infty$ the maxima
of $\psi(z)$ separate tending towards the NOON state, which in the ket
notation used in Eq. (4) is
$\left|\psi\right\rangle=1/\sqrt{2}\big{(}\left|N,0\right\rangle+\left|0,N\right\rangle\big{)}$.

The shape of the potential after the breaking suggests that some of the
properties of the state vector can be extracted by locally approximating the
$V_{\rm eff}(z)$ with two harmonic oscillators located at $\pm z_{0}$, namely
(recall that the fictitious particle has a mass set to unity)

$\displaystyle V_{\rm eff}(z)\simeq\frac{1}{2}\omega^{2}(z\pm
z_{0})^{2}+V_{0},$ (11)

where $\omega=\sqrt{\gamma(1-\gamma^{2})}$ and
$V_{0}=\frac{\gamma^{2}+1}{\gamma}$. In Fig. 2 we show how this approximation
works for $N=500$ qubits and $\gamma=-1.4$. Roughly speaking, the position of
the lowest-lying levels, denoted by the horizontal blue dashed line seems to
be deep in the regime where the approximation works. Nevertheless, this must
be analyzed in more detail, as discussed below.

Note that the harmonic potential (11) is an approximation of an already
approximate model of the Schrödinger-like equation from Eq. (10). In order to
make sure that in this process the errors do not accumulate to some
unacceptable level, we perform the exact diagonalization of the BH Hamiltonian
from Eq. (2) and compare its two lowest-lying energy levels with those coming
from the harmonic approximation (HA), see Eq. (11). In Fig. 3 we show the
normalized energy difference between these two outcomes, expressed in percent
units, namely

$\displaystyle\Delta
E_{i}=\left|\frac{E_{i}^{\rm(BH)}-E_{i}^{\rm(HA)}}{E_{i}^{\rm(BH)}}\right|\times
100\%$ (12)

for $i=0$ (ground state) and $i=1$ (first excited state). The left panel shows
the difference for two $\gamma$’s equal to -1.1 and -1.5 as a function of $N$.
As expected from the procedure described above, the discrepancy diminishes as
$N$ grows and for $N=500$ can be safely kept below 1%. Similarly, the
variation of $\gamma$ for fixed $N$’s (100 and 500), see the right panel,
confirms the satisfactory precision of the harmonic approximation.

Figure 3: Left: The normalized energy difference, $\Delta E_{i}$ see Eq.(12),
for $i=0$ (the ground state, solid lines) and $i=1$ (the first excited state,
dashed lines) of the potentials (9) and (11) for $\gamma=-1.1$ (black) and
$\gamma=-1.5$ (gray) as a function of $N$. Right: Here, $N$ is fixed to
$N=100$ (black) and $N=500$ (gray), while $\gamma$ is varied.

Last but not least, we scrutinize the quality of the harmonic approximation by
calculating the overlap between the eigen-states corresponding to these two
energy levels with the eigen-states of the exact BH Hamiltonian

$\displaystyle\mathcal{F}_{i}=\left|\left\langle\psi_{i}^{({\rm
BH})}|\psi_{i}^{({\rm HA})}\right\rangle\right|^{2}\times 100\%\ ,$ (13)

again with $i=0,1$. In Fig. 4 we show that for the ground state, the fidelity
$\mathcal{F}_{0}$ is above 90% for a wide range of $\gamma$’s and $N$’s. While
the fidelity $\mathcal{F}_{1}$ can be unsatisfactory for smaller $N$’s and in
the direct vicinity of the QPT, we will argue that even in $T\neq 0$ this is
not a concern—the Bell correlator we discuss below is mostly characterized by
the properties of the ground-state.

In this section we have shown how to replace to exact quantum BH model with an
approximate twin-harmonic approximation. It is valid at the negative side of
the quantum critical point and will be our working horse for the derivation of
simple yet powerful analytical formulæ of the many-body Bell correlations.

Figure 4: The quality of the harmonic approximation quantified by the fidelity
$\mathcal{F}_{i}$ of the ground- and first-excited state, expressed in
percent, see Eq. (13) The left panel displays $\mathcal{F}_{i}$ as a function
of the number of particles $N$ for $\gamma=-1.1$ (black) and $\gamma=-1.5$
(gray). The right panel is the reverse: here $N$’s are fixed to $100$ (black)
and 500 (gray) while $\gamma$ changes.

##  3 Many-body Bell correlator

We now briefly review the theory behind the method of detecting many-body Bell
correlators that will be used in this work. For an extensive discussion of its
properties, see References [63, 24, 64, 65, 66].

Consider $m$ objects, each being a subject of local measurements of two binary
quantities $\sigma_{x}^{(k)}=\pm 1$ and $\sigma_{y}^{(k)}=\pm 1$ with
$k=1\ldots m$. These outcomes are combined locally (for each object) into
$\sigma_{+}^{(k)}=\frac{1}{2}(\sigma_{x}^{(k)}+i\sigma_{y}^{(k)})$ and the
following correlator is constructed

$\displaystyle\mathcal{E}_{m}=\left|\langle\sigma_{+}^{(1)}\ldots\sigma_{+}^{(m)}\rangle\right|^{2}.$
(14)

Here, the $\langle\cdot\rangle$ denotes averaging over experimental
repetitions. If the $\mathcal{E}_{m}$ can be reproduced by system that is
consistent with the local hidden variable theory [1, 2, 67, 68], then this
average can be expressed in terms of an integral over this (possibly multi-
variate) variable $\lambda$ distributed with a probability density
$p(\lambda)$, namely

$\displaystyle\mathcal{E}_{m}=\left|\int\\!\\!d\lambda\,\,p(\lambda)\,\sigma_{+}^{(1)}(\lambda)\ldots\sigma_{+}^{(m)}(\lambda)\right|^{2}.$
(15)

Using the Cauchy-Schwarz inequality for complex integrals we obtain

$\displaystyle\mathcal{E}_{m}\leqslant\int\\!\\!d\lambda\,\,p(\lambda)\,\left|\sigma_{+}^{(1)}(\lambda)\right|^{2}\ldots\left|\sigma_{+}^{(m)}(\lambda)\right|^{2}=2^{-m}.$
(16)

Thus $\mathcal{E}_{m}\leqslant 2^{-m}$ is the $m$-body Bell inequality.

This inequality can be tested with quantum systems of $m$ qubits. In such
case, the correlator $\mathcal{E}_{m}$ from Eq. (14) is replaced by its
quantum-mechanical ($q$) equivalent

$\displaystyle\mathcal{E}^{(q)}_{m}=\left|{\rm
Tr}\left[\hat{\varrho}\,\hat{\sigma}_{+}^{(1)}\otimes\ldots\otimes\sigma_{+}^{(m)}\right]\right|^{2}.$
(17)

For the case of $m$-qubit Greenberg-Horne-Zeilinger (GHZ) state [69], i.e.,

$\displaystyle\left|\psi\right\rangle=\frac{1}{\sqrt{2}}(\left|\uparrow\right\rangle^{\otimes
m}+\left|\downarrow\right\rangle^{\otimes m}),\ \ \ \mathrm{where}\
\hat{\sigma}_{+}\left|\downarrow\right\rangle=\left|\uparrow\right\rangle,$
(18)

we obtain $\mathcal{E}^{(q)}_{m}=1/4$, hence the exponential (as a function of
$m$) breaking of the bound (16).

This quantum correlator (17) can be adapted to symmetric systems, where only
collective operations are allowed. In this case, the single-body particle-
resolving operators $\hat{\sigma}_{+}^{(k)}$ are replaced with its (still one-
body) collective equivalent

$\displaystyle\hat{\sigma}_{+}^{(k)}\longrightarrow\hat{J}_{+}=\sum_{k=1}^{N}\hat{\sigma}_{+}^{(k)},$
(19)

in analogy to Eq. (3). Upon replacing the $\sigma_{+}^{(k)}$’s in Eq. (17) by
$\hat{J}_{+}^{m}$ one notices that there are $N!/(N-m)!$ more terms in the
latter case, the consequence of the permutational invariance of the problem.
Hence the proper symmetric (the symmetrization is here denoted by the tilde
symbol) $m$-body Bell inequality is

$\displaystyle\tilde{\mathcal{E}}^{(q)}_{m}\equiv\left|\langle\hat{J}^{m}_{+}\rangle\right|^{2}\leq\left(\frac{N!}{(N-m)!}\right)^{2}2^{-m}$
(20)

and more details on its derivation can be found in Ref. [66]. It is convenient
to normalize the l.h.s. of this inequality by its r.h.s. and introduce [70]

$\displaystyle
Q_{m}=\log_{2}\left[2^{m}\left(\frac{N!}{(N-m)!}\right)^{-2}\tilde{\mathcal{E}}^{(q)}_{m}\right].$
(21)

The Bell inequality then becomes

$\displaystyle Q_{m}\leqslant 0.$ (22)

It is the purpose of the remainder of this work to demonstrate that both at
zero temperature ($T=0$) as well as when $T>0$, there exists an order $m=\mu$
of the Bell correlator which is significantly above the nonlocality bound
$Q_{\mu}>0$ in the vicinity of the QPT. This particular $Q_{\mu}$, as we argue
below, informs that an extensive (i.e., growing linearly with $N$) number of
qubits is Bell-correlated. Moreover, $Q_{\mu}$ turns out to be directly linked
to how far $\gamma$ is from the quantum critical point. Benefiting from the
above-verified harmonic approximation we will connect the correlation order
$\mu$ and the value of the corresponding correlator, $Q_{\mu}$, with $\gamma$
using simple but versatile analytical formulæ, and extend the analysis to non-
zero temperatures.

The concept relies on the following observation. If the operator
$\hat{J}^{m}_{+}$ that determines the correlator $Q_{m}$ acts on a ket
$\left|N-n,n\right\rangle$ it gives

$\displaystyle\hat{J}^{m}_{+}\left|N-n,n\right\rangle=j_{nm}\left|N-n+m,n-m\right\rangle\
\ \ \mathrm{with}\ \ \ j_{nm}=m!\sqrt{\binom{n}{m}\binom{N-n+m}{m}}.$ (23)

Hence the average of this operator is equal to

$\displaystyle\langle\hat{J}^{m}_{+}\rangle=\sum_{n=0}^{N-m}\varrho_{n,n+m}j_{nm},$
(24)

i.e., a coherent sum of all the elements of the density matrix distanced by
$m$ from the diagonal and weighted with $j_{nm}$’s. Here we used the general
expression for the density matrix in the $N$-qubit symmetric subspace, namely

$\displaystyle\hat{\varrho}=\sum_{nn^{\prime}=0}^{N}\varrho_{nn^{\prime}}|N-n,n\rangle\\!\langle
N-n^{\prime},n^{\prime}|.$ (25)

Due to the characteristic twin-peak structure of the ground state, one $m$
stands out, namely such that is directly linked to the separation of the two
maxima $\delta z=2z_{0}$ by means of Eq. (6), i.e.,

$\displaystyle\mu=\frac{N}{2}\lceil\delta z\rceil=N\lceil z_{0}\rceil,$ (26)

see Fig. 5. Here, $\lceil x\rceil$ denotes the integer, which is closest and
bigger equal to $x$ (the “ceil” function). Hence the conjecture: one
particular order of Bell correlator, at a given $\gamma$, stands out, i.e.,
such for which $m=\mu$ holds. Moreover, its value can be reproduced with
excellent precision by reducing the sum in Eq. (24) merely to the contribution
coming from the peak values at these two maximæ positioned at
$n_{\pm}=(N\pm\mu)/2$. In other words, using Eqs (20) and (23), the correlator
is approximated by

$\displaystyle
Q_{\mu}\simeq\log_{2}\left[2^{\mu}\left(\frac{N!}{(N-\mu)!}\right)^{-2}\left|\varrho_{n_{+},n_{-}}j_{n_{+},n_{-}}\right|^{2}\right]$
(27)

We now show, that this very simple formula works exceptionally well for a pure
ground state and allows to determine the critical temperature, above which the
Bell correlator $Q_{\mu}$ drops below the Bell limit.

It should be underlined that $Q_{\mu}$ does not provide a complete information
about Bell nonlocality in the system. At fixed $\gamma$, there are various
correlation orders that give $Q_{m}>0$. Also, there are other ways to
construct the many-body Bell correlator [22], that in principle could yield
more information about the nonlocal correlations in this case. Nevertheless
the analysis presented here has two strong points—it is simple, leading to
analytical predictions, and it allows to lower-bound the extend of Bell
correlations over the system, an important issue for many applicational
aspects, such as quantum-enhanced metrology [26].

Figure 5: The ground state of the Bose-Hubbard Hamiltonian and the effective
potential $V_{\rm eff}(z)$ (both in arbitrary units picked only for
illustration) for $N=500$ qubits and $\gamma=-1.4$. As argued in the text, the
distance between the peaks determines the relevant order of the Bell
correlator.

##  4 Bell correlations at $T=0$

The Harmonic approximation yields the ground state of the system to be

$\displaystyle\left|\psi\right\rangle\simeq\left(\frac{\omega}{2\pi
N}\right)^{\frac{1}{4}}\sum_{n=0}^{N}\left(e^{-(\frac{N-\mu}{2}-n)^{2}\frac{\omega}{N}}+e^{-(\frac{N+\mu}{2}-n)^{2}\frac{\omega}{N}}\right)\left|N-n,n\right\rangle.$
(28)

Hence by taking the peak-values of the state coefficients [i.e., located at
$n_{\pm}=(N\pm\mu)/2$], and substituting this result into Eq. (27), we obtain

$\displaystyle
Q_{\mu}=\log_{2}\left[2^{\mu}\binom{\frac{N+\mu}{2}}{\mu}^{2}\binom{N}{\mu}^{-2}\frac{\omega}{2\pi
N}\right].$ (29)

Keeping the dominant terms that scale with $N$, we obtain a simple but
powerful expression that allows to lower-bound the strength of $\mu$-body Bell
correlations in the vicinity of the quantum critical point, namely

$\displaystyle
Q_{\mu}\simeq\mu\log_{2}\left[\frac{N+\mu}{N-\mu}\right]-\mu-N\log_{2}(\gamma^{2}).$
(30)

Finally, by using Eq. (26), we obtain a compact formula

$\displaystyle Q_{\mu}=Nf(\gamma),\ \ \ \mathrm{where}\ \ \
f(\gamma)=\frac{\sqrt{\gamma^{2}-1}}{|\gamma|}\left[2\log_{2}\left(\sqrt{\gamma^{2}-1}+|\gamma|\right)-1\right]-2\log_{2}|\gamma|.$
(31)

Most importantly, the correlator is extensive in $N$, which has profound
consequences for the quantitative characterization of many-body Bell
nonlocality. Moreover, this form of $Q_{\mu}$ implies that the value of
$\gamma$ at which the many-body Bell correlations emerge is universal (i.e.,
independent of $N$) and equal to $\gamma_{0}\simeq-1.3$.

Note that the value of $Q_{\mu}$ carries information about the nonlocality
depth, i.e., about how many qubits are Bell-correlated. And so, if

$\displaystyle\mu-2-(k+1)<Q_{\mu}\leqslant\mu-2-k,$ (32)

where $k\in\mathbb{N}$, then this correlator can be reproduced with a state of
$m$ qubits, where not more that $k-2$ of them are not Bell-correlated with the
rest [65, 71, 72]. By using Eq. (31) we notice that the depth of Bell
correlations in the ground-state can be lower-bounded as

$\displaystyle k\leqslant N(\lceil z_{0}\rceil-f(\gamma))-2.$ (33)

Hence we conclude that the number of Bell correlated qubits is an extensive
function of $N$, starting from the critical point $\gamma=\gamma_{0}$. This
confirms that $Q_{\mu}$ is a useful tool to lower-bound the nonlocality depth
in this system.

##  5 Quantum Correlations at T>0

Next, we calculate the Bell correlator and the nonlocality depth in non-zero
temperatures, considering the thermal density matrix

$\displaystyle\hat{\varrho}=\frac{1}{\mathcal{Z}}\sum_{n}|\psi_{n}\rangle\\!\langle\psi_{n}|e^{-E_{n}\beta},$
(34)

where, $\beta^{-1}=k_{b}T$, $k_{b}$ is the Boltzmann constant, $\mathcal{Z}$
is the statistical sum and $\hat{H}_{\rm
bh}\left|\psi_{n}\right\rangle=E_{n}\left|\psi_{n}\right\rangle$ with the
Hamiltonian from Eq. (2).

The characteristic temperature scales are set by the eigen-energies (with
respect to the ground state) of the Hamiltonian, hence first we plot a part of
the spectrum for $\gamma$’s of interest and $N=500$ atoms. see Fig. 6. The
energy gaps $\delta_{i}=E_{i}-E_{0}$ of the seven lowest-lying excited states
are all of the order of unity, apart from the first one, i.e., $\delta_{1}$,
which quickly converges to 0 upon passing the QPT point. This observation
allows us to identify the relevant scale of $\beta$’s.

Figure 6: The energy spectrum (8 lowest levels) of the Hamiltonian $H_{\rm
bh}$ for $N=500$ atoms in the vicinity of the critical point. The left panel
shows the energies $E_{i}$ ($i=0\ldots 7$) while the right panel displays the
seven energy gaps with respect to the ground state $E_{0}$.

A following toy model indicates that the temperature should be kept well below
this smallest gap, $\delta_{1}$, to retain high values of $Q_{\mu}$. Let’s
approximate the twin-peak ground state with what is sometimes called a
“Schrödinger kitten” state

$\displaystyle\left|\psi_{0}\right\rangle=\frac{1}{\sqrt{2}}\left(\left|\frac{N+\mu}{2},\frac{N-\mu}{2}\right\rangle+\left|\frac{N-\mu}{2},\frac{N+\mu}{2}\right\rangle\right).$
(35)

The gap $\delta_{1}$ dropping to zero is the result of the fact that this
state is almost degenerate with

$\displaystyle\left|\psi_{1}\right\rangle=\frac{1}{\sqrt{2}}\left(\left|\frac{N+\mu}{2},\frac{N-\mu}{2}\right\rangle-\left|\frac{N-\mu}{2},\frac{N+\mu}{2}\right\rangle\right).$
(36)

This is just an approximation of the true ground and first-excited states at
finite $\gamma$. The exact degeneracy is when $\gamma\rightarrow\infty$ (i.e.,
when $\mu\rightarrow N$). Nevertheless, when the temperature is kept
sufficiently low so that only $\left|\psi_{0}\right\rangle$ and
$\left|\psi_{1}\right\rangle$ are populated, the thermal density matrix is
approximately

$\displaystyle\hat{\varrho}\simeq\frac{1}{1+e^{-\delta_{1}\beta}}\left(|\psi_{0}\rangle\\!\langle\psi_{0}|+|\psi_{1}\rangle\\!\langle\psi_{1}|e^{-\delta_{1}\beta}\right),$
(37)

Substituting this $\hat{\varrho}$ into Eq. (21) gives

$\displaystyle
Q_{\mu}=\log_{2}\left(2^{\mu}\left|\frac{1}{2}\binom{N}{\mu}^{-1}\binom{\frac{N+\mu}{2}}{\mu}\frac{1-e^{-\delta_{1}\beta}}{1+e^{-\delta_{1}\beta}}\right|^{2}\right).$
(38)

Clearly, when $\beta\ll\delta_{1}$, i.e., when $k_{b}T$ is much larger than
the gap, the correlator vanishes. This conjecture is fully confirmed by the
exact diagonalization, that gives the thermal state and the correlator
$Q_{\mu}$, see the right panel of Fig 7. Here. the temperature is picked to be
equal to 10% of the energy gap $\delta_{1}$, value of which is calculated
either using $\gamma=-1.1$ (solid black line) or $\gamma=-1.6$ (solid gray
line). The correlator $Q_{\mu}$ drops drastically when the interaction
strength $\gamma$ is far away from the QPT so that the population of the
excited state becomes significant.

Figure 7: The Bell correlator calculated with the full Bose-Hubbard
Hamiltonian (solid lines) and with the approximate formulæ (dashed lines).
Left: $T=0$ for $N=100$ (gray) and $N=500$ (black) qubits. The approximate
formula is in Eq. (31). Right: the correlator $Q_{\mu}$ calculated with the
thermal state, see Eq. (34) with $N=100$ qubits. The black solid line is for
$k_{b}T_{1}=10\%\times\delta_{1}$, which is the gap taken at $\gamma=-1.1$.
The gray solid line is for $k_{b}T_{2}=10\%\times\delta_{1}$, where the gap is
for $\gamma=-1.6$. Hence $T_{2}>T_{1}$. The vertical dashed lines and red dots
indicate those values of $\gamma$. Clearly, soon upon passing these points,
when the temperature becomes comparable to the corresponding gap, the
correlator drastically drops.

Recall that the many-body Bell correlations detected by $Q_{\mu}$ require
$\gamma<-1.3$, see the discussion below Eq. (31). This observation, together
with the results from the above paragraph, yield the critical temperature,
above which the Bell correlations will not be detected by $Q_{\mu}$. Namely,
if on one hand, $\gamma$ must be below $\gamma_{0}$ and the temperature must
be at least one order of magnitude smaller than the energy gap, the maximal
temperature at which $Q_{\mu}>0$ is set by a fraction of the energy gap
$\delta_{1}$ at $\gamma_{0}$. Since this gap rapidly shrinks with growing $N$,
we conclude that the many-body Bell correlator $Q_{\mu}$ is the more sensitive
to any thermal excitations, the higher $N$ is.

##  6 Conclusion

In this work we have analyzed the character of the many-body Bell correlations
in an interacting multi-qubit systems with particle-exchange symmetry, such as
the fully connected Ising model in a perpendicular magnetic field or the
tight-biding Bose-Hubbard Hamiltonian of a Bose-Einstein condensate in a
double-well potential. Our analysis resides on the observation that such
systems can be mapped onto an effective Schrödinger-like equation. This allows
for detailed analytical calculations, using a harmonic approximation of the
local minima of the effective potential. As our main result, we show that upon
passing a quantum phase transition, the ground state of the system is
characterized by strong Bell correlations that are extensive in $N$. We also
show the existence of the threshold temperature, above which, once the thermal
noise dampens the Bell correlations.

In the expected advent to NISQ devises, such precisely tailored analysis might
help in solving various theoretical and experimental problems.

## Acknowledgements

This project was funded by the National Science Centre, Poland, within the
QuantERA II Programme that has received funding from the European Union’s
Horizon 2020 research and innovation programme under Grant Agreement No
101017733, Project No. 2021/03/Y/ST2/00195.

## ORCID IDs

Danish Ali Hamza: https://orcid.org/0000-0001-9829-8260
Jan Chwedeńczuk: https://orcid.org/0000-0002-9250-4227

## References

## References

  * [1] Einstein A, Podolsky B and Rosen N 1935 Phys. Rev. 47(10) 777–780
  * [2] Bell J S 1964 Physics 1 195
  * [3] Clauser J F, Horne M A, Shimony A and Holt R A 1969 Phys. Rev. Lett. 23(15) 880–884
  * [4] Freedman S J and Clauser J F 1972 Phys. Rev. Lett. 28(14) 938–941
  * [5] Aspect A, Grangier P and Roger G 1981 Phys. Rev. Lett. 47(7) 460–463
  * [6] Aspect A, Dalibard J and Roger G 1982 Phys. Rev. Lett. 49(25) 1804–1807
  * [7] Tittel W, Brendel J, Gisin B, Herzog T, Zbinden H and Gisin N 1998 Phys. Rev. A 57(5) 3229–3232
  * [8] Tittel W, Brendel J, Zbinden H and Gisin N 1998 Phys. Rev. Lett. 81(17) 3563–3566
  * [9] Weihs G, Jennewein T, Simon C, Weinfurter H and Zeilinger A 1998 Phys. Rev. Lett. 81(23) 5039–5043
  * [10] Pan J W, Bouwmeester D, Daniell M, Weinfurter H and Zeilinger A 2000 Nature 403 515–519
  * [11] Kielpinski D, Meyer V, Sackett C A, Itano W M, Monroe C and Wineland D J 2001 Nature 409 791–794
  * [12] Gröblacher S, Paterek T, Kaltenbaek R, Brukner Č, Żukowski M, Aspelmeyer M and Zeilinger A 2007 Nature 446 871–875
  * [13] Salart D, Baas A, van Houwelingen J A W, Gisin N and Zbinden H 2008 Phys. Rev. Lett. 100(22) 220404
  * [14] Ansmann M, Wang H, Bialczak R C, Hofheinz M, Lucero E, Neeley M, O’Connell A D, Sank D, Weides M, Wenner J, Cleland A N and Martinis J M 2009 Nature 461 504–506
  * [15] Giustina M, Mech A, Ramelow S, Wittmann B, Kofler J, Beyer J, Lita A, Calkins B, Gerrits T, Nam S W, Ursin R and Zeilinger A 2013 Nature 497 227–230
  * [16] Hensen B, Bernien H, Dréau A, Reiserer A, Kalb N, Blok M, Ruitenberg J, Vermeulen R, Schouten R, Abellán C et al. 2015 Nature 526 682–686
  * [17] Schrödinger E 1935 Naturwissenschaften 23 807–812
  * [18] Horodecki R, Horodecki P, Horodecki M and Horodecki K 2009 Rev. Mod. Phys. 81(2) 865–942 URL https://link.aps.org/doi/10.1103/RevModPhys.81.865
  * [19] Reid M D, Drummond P D, Bowen W P, Cavalcanti E G, Lam P K, Bachor H A, Andersen U L and Leuchs G 2009 Rev. Mod. Phys. 81 1727–1751
  * [20] Yadin B, Fadel M and Gessner M 2021 Nature communications 12 1–8
  * [21] Mermin N D 1990 Phys. Rev. Lett. 65 1838
  * [22] Żukowski M and Brukner Č 2002 Phys. Rev. Lett. 88 210401
  * [23] He Q, Cavalcanti E, Reid M and Drummond P 2010 Phys. Rev. A 81 062106
  * [24] He Q, Drummond P and Reid M 2011 Phys. Rev. A 83 032120
  * [25] Cavalcanti E G, Jones S J, Wiseman H M and Reid M D 2009 Phys. Rev. A 80 032112
  * [26] Niezgoda A and Chwedeńczuk J 2021 Phys. Rev. Lett. 126(21) 210506 URL https://link.aps.org/doi/10.1103/PhysRevLett.126.210506
  * [27] Kitagawa M and Ueda M 1993 Phys. Rev. A 47 5138–5143
  * [28] Wineland D, Bollinger J, Itano W and Heinzen D 1994 Phys. Rev. A 50 67
  * [29] Esteve J, Gross C, Weller A, Giovanazzi S and Oberthaler M 2008 Nature 455 1216–1219
  * [30] Appel J, Windpassinger P J, Oblak D, Hoff U B, Kjærgaard N and Polzik E S 2009 PNAS 106 10960–10965
  * [31] Pezzé L and Smerzi A 2009 Phys. Rev. Lett. 102 100401
  * [32] Riedel M F, Böhi P, Li Y, Hänsch T W, Sinatra A and Treutlein P 2010 Nature 464 1170–1173 ISSN 1476-4687 URL http://dx.doi.org/10.1038/nature08988
  * [33] Gross C, Zibold T, Nicklas E, Estève J and Oberthaler M K 2010 Nature 464 1165–1169 ISSN 1476-4687 URL http://dx.doi.org/10.1038/nature08919
  * [34] Berrada T, van Frank S, Bücker R, Schumm T, Schaff J F and Schmiedmayer J 2013 Nat. Commun. 4
  * [35] Pezze L, Smerzi A, Oberthaler M K, Schmied R and Treutlein P 2018 Reviews of Modern Physics 90 035005
  * [36] Blatt R and Roos C F 2012 Nat. Phys. 8 277–284 URL https://doi.org/10.1038/nphys2252
  * [37] Zhang J, Pagano G, Hess P W, Kyprianidis A, Becker P, Kaplan H, Gorshkov A V, Gong Z X and Monroe C 2017 Nature 551 601–604
  * [38] Monroe C, Campbell W C, Duan L M, Gong Z X, Gorshkov A V, Hess P W, Islam R, Kim K, Linke N M, Pagano G, Richerme P, Senko C and Yao N Y 2021 Rev. Mod. Phys. 93(2) 025001 URL https://link.aps.org/doi/10.1103/RevModPhys.93.025001
  * [39] Joshi M K, Kranzl F, Schuckert A, Lovas I, Maier C, Blatt R, Knap M and Roos C F 2022 Science 376 720–724
  * [40] Morong W, Liu F, Becker P, Collins K, Feng L, Kyprianidis A, Pagano G, You T, Gorshkov A and Monroe C 2021 Nature 599 393–398
  * [41] Dumitrescu P T, Bohnet J G, Gaebler J P, Hankin A, Hayes D, Kumar A, Neyenhuis B, Vasseur R and Potter A C 2022 Nature 607 463–467
  * [42] Feng L, Katz O, Haack C, Maghrebi M, Gorshkov A V, Gong Z, Cetina M and Monroe C 2023 Nature 623 713–717
  * [43] Jo G B, Shin Y, Will S, Pasquini T, Saba M, Ketterle W, Pritchard D, Vengalattore M and Prentiss M 2007 Phys. Rev. Lett. 98 030407
  * [44] Hofstetter W and Qin T 2018 Journal of Physics B: Atomic, Molecular and Optical Physics 51 082001 URL https://doi.org/10.1088%2F1361-6455%2Faaa31b
  * [45] Hernández Yanes T, Płodzień M, Mackoit Sinkevičienė M, Žlabys G, Juzeliūnas G and Witkowska E 2022 Phys. Rev. Lett. 129(9) 090403 URL https://link.aps.org/doi/10.1103/PhysRevLett.129.090403
  * [46] Schumm T, Hofferberth S, Andersson L M, Wildermuth S, Groth S, Bar-Joseph I, Schmiedmayer J and Krüger P 2005 Nat. Phys. 1 57–62
  * [47] Dziurawiec M, Hernández Yanes T, Płodzień M, Gajda M, Lewenstein M and Witkowska E 2023 Phys. Rev. A 107(1) 013311 URL https://link.aps.org/doi/10.1103/PhysRevA.107.013311
  * [48] Pawłowski K, Fadel M, Treutlein P, Castin Y and Sinatra A 2017 Physical Review A 95 ISSN 2469-9934 URL http://dx.doi.org/10.1103/PhysRevA.95.063609
  * [49] Gati R, Hemmerling B, Fölling J, Albiez M and Oberthaler M K 2006 Phys. Rev. Lett. 96 130404
  * [50] Dziarmaga J, Smerzi A, Zurek W and Bishop A 2002 Phys. Rev. Lett. 88 167001
  * [51] Trenkwalder A, Spagnolli G, Semeghini G, Coop S, Landini M, Castilho P, Pezze L, Modugno G, Inguscio M, Smerzi A et al. 2016 Nat. Phys. 12 826–829
  * [52] Invernizzi C, Korbman M, Campos Venuti L and Paris M G A 2008 Phys. Rev. A 78(4) 042106 URL https://link.aps.org/doi/10.1103/PhysRevA.78.042106
  * [53] Zanardi P, Paris M G A and Campos Venuti L 2008 Phys. Rev. A 78(4) 042105 URL https://link.aps.org/doi/10.1103/PhysRevA.78.042105
  * [54] Salvatori G, Mandarino A and Paris M G A 2014 Phys. Rev. A 90(2) 022111 URL https://link.aps.org/doi/10.1103/PhysRevA.90.022111
  * [55] Bina M, Amelio I and Paris M G A 2016 Phys. Rev. E 93(5) 052118 URL https://link.aps.org/doi/10.1103/PhysRevE.93.052118
  * [56] Garbe L, Bina M, Keller A, Paris M G A and Felicetti S 2020 Phys. Rev. Lett. 124(12) 120504 URL https://link.aps.org/doi/10.1103/PhysRevLett.124.120504
  * [57] Piga A, Aloy A, Lewenstein M and Frérot I 2019 Phys. Rev. Lett. 123(17) 170604 URL https://link.aps.org/doi/10.1103/PhysRevLett.123.170604
  * [58] Gühne O, Tóth G and Briegel H J 2005 New Journal of Physics 7 229 URL https://dx.doi.org/10.1088/1367-2630/7/1/229
  * [59] Tóth G and Apellaniz I 2014 Journal of Physics A: Mathematical and Theoretical 47 424006 URL https://dx.doi.org/10.1088/1751-8113/47/42/424006
  * [60] Iglói F and Tóth G 2023 Phys. Rev. Res. 5(1) 013158 URL https://link.aps.org/doi/10.1103/PhysRevResearch.5.013158
  * [61] Ziń P, Chwedeńczuk J, Oleś B, Sacha K and Trippenbach M 2008 Europhysics Letters 83 64007
  * [62] Chwedenczuk J, Hyllus P, Piazza F and Smerzi A 2012 New J. Phys. 14 093001
  * [63] Cavalcanti E G, Foster C J, Reid M D and Drummond P D 2007 Phys. Rev. Lett. 99 210405
  * [64] Cavalcanti E, He Q, Reid M and Wiseman H 2011 Phys. Rev. A 84 032115
  * [65] Niezgoda A, Panfil M and Chwedeńczuk J 2020 Phys. Rev. A 102(4) 042206 URL https://link.aps.org/doi/10.1103/PhysRevA.102.042206
  * [66] Chwedeńczuk J 2022 SciPost Physics Core 5 025
  * [67] Bell J S 1966 Rev. Mod. Phys. 38(3) 447–452
  * [68] Brunner N, Cavalcanti D, Pironio S, Scarani V and Wehner S 2014 Rev. Mod. Phys. 86(2) 419–478
  * [69] Greenberger D M, Horne M A and Zeilinger A 1989 69–72
  * [70] Płodzień M, Wasak T, Witkowska E, Lewenstein M and Chwedeńczuk J 2023 arXiv preprint arXiv:2306.06173
  * [71] Płodzień M, Lewenstein M, Witkowska E and Chwedeńczuk J 2022 Phys. Rev. Lett. 129(25) 250402 URL https://link.aps.org/doi/10.1103/PhysRevLett.129.250402
  * [72] Płodzień M, Wasak T, Witkowska E, Lewenstein M and Chwedeńczuk J 2023 Generation of scalable many-body bell correlations in spin chains with short-range two-body interactions (Preprint 2306.06173)