Building upon this relation, at each step $n\in\\{0,\ldots,N-1\\}$, we propose
to approximate $\nabla\log p^{\alpha}_{t_{0}}(Y^{(n)})$ in the recursion (134)
with a MCMC-based estimator of $u^{\alpha}_{t_{0}}(Y^{(n)})$ obtained by
sampling from $q^{\alpha}_{t_{0}}(\cdot|Y^{(n)})$. This results in the
Langevin-within-Langevin procedure summarized in Algorithm 2. We finally
highlight three main choices in the design of this algorithm: (i) we choose
the ULA initialization to be
$Y^{(0)}\sim\mathsf{N}(0,\sigma^{2}t_{0}\,\mathrm{I}_{d})$, relying on the
study conducted in Section D.3; (ii) consequently, we set the step-size
$\lambda$ to be slightly smaller than the variance of this distribution, i.e.,
$\lambda=\sigma^{2}t_{0}/2$; (iii) the MALA-based posterior sampling is
initialized with the best estimation of $u^{\alpha}_{t_{0}}(Y^{(n)})$ that we
have, i.e., $Y^{(n)}/\alpha(t_{0})$, recalling that $Y^{\alpha}_{t}/\alpha(t)$
and $u^{\alpha}_{t}(Y^{\alpha}_{t})$ have the same localization behaviour, see
Section B.2.

#### Algorithmic technicalities.

We now present three algorithmic subtleties featured in SLIPS.

  1. (a)

The step-size of MALA used to sample from the posterior is adapted. Using the
acceptance ratio of the Metropolis-Hasting filter, we geometrically decrease
(respectively increase) the step-size when the acceptance ratio is below
(respectively above) a target ratio of 75%;

  2. (b)

The last MCMC samples obtained when sampling the posterior at step $k$ will be
the first samples when sampling at step $k+1$. This persistent trick is
motivated by the fact that the posterior is expected to change little from one
iteration to the next. This behaviour is notably illustrated in Figure 3
(bottom row).

  3. (c)

Lastly, as mentioned and justified in Section B.2, we use the estimated
denoiser $\tilde{U}^{\alpha}_{T}$ at the final integration timestep $T$ as an
approximate sample from $\pi$ rather than $Y^{\alpha}_{T}/\alpha(T)$ to align
with related work.

###  F.2 Ablation study

In this section, we investigate the impact of the different hyper-parameters
of SLIPS. Similarly to Appendix D, we run our study by applying our algorithm
on the unbalanced bimodal Gaussian mixture defined in (77) with $d=10$. This
section will clarify three points:

  1. (a)

The behaviour of SLIPS within the assumptions of Theorem 3;

  2. (b)

The existence of a “sweet spot" for $t_{0}$ outside the assumptions of Theorem
3 and the impact of the estimation of $R_{\pi}$;

  3. (c)

The tuning of hyper-parameters in the Langevin-within-Langevin initialization
(see Algorithm 2).

#### Behaviour of SLIPS within the restrictive assumptions.

In this part only, we consider the case of a mixture of two Gaussian
distributions in dimension 5 given by
$\mathsf{N}(-(2/3)a\,\mathbf{1}_{5},\tau^{2}\mathrm{I}_{5})$ and
$\mathsf{N}(+(4/3)a\,\mathbf{1}_{5},\tau^{2}\mathrm{I}_{5})$, where
$\tau^{2}=0.1$ and $a>0$ is _not fixed_ , with respective weights $2/3$ and
$1/3$. By varying $a$, we vary $R$, recalling that $R=(4/3)a$, see Appendix C.
In Figure 9, we run SLIPS with different schemes while setting $t_{0}$ as the
mean value between $t_{p}$ and $t_{q}$ derived in Appendix E, in accordance
with Theorem 3 (_Left_), or derived by combining the results of Appendix E and
Appendix E, which is less restrictive (_Right_). When $p^{\alpha}_{t}$ and
$q^{\alpha}_{t}$ are both log-concave independently of $t$, i.e., when
$dR^{2}/\tau^{2}\leq 1$ for the general case (see Appendix E) or
$dR^{2}/\tau^{2}\leq 16/9$ in the refined Gaussian mixture setting (see
Appendix E), we always set $t_{0}=10^{-2}$. This choice ensures that the
Gaussian approximation which corresponds to the initialization of the
Langevin-within-Langevin algorithm is accurate (see Figure 7 for details).
This figure shows that choosing $t_{0}$ according to Appendix E leads to
accurate sampling below the threshold and degraded (or equal) performance
above this threshold.

#### Existence of a sweet spot for duality outside of the restrictive
assumptions.

Given the target distribution described in (77) with $d=10$, we have that
$R=4/3$ and $\tau^{2}=0.05$, which means that $R^{2}d/\tau^{2}\approx 350$.
Hence, the distribution at stake does not fit the additional assumption of
Theorem 3, or even the refined assumption of Appendix E. However, Figure 10
shows that there still exists a sweet spot for $t_{0}$ in this case. Moreover,
we observe in Figure 11 that a poor estimation of $R_{\pi}$ shifts the sweet
spot. This makes sense as the sweet spot likely corresponds to a SNR level.
Recalling that
$\log\operatorname{SNR}(t)=2\log(g(t))+\log(R_{\pi}^{2}/(\sigma^{2}d))$, if
$\sigma>R_{\pi}/\sqrt{d}$ (i.e., we overestimate $R_{\pi}$) the log-SNR is
shifted downwards and the resulting $t_{0}$ should be larger (as we see in the
right column of Figure 11); if $\sigma<R_{\pi}/\sqrt{d}$ (i.e., we
underestimate $R_{\pi}$), the resulting $t_{0}$ should be smaller (as we see
in the left column of Figure 11).

#### Choice of hyper-parameters in the Langevin-within-Langevin
initialization.

Figure 12 shows that only a few steps in the Langevin-within-Langevin
algorithm (see Algorithm 2) are needed to reach a stationary sampling error.

Figure 9: Sliced Wasserstein distance depending on the ratio
$R^{2}d/\tau^{2}$. Left: $R^{2}d/\tau^{2}=2$ is the threshold of duality of
log-concavity in Theorem 3 (_general case_). Right: $R^{2}d/\tau^{2}=32/9$ is
the threshold of duality of log-concavity in Appendix E (_refinement_). In
both cases, the grey line corresponds to the threshold where the log-concavity
conditions start to be restrictive on the marginal distribution of the
observation process and the corresponding posterior distribution. Figure 10:
Sliced Wasserstein distance and estimation error on the relative weights
depending on $t_{0}$ with different schemes - (Left) Standard (Middle)
Geom(1,1) (Right) Geom(2,1). Figure 11: Sliced Wasserstein distance depending
on $t_{0}$ for different values of $\sigma$. We denote
$\sigma_{r}=R_{\pi}/\sqrt{d}$. Figure 12: Sliced Wasserstein distance
depending on the number of Langevin-within-Langevin steps.

##  Appendix G Details on related work

In this appendix, we provide further details on related works, highlighting
similarities and identifying the main limitations.

###  G.1 Reverse Diffusion Monte Carlo (Huang et al., 2024)

Given $T>0$, the authors consider the denoising process $(Y_{t})_{t\in[0,T]}$
that is solution to the following SDE

$\displaystyle\mathrm{d}Y_{t}=\\{Y_{t}+2\nabla\log
p_{T-t}(Y_{t})\\}+\sqrt{2}\mathrm{d}B_{t}\;,\;Y_{0}\sim p_{T}\;,$ (136)

where $(B_{t})_{t\geq 0}$ is a Brownian motion in $\mathbb{R}^{d}$ and $p_{s}$
is the _intractable_ marginal distribution defined for any $s>0$ and any
$y\in\mathbb{R}^{d}$ by
$p_{s}(y)=\int_{\mathbb{R}^{d}}\mathsf{N}(y;\mathrm{e}^{-s}x,(1-\mathrm{e}^{-2s})\,\mathrm{I}_{d})\mathrm{d}\pi(x)$.
Under mild conditions on $\pi$ (Cattiaux et al., 2023), $(Y_{t})_{t\in[0,T]}$
is the time-reversed process of the standard Ornstein-Uhlenbeck process
defined in (37) and corresponds to a Variance-Preserving diffusion model (Song
et al., 2021). For any $t\in[0,T]$, given $X\sim\pi$, we therefore have

$\displaystyle Y_{t}=\mathrm{e}^{-(T-t)}X+\sqrt{1-\mathrm{e}^{-2(T-t)}}Z\;,$
(137)

where $Z$ is distributed according to the standard centered Gaussian
distribution. In this case, in a similar fashion to the SL posterior given in
(10), the conditional density of $X$ given $Y_{t}=y\in\mathbb{R}^{d}$ can be
defined as

$\displaystyle
q_{t}(x|y)\propto\pi(x)\mathsf{N}(x;\mathrm{e}^{T-t}y,(\mathrm{e}^{2(T-t)}-1)\,\mathrm{I}_{d}),\;t\in[0,T)\;,$
(138)

and the denoiser function as
$u_{t}(y)=\int_{\mathbb{R}^{d}}xq_{t}(x|y)\mathrm{d}x$. Similarly to our
framework, the random vector $u_{t}(Y_{t})$ can be seen as the denoiser of
$Y_{t}$. Following Tweedie’s formula given in Appendix A, the denoiser
function is related to the score of $p_{T-t}$ for any $t\in[0,T)$ by

$\displaystyle\nabla\log
p_{T-t}(y)=-\frac{y}{1-\mathrm{e}^{-2(T-t)}}+\frac{\mathrm{e}^{-(T-t)}}{1-\mathrm{e}^{-2(T-t)}}u_{t}(y)\;.$
(139)

Therefore, the SDE (136) describing the denoising process is strictly
equivalent to the following SDE

$\displaystyle\mathrm{d}Y_{t}=\left\\{\frac{\mathrm{e}^{-2(T-t)}+1}{\mathrm{e}^{-2(T-t)}-1}Y_{t}+\frac{2\mathrm{e}^{-(T-t)}}{1-\mathrm{e}^{-2(T-t)}}u_{t}(Y_{t})\right\\}\mathrm{d}t+\sqrt{2}\mathrm{d}B_{t},\;Y_{0}\sim
p_{T}\;.$ (140)

This last SDE has to be compared with the SDE describing our observation
process in (8): in both cases, the drift involves the denoiser function, which
is not tractable in practice. While $Y^{\alpha}_{T}/\alpha(T)$ is an
approximate sample from $\pi$ in our framework, here, it is given by $Y_{T}$.

The approach of (Huang et al., 2024) to sample from the SDE (140), while
handling the estimation of the denoiser $u_{t}$, is close to ours. Indeed,
after discretizing this SDE with a certain time discretization
$(t_{k})_{k=1}^{K}$ of $[0,T]$, they propose to estimate the denoiser at time
$t_{k}$ with a Markov Chain Monte Carlo estimation. To do so, at each step
$t_{k}$, they approximately sample from the random posterior
$q_{t_{k}}(\cdot|\tilde{Y}_{t_{k}})$ with ULA (while we use MALA), where
$\tilde{Y}_{t_{k}}$ is the $k$-th realization of their discretized process.
Besides this, they also propose a Langevin-within-Langevin initialization to
approximately sample from $p_{T}$, that involves the corresponding posterior
$q_{0}$. Although they briefly discuss the difficulty of this initialization
by considering the log-Sobolev constants of $p_{T}$ and $q_{0}$, their
analysis lacks readability and does not emphasize any crucial trade-off on $T$
(which corresponds to our $t_{0}$).

Their scheme can nonetheless be analyzed as a localization scheme under the
scope of what we call “duality of log-concavity", see Section 4.3. Assume that
$\pi$ verifies A and is not log-concave. Since we have
$q_{0}(x|y)\propto\pi(x)\mathsf{N}(x;\mathrm{e}^{T}y,(\mathrm{e}^{2T}-1)\,\mathrm{I}_{d})$,
one can show that

  1. (a)

if $T$ is large: $q_{0}$ will be close to $\pi$ (not log-concave), while
$p_{T}$ will be close to $\mathsf{N}(0,\mathrm{I}_{d})$ (log-concave),

  2. (b)

if $T$ is small: $q_{0}$ will localize as a Dirac mass (log-concave), while
$p_{T}$ will be close to $\pi$ (not log-concave).

Therefore, the setting of $T$ is crucial in RDMC to ensure a proper
initialization via Langevin-within-Langevin algorithm. However, the authors
claim that setting $T$ large would only incur computational waste and that
RDMC shows insensitivity to $T$, aside from the case where $T$ is too close to
0121212See (Huang et al., 2024, Appendix F.4).. As explained above, these
claims are not realistic, since sampling from the posterior $q_{0}$ becomes as
hard as sampling from $\pi$ when $T$ is large. We illustrate this fundamental
limitation in Figure 13.

Figure 13: Sampling the 8-Gaussians distribution via RDMC with different
values of $T$.

In contrast, we pay a particular attention to point out the importance of our
hyper-parameter $t_{0}$ throughout our analysis, with practical and
theoretical perspectives. Finally, we highlight that RDMC relies on a uniform
time discretization (while we propose a SNR-adapted discretization), and do
not use persistent initialization of their Markov chains in posterior
sampling.

###  G.2 Multi-Noise Measurements sampling (Saremi & Srivastava, 2022; Saremi
et al., 2023, 2024)

Given $M\geq 1$, the Multi-Noise Measurements (MNM) model introduced by
(Saremi & Srivastava, 2022) defines a sequence of random variables
$(Y^{m})_{m=1}^{M}$ as

$\displaystyle Y^{m}=X+\sigma Z^{m}\;,$ (141)

where $X\sim\pi$, $\sigma>0$ and $(Z^{m})_{m=1}^{M}$ are independently
distributed according to the standard centered Gaussian distribution. This
_non-Markovian_ stochastic process can be seen as a denoising process. The
rationale behind this formulation is that simulating an increasing number of
_equally_ noised measurements of a sample from $\pi$ helps to obtain more
information on this sample, and finally approximate it. We now explain how it
can be interpreted as a non-Markovian analog to the standard localization
scheme given in (1).

For any $m\in\\{1,\ldots,M\\}$, the conditional density of $X$ given the
$m$-tuple $Y^{1:m}=y_{1:m}\in(\mathbb{R}^{d})^{m}$ is defined as

$\displaystyle
q_{m}(x|y_{1:m})\propto\pi(x)\mathsf{N}(x;\bar{y}_{1:m},\sigma^{2}/m\,\mathrm{I}_{d})\;,$
(142)

where $\bar{y}_{1:m}$ denotes the empirical mean of the noised measurements,
i.e., $\bar{y}_{1:m}=(1/m)\sum_{i=1}^{m}y_{i}$. By defining
$u_{m}(y_{1:m})=\int_{\mathbb{R}^{d}}xq_{m}(x|y_{1:m})\mathrm{d}x$, one can
derive the Bayes estimator of $X$ given $Y^{1:m}$ as the random vector
$u_{m}(Y^{1:m})$ (Saremi & Srivastava, 2022). Interestingly, this denoiser
shares the same rate of convergence as the denoiser of standard stochastic
localization given in Section B.2 when $t=m$. Indeed, (Saremi et al., 2024,
Proposition 2) states that for any target distribution $\pi$, we have
$W_{2}(\pi,\tilde{\pi}_{m})\leq\sigma\sqrt{d/m}$, where $\tilde{\pi}_{m}$
denotes the distribution of $u_{m}(Y^{1:m})$. In other words, the MNM denoiser
localizes to $X$ with the same localization rate as the standard SL denoiser,
defined in Section 2, by taking $T=M$.

To sample from $\pi$, the MNM approach consists in first sampling $Y^{1:M}$,
and then computing the denoiser $u_{M}(Y^{1:M})$, in the same fashion as in
the SL framework. Recently, (Saremi et al., 2024) proposed to tackle the
sampling of the $M$-tuple $Y^{1:M}$ by first simulating $Y^{1}$ and then
sequentially sampling $Y^{m}$ given $Y^{1:m-1}$ for $m\in\\{2,\ldots,M\\}$
using a Monte Carlo Markov Chain method - in this case, Underdamped Langevin
Algorithm (Sachs et al., 2017). At step $m$, the Langevin procedure then
involves pointwise evaluations of the conditional score of the distribution of
$Y^{m}$ given $Y^{1:m-1}$ (or simply the score of the distribution of $Y_{1}$
when $m=1$). This introduces the _Once-At-a-Time_ (OAT) algorithm. The
interest of such method lies in the fact that, under Assumption A that may
include multi-modal distributions, the distributions to sample from are
increasingly log-concave with $m$, see (Saremi et al., 2024, Theorem 1).
Hence, the most challenging step of this approach lies in the sampling of
$Y^{1}$, in the same spirit as in SLIPS.

In their work, (Saremi et al., 2024) mainly consider the case where the scores
are analytically available. To handle realistic settings, they implement an IS
estimator, see (Saremi et al., 2024, Section 4.2.1), but their numerical
results show that its performance significantly degrades compared to setting
of perfect knowledge of the score, see (Saremi et al., 2024, Appendix H). We
include this approach in our numerical benchmark. Alternatively, for high-
dimensional settings, they propose an estimator based on posterior sampling,
see (Saremi et al., 2024, Appendix F.2.), in the same fashion as in SLIPS.
However, they do not pay attention to the limitation of this approach, in
particular at the initialization of their algorithm. Indeed, while (Saremi et
al., 2024, Theorem 1) suggests to take $\sigma$ very large to ensure that the
distribution of $Y^{1}$ is approximately log-concave (i.e., taking $t_{0}$
small in SLIPS), the corresponding posterior becomes unfortunately closer to
$\pi$, and then hard to sample with standard MCMC methods when $\pi$ is not
log-concave. This fundamental constraint is well illustrated in (Saremi et
al., 2024, Figure 7, left), where the posterior sampling approach is shown to
systematically fail for an arbitrary choice of $\sigma$, independently of the
computational budget allocated to the MCMC sampling. Therefore, the
combination of OAT with posterior sampling requires a trade-off on the hyper-
parameter $\sigma$, similarly to $t_{0}$ in SLIPS, that reflects once again
the “duality of log-concavity".

##  Appendix H Details on numerical experiments

###  H.1 Algorithms and hyper-parameters

#### Using the information from $R_{\pi}$ and adapting the algorithms.

As our definition of $\sigma$ in SLIPS requires to have an estimate of
$R_{\pi}$ (see Section 3.2), we informed other algorithms with such knowledge
when possible. More precisely, we use A and consider the approximation
$R_{\pi}^{2}=d(R^{2}+\tau^{2})$ where $R$ and $\tau$ are assumed to be known.

  * •

For SMC and AIS, we use a starting distribution $\rho_{0}$ as the centered
Gaussian with variance $R^{2}d\,\mathrm{I}_{d}$;

  * •

For OAT, we use $\sigma$ as per (Saremi et al., 2024, Equation 4.3) by setting
$\sigma^{2}=R^{2}d-\tau^{2}$ (ensuring log-concavity);

  * •

For RDMC, we do not include this information as the authors do not give any
heuristic on $T$ with respect to the variance of the target distribution.

Moreover, the implementations of the algorithms were slightly adjusted from
their general definition. SMC and AIS use a MALA kernel (which leaves the
target distribution invariant) as transition kernel. Additionally, RDMC also
uses a MALA kernel for posterior sampling instead of ULA. This reduces the
estimation bias and also enables automatic tuning of the Langevin step-size by
leveraging the acceptance ratio (here we adapt the step-size to maintain the
acceptance ratio at 75%). Still on RDMC, we drop the first 50% of the MCMC
samples to ignore the warm-up period in the estimation. We also ignore the
warm-up period with the same proportion in SLIPS. In SLIPS, we reuse the step-
size from the MCMC sampling of $q_{t_{k}}$ when sampling $q_{t_{k+1}}$ and we
also initialize the chains of the later with the last samples of the chain
from the former. This is an intuitive initialization when looking at the
bottom row of Figure 3 as the modes of the posterior seem to be stable and
$q_{t_{k}}$ is expected to be close to $q_{t_{k+1}}$.

#### Estimating the scores in OAT.

In OAT, denoting $p(y)$ the likelihood of the measurement process $Y=X+\sigma
Z$ with $X\sim\pi$ and $Z\sim\mathsf{N}(0,\mathrm{I}_{d})$, the score can be
written using the following identity from (Saremi et al., 2024, Appendix F) :
$\nabla\log p(y)=\sigma^{-1}\mathbb{E}_{Z\sim q(\cdot|y;\sigma)}[Z]$ where the
posterior $q(\cdot|y)$ is defined by $q(z|y)\propto\pi(y+\sigma
z)\mathsf{N}(z;0,\mathrm{I}_{d})$. This means that the score $\nabla\log p(y)$
can be estimated with the following IS estimator $\nabla\log
p(y)\approx\sigma^{-1}\sum_{i=1}^{N}w_{i}Z_{i}$ where
$Z_{i}\overset{i.i.d.}{\sim}\mathsf{N}(0,\mathrm{I}_{d})$ and
$w_{i}=\pi(y+\sigma Z_{i})/\sum_{j=1}^{N}\pi(y+\sigma Z_{j})$. We reduce the
variance of this estimator by applying the antithetic trick.

#### Selecting the hyper-parameters of the algorithms.

Table 3: Hyper-parameter grids used for SLIPS on different targets.
|  | Mixture of Gaussians and Bayesian |  $\phi^{4}$ model | Others
---|---|---|---|---
Standard | $\eta$ | $\\{5.0\\}$ | $\\{5.7,6.1\\}$ | $\\{5.0,5.7\\}$
$t_{0}$ | $\\{0.03,0.05,0.1,0.2,0.4\\}$ | $\\{0.8,1.0,1.2,1.4,1.8\\}$ | $\\{0.1,0.2,0.4,1.0,1.2\\}$
Geom(1,1) | $\eta$ | $\\{5.0\\}$ | $\\{5.7,6.1\\}$ | $\\{4.6,5.0\\}$
$t_{0}$ | $\\{0.03,0.05,0.1,0.15,0.25\\}$ | $\\{0.30,0.35,0.40,0.45\\}$ | $\\{0.1,0.15,0.20\\}$
Geom(2,1) | $\eta$ | $\\{5.0\\}$ | $\\{5.7,6.1\\}$ | $\\{4.6,5.0\\}$
$t_{0}$ | $\\{0.15,0.20,0.25,0.35,0.45\\}$ | $\\{0.40,0.45,0.50,0.55\\}$ | $\\{0.30,0.35,0.45\\}$

For each algorithm, we search its hyper-parameters within a predetermined
grid. The selection is based on the metrics which will be later detailed. The
metrics were computed by comparing $4096$ samples against true samples. We
globally fixed the computational budget by setting the SDE discretization of
SLIPS as $K=1024$. Moreover, we define the number of MCMC steps of SLIPS,
denoted by $L$, to be equal to $32$ except with the mixture of Gaussian in
high dimensions where it is equal to $48$, $64$ and $96$ in dimensions $32$,
$64$ and $128$ respectively and in the $\phi^{4}$ experiments were it is set
to $64$. The selected hyper-parameters for each algorithm are summarized in
Table 4. Below, we detail how the grids were built for each one.

  * •

The SMC and AIS algorithms define a sequence of annealed distributions
$\rho_{k}$ for $k\in\\{0,\ldots,K\\}$ from $\rho_{0}$ (defined above) to
$\rho_{K}=\pi$ as
$\rho_{k}\propto\exp((1-\beta_{k})\log\rho_{0}+\beta_{k}\log\pi))$ where
$(\beta_{k})_{k=0}^{K}$ is a linear schedule of size $K$ between $0$ and $1$.
Both algorithms used $N=4096$ particles;

  * •

The RDMC algorithm has a single hyper-parameter $T$ which we search within the
grid $T\in\\{-\log 0.99,-\log 0.95,-\log 0.9,-\log 0.8,-\log 0.7\\}$ given in
(Huang et al., 2024, Appendix F.1). However, we also find that large $T$ may
not work systematically (see Figure 13). We use $16$ steps of (Huang et al.,
2024, Algorithm 3) with $4$ MCMC chains. The chains are initialized with an IS
approximation of the posterior powered by $128$ particles. The SDE is
discretized over $K$ steps and the expectation leading to the drift is
estimated with $L$ Langevin steps. The initial sample is distributed according
to $\mathsf{N}(0,(1-\exp(-2T))\,\mathrm{I}_{d})$, the best estimation of
$p_{T}$ that we have, and the initial step-size is taken according to its
variance ;

  * •

The OAT algorithm has its parameter $M$ set as
$\left\lfloor(K/8)\times(L/32)\right\rfloor$ and uses $8$ MCMC steps per noise
level. Those specific choices ensure that the computational complexity of OAT
is on par with the other algorithms. The Langevin steps are done using the
underdamped Langevin algorithm as suggested by the authors. Its step-size is
searched in $\\{0.03,1.0\\}$ and the efficient friction is searched in
$\\{0.0625,0.05\\}$ as recommended by the authors. These prescriptions were
extracted from (Saremi et al., 2024, Appendix G). The OAT algorithm has
slightly shorter grid sizes because of its prohibitive computational cost;

  * •

The SLIPS algorithm is declined in three flavors depending on the choice of
the schedule $\alpha$. Since the choice of $t_{0}$ is sensitive to the
accuracy of the estimation of $R_{\pi}$ (different values of $R_{\pi}$ will
shift the log-SNR upwards or downwards), we decided to search the hyper-
parameters in different areas depending on the target distribution. Those
grids can be found in Table 3. The values for $t_{0}$ were chosen by
approximate equal log-SNR spacing in $[-3.5,-1.0]$ for Gaussian mixtures and
Bayesian logistic regression, in $[-1.0,0.2]$ for $\phi^{4}$ and $[-2.0,0.0]$
for the others. The values for $\eta$ were chosen to be around $5.0$.

Table 4: Hyper-parameters selected for the experiments for each algorithm and
target density.
Target | RDMC | OAT |  SLIPS Standard |  SLIPS Geom(1,1) |  SLIPS Geom(2,1)
---|---|---|---|---|---
|  Time $T$ |  step-size $\delta$ |  Frict. $\delta\gamma$ | $\eta$ | $t_{0}$ | $\eta$ | $t_{0}$ | $\eta$ | $t_{0}$
8-Gaussians | $-\log(0.80)$ | 0.05 | 0.03 | 5.7 | 0.60 | 5.7 | 0.35 | 5.0 | 0.35
Rings | $-\log(0.80)$ | 0.05 | 0.03 | 4.6 | 1.20 | 4.6 | 0.10 | 4.6 | 0.30
Funnel | $-\log(0.90)$ | 0.05 | 1.0 | 5.0 | 1.00 | 4.6 | 0.30 | 4.6 | 0.40
Mixture ($d=8$) | $-\log(0.70)$ | 0.0625 | 0.03 | 5.0 | 0.40 | 5.0 | 0.25 | 5.0 | 0.45
Mixture ($d=16$) | $-\log(0.70)$ | 0.05 | 0.03 | 5.0 | 0.20 | 5.0 | 0.15 | 5.0 | 0.35
Mixture ($d=32$) | $-\log(0.70)$ | 0.05 | 0.03 | 5.0 | 0.10 | 5.0 | 0.10 | 5.0 | 0.25
Mixture ($d=64$) | $-\log(0.70)$ | 0.0625 | 0.03 | 5.0 | 0.05 | 5.0 | 0.05 | 5.0 | 0.20
Ionosphere | $-\log(0.95)$ | 0.05 | 0.03 | 5.0 | 0.03 | 5.0 | 0.03 | 5.0 | 0.15
Sonar | $-\log(0.95)$ | 0.0625 | 0.03 | 5.0 | 0.03 | 5.0 | 0.03 | 5.0 | 0.15
$\phi^{4}$ ($b=0$) | $-\log(0.95)$ | 0.0625 | 0.03 | 5.7 | 0.80 | 5.7 | 0.30 | 5.7 | 0.40
$\phi^{4}$ ($b=0.025$) | $-\log(0.95)$ | 0.05 | 0.03 | 5.7 | 1.80 | 5.7 | 0.35 | 6.1 | 0.45
$\phi^{4}$ ($b=0.05$) | $-\log(0.70)$ | 0.05 | 0.03 | 6.1 | 1.00 | 5.7 | 0.30 | 5.7 | 0.40
$\phi^{4}$ ($b=0.075$) | $-\log(0.70)$ | 0.0625 | 0.03 | 5.7 | 1.80 | 5.7 | 0.35 | 5.7 | 0.40
$\phi^{4}$ ($b=0.1$) | $-\log(0.90)$ | 0.0625 | 0.03 | 5.7 | 1.40 | 5.7 | 0.45 | 5.7 | 0.40

###  H.2 Target distributions and metrics

#### 8 Gaussians, Rings and Funnel distributions and their respective metrics.

  1. (a)

The 8 Gaussians distribution consists of 8 equally weighted Gaussian
distributions with mean $\mathbf{m}_{i}=10\times(\cos(2\pi i/8),\sin(2\pi
i/8))$ for $i\in\\{0,\ldots,7\\}$ and covariance $0.7\,\mathrm{I}_{2}$. This
distribution satisfies A with $R=10/\sqrt{2}$ and $\tau^{2}=0.7$.

  2. (b)

The Rings distribution is the inverse polar reparameterization of a
distribution $p_{z}$ which has itself a decomposition into two univariate
marginals $p_{r}$ and $p_{\theta}$: $p_{r}$ is a mixture of 4 Gaussian
distributions $\mathsf{N}(i+1,0.15^{2})$ with $i\in\\{0,\ldots,3\\}$
describing the radial position and $p_{\theta}$ is a uniform distribution over
$[0,2\pi]$, which describes the angular position of the samples. This
distribution satisfies A with $R=4/\sqrt{2}$ and $\tau=0.15$.

  3. (c)

The Funnel distribution (Neal, 2003) has its density given by
$x_{1:10}\mapsto\mathsf{N}(x_{1};0,\sigma^{2}\mathrm{I}_{10})\mathsf{N}(x_{2:10};0,\exp(x_{1})\mathrm{I}_{10})$
where $\sigma^{2}=9$. We approximate $R$ by the maximum scalar standard
deviation of the distribution, i.e., $R=2.12$, and set $\tau=0$.

For 8 Gaussians and Rings, we evaluate the quality of sampling by
approximating the entropy-regularized 2-Wasserstein distance defined in
Appendix A, with regularization $\varepsilon=0.05$ via POT library (Flamary et
al., 2021). Regarding the Funnel, we evaluate the quality of sampling using
the Kolmogorov-Smirnov distance, see Appendix A, and adopt the sliced version
from (Grenioux et al., 2023b, Appendix D.1).

#### Bayesian Logistic Regression models.

Consider a training dataset $\mathcal{D}=\\{(x_{j},y_{j})\\}_{j=1}^{M}$ where
$x_{j}\in\mathbb{R}^{d}$ and $y_{j}\in\\{0,1\\}$ for all
$j\in\\{1,\ldots,M\\}$. We evaluate the likelihood of a pair $(x,y)$ as given
by $p(y|x;w,b)=\operatorname{Bernoulli}(y;\sigma(x^{T}x+b))$ where
$w\in\mathbb{R}^{d}$ is a weight vector, $b\in\mathbb{R}$ is an intercept and
$\sigma$ is the sigmoid function. Given a prior distribution $p(w,b)$, we
sample from the posterior distribution $p(w,b|\mathcal{D})\propto
p(\mathcal{D}|w,b)p(w,b)=\prod_{j=1}^{M}p(y_{j}|x_{j};w,b)p(w,b)$. The prior
is built as $p(w,b)=\mathsf{N}(w;0,\mathrm{I}_{d})\mathsf{N}(b;0,(2.5)^{2})$.
In our experiments, we approximate $R$ by the maximum scalar standard
deviation of the prior distribution, i.e., $R=2.5/\sqrt{d+1}$ and set
$\tau=0$. The quality of the samples are obtained by computing the mean
predictive log-likelihood (i.e., computing $p(w,b|\mathcal{D}_{\text{test}})$
with $\mathcal{D}_{\text{test}}$ a separate test dataset).

#### High dimension Gaussian mixture.

We consider the mixture of two Gaussian distributions defined in (77).
Following the computations of Appendix C, we set $R=4/3$ and $\tau^{2}=0.05$.
Complementary to the estimation error on the relative weights given in the
main paper, we also report the Sliced Wasserstein distance (Bonneel et al.,
2015; Nadjahi et al., 2019) on Figure 14 (Left). This figure shows that SLIPS
also recovers the local structure of the target distribution quite well.
Moreover, Figure 14 (Right) shows that SMC and AIS are unable to scale with
dimension and collapse into a single mode.

Figure 14: (Left) Sliced Wasserstein distance computed on bimodal Gaussian
mixtures with increasing dimension. (Right) Estimation error when computing
the relative weights of the corresponding two modes.

#### The $\phi^{4}$ distribution.

The $\phi^{4}$ model is a toy model defined as a continuous relaxation of the
Ising model that serves the study of phase transitions in statistical
mechanics. Following (Gabrié et al., 2022), we consider a version of the model
discretized on a 1-dimensional grid of size $d=100$. One configuration is
therefore a $d$-dimensional vector $(\phi_{i})_{i=1}^{d}$. We additionally
clip the field to 0 at both extremities by defining the extra
$\phi_{0}=\phi_{d+1}=0$. The negative log-density of the distribution writes

$\displaystyle\ln\pi_{h}(\phi)=-\beta\left(\frac{ad}{2}\sum_{i=1}^{d+1}(\phi_{i}-\phi_{i-1})^{2}+\frac{1}{4ad}\sum_{i=1}^{d}(1-\phi_{i}^{2})^{2}+h\phi_{i}\right)\;.$
(143)

We chose parameter values for which the system is bimodal, $a=0.1$ and inverse
temperature $\beta=20$, and vary the value of $h$. We denote by $w_{+}$ the
statistical occurrence of configurations such that $\phi_{d/2}>0$ and $w_{-}$
the statistical occurrence of configurations such that $\phi_{d/2}<0$. At
$h=0$, the measure is invariant under the symmetry $\phi\rightarrow-\phi$,
such that we expect $w_{+}=w_{-}$. For $h>0$, the negative mode dominates.

When $d$ is large, the relative probability of the modes can be estimated by a
Laplace approximations at 0-th and 2-nd order. Denoting by $\phi_{+}^{h}$ and
$\phi_{-}^{h}$ the local maxima of (143), these approximations yield
respectively,

$\displaystyle\frac{w_{-}}{w_{+}}\approx\frac{\pi_{h}(\phi_{-}^{h})}{\pi_{h}(\phi_{+}^{h})}\;,\quad\frac{w_{-}}{w_{+}}\approx\frac{\pi_{h}(\phi_{-}^{h})\times|\det
H_{h}(\phi_{-}^{h})|^{-1/2}}{\pi_{h}(\phi_{+}^{h})\times|\det
H_{h}(\phi_{+}^{h})|^{-1/2}}\;,$ (144)

where $H_{h}$ is the Hessian of the function
$\phi\rightarrow\ln\pi_{h}(\phi)$. In our experiments, we considered $R=4.5$
and $\tau=10^{-2}$ by running MALA chains started in each mode : $\tau^{2}$ is
set as the variance within the modes while $R$ corresponds to the distance
between the modes and $0_{100}$. In this setting, we observed that SMC and AIS
suffered from mode collapse, while OAT and RDMC produced degenerate samples;
in contrast, our methodology SLIPS recovers accurate samples from the target
distribution, with correct relative weights, as shown in Figure 5.


the other notations being obvious. One has also defined :

$\begin{array}[]{llll}\displaystyle
G(\delta_{P})=\left[\frac{\Delta_{A}}{2}+\epsilon\sin{\delta_{P}}-\epsilon^{\prime}\cos{\delta_{P}}\right]+\frac{\lambda_{3}}{2\sqrt{3}}(\sqrt{3}\lambda_{3}+\sqrt{2}\lambda_{0}+\lambda_{8})\;.\end{array}$
(136)

Three pieces are coming from the $VPPP$ :

$\begin{array}[]{ll}\displaystyle T_{VPPP}=&\displaystyle
C_{VPPP}\left[\sum_{V=\rho,\omega,\phi}\frac{F^{e}_{V\gamma}(s)}{D_{V}(s)}g_{V\pi}^{R}(s)\right]\epsilon_{\mu\nu\alpha\beta}~{}~{}\epsilon_{\mu}(\gamma)p_{0}^{\nu}p_{-}^{\alpha}p_{+}^{\beta}\,,\end{array}$
(137)

where the renormalized vector couplings $g_{V\pi}^{R}(s)$ to 3 pions have been
derived using the vector relation :

$g_{V\pi}^{R}(s)={\cal R}(s)g_{V\pi}^{I}~{},$ (138)

${\cal R}(s)$ being the matrix given in Eq. (38) or in Eq. (43) of the BHLS2
companion paper [ExtMod7] depending on whether the Primordial Breaking is
discarded or not. The components of the $g_{V\pi}^{I}$ vector which refers to
the coupling of the ideal vector field combinations are :

$\left\\{\begin{array}[]{ll}\displaystyle
g_{\rho\pi}^{I}=\frac{1}{4}\left[\frac{\Delta_{A}}{2}+\epsilon\sin{\delta_{P}}-\epsilon^{\prime}\cos{\delta_{P}}+\frac{\lambda_{3}}{2\sqrt{3}}(\sqrt{2}\lambda_{0}+\lambda_{8})\right]\,,\\\\[14.22636pt]
\displaystyle
g_{\omega\pi}^{I}=-\frac{3}{4}\left[1+\frac{2\xi_{0}+\xi_{8}}{3}-\frac{\lambda_{3}^{2}}{2}\right]\,,\\\\[14.22636pt]
\displaystyle
g_{\phi\pi}^{I}=-\frac{\sqrt{2}}{4}(\xi_{0}-\xi_{8})~{}~{}.\end{array}\right.$
(139)

The $V-\gamma$ amplitudes $F^{e}_{V\gamma}(s)$ and the inverse $\rho$
propagators have been constructed in Section 11 of [ExtMod7]. The inverse
propagators for the $\omega$ and $\phi$ mesons have been discussed and defined
in Section 9 of the same Reference. One has also defined :

$\displaystyle
C_{VPPP}=-\frac{3ige}{4\pi^{2}f_{\pi}^{3}}\left[c_{1}-c_{2}-c_{3}\right]\;.$
(140)

The $VVP$ Lagrangian piece in Equation (117) given in terms of ideal vector
fields has to be reexpressed in terms of their renormalized partners as
developped in Section 12 of [ExtMod7] – see Equations (70-75) therein. The
simplest way to write down $T(\gamma^{*}\rightarrow\pi^{0}\pi^{+}\pi^{-})$ in
a way easy to code within our global fit procedure is displayed just below.

One first defines the $H_{i}(s)$ functions :

$\left\\{\begin{array}[]{llll}\displaystyle H_{0}(s)=1\,,&\displaystyle
H_{1}(s)=\frac{1}{D_{\rho}(s_{+-})}+\frac{1}{D_{\rho}(s_{0+})}+\frac{1}{D_{\rho}(s_{0-})}\,,{}\\\\[14.22636pt]
\displaystyle
H_{2}(s)=\displaystyle\frac{1}{D_{\rho}(s_{+-})}~{}~{},&\displaystyle
H_{3}(s)=\displaystyle\widetilde{\alpha}(s_{+-})\left[\frac{1}{D_{\rho}(s_{+-})}-\frac{1}{D_{\omega}(s_{+-})}\right]\,,\end{array}\right.$
(141)

where $s$ is the incoming squared energy and the $s_{ij}$’s indicate the
invariant mass squared of the corresponding outgoing $(i,j)$ pairs; the tilde
mixing angles are those defined by Equation (43) in [ExtMod7]. $T_{VVP}$
depends on the 3 functions $(H_{i}(s),~{}i=1~{}\cdots 3)$ with the
$s$-dependent coefficients $F_{i}(s)$ given below.

Collecting all terms, the full amplitude writes :

$\displaystyle
T(\gamma^{*}\rightarrow\pi^{0}\pi^{+}\pi^{-})=\left[F_{0}(s)H_{0}(s)+C_{VVP}\sum_{i=1\cdots
3}F_{i}(s)H_{i}(s)\right]\epsilon_{\mu\nu\alpha\beta}~{}~{}\epsilon_{\mu}(\gamma)p_{0}^{\nu}p_{-}^{\alpha}p_{+}^{\beta}\,,$
(142)

with :

$\displaystyle
C_{VVP}=-i\frac{3egm^{2}}{8\pi^{2}f_{\pi}^{3}}(1+\Sigma_{V})c_{3}\;.$ (143)

In this way, to write down the full amplitude, the various $F_{i}(s)$
functions only depend on the incoming off-shell photon energy squared $s$; the
dependence upon the various sub-energies $s_{ij}$ is, instead, only carried by
the $H_{i}(s)$ functions as clear from Equations (141). One has :

$\left\\{\begin{array}[]{llll}\displaystyle
F_{0}(s)=C_{APPP}\left[1-G(\delta_{P})\right]+C_{VPPP}\left[\frac{F^{R}_{\rho\gamma}(s)}{D_{\rho}(s)}g_{\rho\pi}^{R}(s)+\frac{F^{R}_{\omega\gamma}(s)}{D_{\omega}(s)}g_{\omega\pi}^{R}(s)+\frac{F^{R}_{\phi\gamma}(s)}{D_{\phi}(s)}g_{\phi\pi}^{R}(s)\right]\,,\\\\[14.22636pt]
\displaystyle
F_{1}(s)=\displaystyle\widetilde{\alpha}(s)\frac{F^{R}_{\rho\gamma}(s)}{D_{\rho}(s)}+\left[1-\frac{\lambda_{3}^{2}}{2}\right]\left[1+\frac{2\xi_{0}+\xi_{8}}{3}\right]\frac{F^{R}_{\omega\gamma}(s)}{D_{\omega}(s)}+\left[\frac{\sqrt{2}}{3}(\xi_{0}-\xi_{8})+\widetilde{\gamma}(s)\right]\frac{F^{R}_{\phi\gamma}(s)}{D_{\phi}(s)}\,,\\\\[14.22636pt]
\displaystyle
F_{2}(s)=\displaystyle\left[\epsilon^{\prime}\cos{\delta_{P}}-\epsilon\sin{\delta_{P}}-\frac{\Delta_{A}}{2}-\frac{\lambda_{3}}{2\sqrt{3}}(\sqrt{2}\lambda_{0}+\lambda_{8})\right]\frac{F^{R}_{\rho\gamma}(s)}{D_{\rho}(s)}+2\xi_{3}\frac{F^{R}_{\omega\gamma}(s)}{D_{\omega}(s)}\,,\\\\[14.22636pt]
\displaystyle
F_{3}(s)=\displaystyle\frac{F^{R}_{\rho\gamma}(s)}{D_{\rho}(s)}\,,\end{array}\right.$
(144)

where $\widetilde{\alpha}(s)=\alpha(s)+\psi_{\omega}$ and
$\widetilde{\gamma}(s)=\gamma(s)+\psi_{0}$, to possibly keep track of the
Primordial Breaking [ExtMod7], $\alpha(s)$ and $\beta(s)$ being the angles
generated by the Dynamical Mixing mechanism [ExtMod3, ExtMod7].

A global fit to all cross-sections but
$e^{+}e^{-}\rightarrow\pi^{0}\pi^{+}\pi^{-}$ allows to yield the relevant
parameters with a good approximation; then, having at hand all ingredients
defining the $F_{i}(s)$’s, a first minimization run [ExtMod5] including the
$e^{+}e^{-}\rightarrow\pi^{0}\pi^{+}\pi^{-}$ cross-section can be performed to
also derive a first estimate for $c_{1}-c_{2}$. The output of this minuit
minimization run is then used as input for a next minimization step. This
initiates an iteration procedure involving all cross-sections which is carried
on up to convergence – when some criterion is met, e.g. generally,
$\Delta\chi^{2}\leq 0.3$ – for the global $\chi^{2}$ of all the processes
involved in the EBHLS2 procedure.

This method converges in a couple of minimization steps [ExtMod5]. What makes
such a minimization procedure unavoidable is that the
$e^{+}e^{-}\rightarrow\pi^{0}\pi^{+}\pi^{-}$ cross-section expression implies
to integrate over the $s_{ij}$’s at each step. This is obviously prohibitively
computer time consuming for a negligible gain. Hence, at each step, one starts
by tabulating coefficient functions, exhibited in the next expression between
curly brackets :

$\left.\begin{array}[]{ll}\displaystyle\sigma(e^{+}e^{-}\rightarrow\pi^{0}\pi^{+}\pi^{-},s)=&\displaystyle\frac{\alpha_{\rm
em}~{}s^{2}}{192\pi^{2}}\displaystyle\times\left[\left\\{\int
Gdxdy\right\\}|F_{0}(s)|^{2}\right.\\\\[8.5359pt]
&\left.+C_{VVP}^{2}\sum_{i,j=1\cdots 3}F_{i}(s)F_{j}^{*}(s)\left\\{\int
GH_{i}H_{j}^{*}dxdy\right\\}\right.\\\\[8.5359pt]
{}{}{}&\displaystyle\left.+C_{VVP}\sum_{i=1\cdots
3}\left(F_{0}(s)F_{i}^{*}(s)\left\\{\int
GH_{i}^{*}dxdy\right\\}\displaystyle+F_{0}^{*}(s)F_{i}(s)\left\\{\int
GH_{i}dxdy\right\\}\right)\right]\end{array}\right.$ (145)

and these tables are used all along the next step. Equation (145) uses the
Kuraev-Silagadze parametrization [Kuraev] and its kernel $G(x,y)$ function;
these are reminded in Appendix H of [ExtMod3].

##  Appendix D Energy Shift of Resonances and Secondary ISR Photons

For the BESIII data sample [BESIII_3pi] of concern here, photon radiation
effects have been unfolded by utilizing the PHOKHARA code [Czyz:2013xga],
which also includes second-order photon emission707070For a discussion of the
radiative correction tools see [Rodrigo:2001kf, Actis:2010gg].. In the leading
ISR process, $\sqrt{s}=3.773$ GeV is the energy at which the primary ISR
photon is emitted. Here we address the emission of a second photon near the
$\phi$ and $\omega$ resonances, where the corresponding radiation effect
appears resonance-enhanced.

The effect of a photon radiated off a Breit-Wigner (BW) resonance is well
known from the $Z$ resonance physics at LEP I [Barroso:1987ae, Bardin:1989qr,
Beenakker:1989ec, Jegerlehner:1991ed]. A specific analysis for the process
$e^{+}e^{-}\to\pi^{+}\pi^{-}$ may also be found in [Hoefer:2001mx,
Gluza:2002ui]. The leading effect is due to initial state radiation, where the
observed cross-section is given by a convolution717171Here one only considers
the ISR part at order $\cal{O}(\alpha)$ which yields the main shift between
the observed and the physical peak cross-sections. For a more complete
discussion of photon radiation effects see e.g. Eq. (39) ff in
[Hoefer:2001mx]. The final state radiation piece contributing to the
$\pi^{+}\pi^{-}\pi^{0}$ cross-section has been taken into account in the BES
III measurement. :

$\sigma^{\rm obs}(s)=\int_{0}^{1}dk\rho_{ini}(k)\;\sigma_{phys}(s(1-k))\>,$
(146)

where $k=E_{\gamma}/E_{e}$ is the energy emitted by the photon in units of the
electron or positron energy $E_{e}$ at which the resonance is formed and
$s=4E_{e}^{2}$. The photon spectral function is

$\rho_{ini}(k)=\beta
k^{\beta-1}\>(1+\delta_{1}^{v+s}+\dots)+\delta_{1}^{h}+\cdots\>,$ (147)

with $\beta=\frac{2\alpha}{\pi}\>(L-1)$, $L=\ln\frac{s}{m_{e}^{2}}$ and the
photon radiation corrections for the virtual + soft (v+s) and the hard (h)
parts read :

$\displaystyle\delta_{1}^{v+s}=\frac{\alpha}{\pi}\left(\frac{3}{2}\,L+\frac{\pi^{2}}{3}-2\right)\;;\delta_{1}^{h}=\frac{\alpha}{\pi}\>(1-L)(2-k)\>.$
(148)

Let us consider the narrow resonances $\omega$ and $\phi$ which are well
parametrized by the Breit-Wigner (BW) formula :

$\sigma_{\rm
BW}(s)=\frac{12\pi\Gamma_{e}\Gamma_{f}}{M_{R}^{2}}\frac{s}{(s-M_{R}^{2})^{2}+M_{R}^{2}\Gamma_{R}^{2}}\>,$
(149)

which has its maximum at $\sqrt{s}=M_{R}\>(1+\gamma^{2})^{\frac{1}{4}}$ with
$\gamma=\Gamma_{R}/M_{R}$. The peak value is
$\sigma_{max}=\frac{12\pi\Gamma_{e}\Gamma_{f}}{M_{R}^{2}\Gamma_{R}^{2}}\>(1+\frac{1}{4}\gamma^{2})$
and the half maximum locations read
$(\sqrt{s})_{\pm}=M_{R}\>(1+\frac{3}{8}\gamma^{2})\pm\frac{\Gamma_{R}}{2}\>(1-\frac{1}{8}\gamma^{2})$.
The leading photon radiation modifies these resonance parameters to

$\displaystyle\sqrt{s}_{\rm max}$ $\displaystyle=$ $\displaystyle
M_{R}+\frac{\beta\pi\Gamma_{R}}{8}-\frac{1}{4}\gamma^{2}\>M_{R}\>,$ (150)
$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\frac{\sigma^{\rm
obs}}{\sigma_{\rm
BW}}=\left(\frac{\Gamma_{R}}{M_{R}}\right)^{\beta}\>(1+\delta_{1}^{v+s})\>,$
(151) $\displaystyle\sqrt{s}_{+}-\sqrt{s}_{-}$ $\displaystyle=$
$\displaystyle\Gamma_{R}\>\left[1-\frac{\pi}{2}\beta\gamma-\frac{5}{8}\gamma^{2}+(\frac{\pi}{4}\beta+\frac{\beta}{2}\ln
2)(1+\beta+\frac{\pi}{4}\beta)\right]\>.$ (152)

The physical resonance appears wider with a reduced peak cross-section and the
position of the observed peak is shifted towards higher energies. Considering
the emission of a second photon, $s^{\prime}=s(1-k)$ in the convolution
formula (146) turns into $s^{\prime\prime}=s^{\prime}(1-k_{2})$ where
$s^{\prime}=s(1-k_{1})\sim M_{R}^{2}$, i.e. $s=M_{R}^{2}$ in the relations
listed before. This second photon shifts the peak energy of the $\omega$ and
$\phi$ by $\delta E$ as given in Table 14.

The last data column here displays, correspondingly, the values for the shifts
$\delta E^{\omega}_{BESIII}$ and $\delta E^{\phi}_{BESIII}$ returned by a
global EBHLS2 fit involving all data and, especially for the 3$\pi$ channel,
both the NSK and the BESIII samples; these fit values differ from the
corresponding calculated ones by only $1\sigma$.

Table 14:  The energy shifts in keV of the resonance locations by local ISR photon emission at $s=M_{R}^{2}$. We adopted PDG resonance parameter values (in MeV). The values for $(\delta E)^{\rm fit}$ are obtained in a global fit involving all 3$\pi$ data samples (BESIII+NSK).  | $M_{R}$ | $\Gamma_{R}$ | $\delta E$ | $(\delta E)^{\rm fit}$
---|---|---|---|---
$\omega$ | 782.65 | 8.49 | -403 | -486 $\pm$ 72
$\phi$ | 1019.46 | 4.26 | -213 | -135 $\pm$ 59

Interestingly the calculated leading order shifts quite well reproduce the
shifts found by optimizing the fits of the BESIII data (about a $1\sigma$
distance only for both signals). Therefore, the question is whether BESIII has
included (or not) the corresponding corrections when unfolding the raw data
from radiation effects.


#  Nucleon Energy Correlators for the Color Glass Condensate

Hao-Yu Liu  Xiaohui Liu<EMAIL_ADDRESS>Center of Advanced Quantum Studies,
Department of Physics, Beijing Normal University, Beijing, 100875, China  Ji-
Chen Pan  Institute of High Energy Physics and School of Physics, Chinese
Academy of Sciences, Beijing 100049, China  Feng Yuan<EMAIL_ADDRESS>Nuclear
Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA  Hua Xing Zhu<EMAIL_ADDRESS>Zhejiang Institute of Modern Physics,
Department of Physics, Zhejiang University, Hangzhou, 310027, China

###### Abstract

We demonstrate the recently proposed nucleon energy-energy correlator (nucleon
EEC) $f_{\rm EEC}(x,\theta)$ can unveil the gluon saturation in the small-$x$
regime in $eA$ collisions. The novelty of this probe is that it is fully
inclusive just like the deep-inelastic scattering (DIS), with no requirements
of jets or hadrons, but still provides an evident portal to the small-$x$
dynamics through the shape of the $\theta$-distribution. We find that the
saturation prediction is significantly different from the expectation of the
collinear factorization.

_Introduction._ Small-$x$ gluon saturation [1, 2, 3, 4, 5, 6] has been one of
the central focuses in nuclear physics community in recent years and will be a
major research area in the future Electron Ion Collider (EIC) [7, 8, 9]. An
effective field theory called color-glass-condensate (CGC) [4, 5, 6] has been
established to compute the hadronic and nuclear structure functions in deep
inelastic scattering (DIS) at small values of Bjorken-$x_{B}$ [10, 11]. The
CGC predicts the gluon saturation with a characteristic scale $Q_{s}$, as a
consequence of the small-$x$ nonlinear dynamics governed by the BK-JIMWLK
equation [12, 13, 14, 15, 16, 17]. The saturation scale $Q_{s}$ represents the
typical size of the gluon transverse momentum inside the nucleus and grows as
the momentum fraction $x\to 0$. For large nucleus and small-$x$, typically
$Q_{s}>\Lambda_{\rm QCD}$.

Previous experiments from DIS in $ep$ collisions at HERA and hadron
productions in $pA$ collisions at RHIC and LHC have shown some evidence of
gluon saturation at small-$x$. With the planned EIC in the horizon, this
physics will be explored in a systematic manner with unprecedented precision
[7, 8, 9]. Extensive studies have been carried out for the EIC experiments,
including the inclusive DIS structure functions at small-$x_{B}$ [18, 19, 20]
and the azimuthal correlations of di-jet/di-hadron/photon-jet/lepton-jet in
the inclusive or diffractive processes [21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52, 53, 54, 55, 56]. These processes are considered as promising
channels to look for the gluon saturation in $eA$ collisions.

In this manuscript, we present a novel approach to probe the gluon saturation
in $eA$ collisions in terms of the nucleon energy-energy correlator (nucleon
EEC) recently proposed in Ref. [57], which is an extension of the EEC [58, 59]
to the nucleon case. The EEC is the vacuum expectation of a set of final state
correlators to reformulate jet substructures [60, 61, 62, 63, 64, 65, 66, 67,
68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84], while the
nucleon EEC is the nucleon expectation of the initial-final state correlator.
The latter encodes the partonic angular distribution induced by the intrinsic
transverse momentum within the nucleon [57]. Therefore we expect the features
of the gluon saturation, especially the saturation scale $Q_{s}$ that measures
the size of the intrinsic transverse momentum, should be naturally imprinted
in the nucleon EEC. Our numeric results in Figs. 5 and 6 will show that the
saturation predictions have distinguished behaviors as compared to those from
the collinear factorization. From this comparison, we can further deduce the
saturation scales in $ep$ and $eA$ collisions, respectively.

The quark contribution to the nucleon EEC in the momentum space is defined as

$\displaystyle f_{q,\rm EEC}(x,\theta)=\\!\int\frac{dy^{-}}{4\pi
E_{A}}e^{-ixP\,y^{-}}\gamma^{+}\langle A|{\bar{\chi}}(y^{-})\,{\cal
E}(\theta)\chi(0)|A\rangle\,,$ (1)

where $x$ is the momentum fraction that initiates a scattering process,
meanwhile we measure the energy deposit in a detector at a given angle
$\theta$ from the initial state radiation and the remnants through the energy
operator ${\cal
E}(\theta)=\lim_{r\to\infty}\int_{0}^{\infty}dtT_{0\vec{n}}(t,\vec{n}r)r^{2}$
[90, 91, 92, 93], ${\cal E}(\theta)|X\rangle=\sum_{i\in
X}E_{i}\delta(\theta_{i}^{2}-\theta^{2})|X\rangle$ . The measured energy
deposit is normalized to the energy $E_{A}$ carried by the nucleus $A$. Here,
$\chi$ is the gauge invariant collinear quark field [85, 86, 87, 88, 89]. The
gluon EEC can be defined similarly. When $\theta E_{A}\sim\Lambda_{\rm QCD}$,
the $f_{\rm EEC}$ probes the intrinsic transverse dynamics of the nucleus $A$
through the operator ${\cal E}(\theta)$.

In the collinear factorization, it has been shown that when $\theta
E_{A}\gg\Lambda_{\rm QCD}$, the $f_{q,\rm EEC}(x,\theta)$ can be further
factorized as [57]

$\displaystyle f_{i,\rm
EEC}(x,\theta)=\int\frac{d\xi}{\xi}I_{ij}\left(\frac{x}{\xi},\theta\right)\left[\xi
f_{j/A}\left(\xi\right)\right]\,,$ (2)

where $f_{j/A}(\xi)$ is the collinear PDF, and $I_{ij}$ is the matching
coefficient found to be solely determined by the vacuum collinear splitting
functions [57].

Figure 1: The $f_{EEC}(x,\theta)$ in the collinear factorization (left) and
the CGC framework (right). Here $Q$ represents the center of mass energy of
the partonic cross section.

As the values of $x$ decreases, the $f_{q,\rm EEC}$ receives dramatically
enhanced contributions from the low $x$ gluon. In this regime, the non-linear
small-$x$ dynamics becomes important. Consequently, if compared to the
collinear factorization in which the distribution is determined by vacuum
collinear splitting, the shape of the $\theta$-distribution will be modified,
due to a sizable initial transverse momentum $q_{t}$ of order the saturation
scale $Q_{s}$, see, the illustrations in Fig. 1. Therefore, the nucleon EEC
can be used to probe the gluon saturation phenomenon and the small-$x$
dynamics, as we will show in the rest of this manuscript.

_The measurement and the factorization theorem._ We follow [57] to consider
the unpolarized DIS process $l+A\to l^{\prime}+X$ in the Breit frame. We
assume the nucleus is moving along the $+z$-direction. We measure the Bjorken
$x_{B}=\frac{-q^{2}}{2P\cdot q}$, the photon virtuality $Q^{2}=-q^{2}$ and the
energy $\sum_{i}E_{i}$ that deposits in a calorimeter at an angle $\theta$
with respect to the beam, as shown in Fig. 2. Here $q=l^{\prime}-l$ is the
momentum carried by the virtual photon. We then measure the weighted cross
section $\Sigma(Q^{2},x_{B},\theta)$ defined as

$\displaystyle\Sigma(Q^{2},x_{B},\theta)=\sum_{i}\int
d\sigma(x_{B},Q^{2},p_{i})\,\frac{E_{i}}{E_{A}}\,\delta(\theta^{2}-\theta_{i}^{2})\,,$
(3)

where $E_{A}$ is the energy carried by the incoming nucleus. We note that the
energy weight suppresses the soft contributions, which is an important feature
of the proposed measurement and its resulting nucleon EEC.

Figure 2: The $x_{B}$ and $Q^{2}$ measurement in DIS with a forward detector
that records the energy flow $\sum_{i}E_{i}$ at the angle $\theta$.

In order to probe the small-$x$ dynamics, we are particularly interested in
the scenario in which $x_{B}\ll 0.1$, and we place the detector in the far-
forward region such that $Q\theta\ll Q$ while $Q\theta\sim
Q_{s}\gg\Lambda_{\rm QCD}$. At this point, we emphasize that the measurement
involves neither additional hadron tagging nor jet clustering, and in contrast
to the TMD which restricts the events in the small $q_{t}$ region, this
approach is inclusive and does not veto events. It weights the full cross
section by the energy recorded at a certain angle $\theta$, therefore the
probe is as inclusive as the DIS but with additional control via $\theta$.

When $\theta Q\gg\Lambda_{\rm QCD}$, the weighted cross section can be
calculated perturbatively in the collinear factorization. More interestingly,
when $Q\theta\ll Q$, it has been shown that the $\Sigma(Q^{2},x_{B},\theta)$
fulfils the factorized form [57]

$\displaystyle\Sigma(Q^{2},x_{B},\theta)=\int\frac{dx}{x}{\hat{\sigma}}_{i,{\rm
DIS}}\left(\frac{x_{B}}{x},Q\right)f_{i,{\rm EEC}}(x,\theta)\,,$ (4)

where $\hat{\sigma}_{i,{\rm DIS}}$ is the fully inclusive partonic DIS cross
section. $f_{i,{\rm EEC}}$ is the nucleon EEC in Eq. (1). The
$\theta$-dependence enters entirely through the nucleon EEC $f_{\rm
EEC}(x,\theta)$, and therefore the $\theta$ distribution of the
$\Sigma(Q^{2},x_{B},\theta)$ probes the nucleon EEC when $\theta$ is small. We
note that $f_{\rm EEC}$ satisfies the same collinear evolution as the
collinear PDFs [57]

$\displaystyle\frac{df_{i,{\rm EEC}}(x,\theta)}{d\ln\mu}=P_{ij}\otimes
f_{j,{\rm EEC}}\,,$ (5)

as required by $d\Sigma/d\ln\mu=0$, and since $d{\hat{\sigma}}_{i,{\rm
DIS}}/d\ln\mu=-P_{ji}\otimes\hat{\sigma}_{j,{\rm DIS}}$. Here the convolution
in the momentum fraction is defined as $f\otimes
g(x)\equiv\int_{x}^{1}\frac{dz}{z}f\left(\frac{x}{z}\right)g(z)$. It is clear
from the evolution that there is no perturbative Sudakov suppression in
$f_{\rm EEC}$, due to the absence of the soft contribution in the collinear
factorization eliminated by the energy weight [57].

Figure 3:  The collinear splitting that initiates the DIS process and a
daughter parton that hits the detector at $\theta\ll 1$. The momentum
fractions are shown. We abbreviates $P^{+}$ with $P$ in this work for
simplicity notation.

The factorization theorem in Eq. (4) and Eq. (2) can be easily understood by
considering the leading contribution shown in Fig. 3, where a parton out of
the nucleus $A$ with momentum $\xi P$ splits into a parton with momentum
fraction $(1-z)\xi$ that hits the detector at $\theta$, and an internal line
with fraction $z\xi$ and virtuality $t=-\frac{\vec{k}_{t}^{2}}{1-z}$ that
initiates the partonic inclusive DIS process. Here $z$ is the momentum
fraction with respect to the incoming parton and
$k_{t}=\frac{1}{2}\xi(1-z)P\,\theta$ is the transverse momentum of the final
state parton. In the vacuum, the splitting is described by the leading order
vacuum collinear splitting kernel $\frac{1}{t}P_{ij}^{(0)}$. Since $\theta
Q\ll Q$, the chance for the radiations from the hard interaction to reach the
calorimeter vanishes as $\theta\to 0$. It is then found that in the small
$\theta$ limit,

$\displaystyle\Sigma(Q^{2},x_{B},\theta)=\int\frac{dx_{i}}{x_{i}}\hat{\sigma}_{i,{\rm
DIS}}\left(Q^{2},\frac{x_{B}}{x_{i}}\right)$ $\displaystyle\times\int d\xi
dz\frac{1}{\theta^{2}}\delta(x_{i}-\xi z)(1-z)\xi
P^{(0)}_{ij}\left(z\right)f_{j/A}(\xi)\,,\qquad$ (6)

which, aftet performing the $z$ integration, gives

$\displaystyle\Sigma(Q^{2},x_{B},\theta)=\int\frac{dx_{i}}{x_{i}}\hat{\sigma}_{i,{\rm
DIS}}\left(Q^{2},\frac{x_{B}}{x_{i}}\right)$
$\displaystyle\times\frac{1}{\theta^{2}}\int\frac{d\xi}{\xi}\left(1-\frac{x_{i}}{\xi}\right)P^{(0)}_{ij}\left(\frac{x_{i}}{\xi}\right)\left[\xi
f_{j/A}(\xi)\right].$ (7)

This produces the factorized form in Eq. (4) and Eq. (2) by identifying the
leading order matching coefficient
$I_{ij}^{(0)}(\xi,\theta)=\frac{1}{\theta^{2}}(1-\xi)P^{(0)}_{ij}(\xi)$.

If $x_{B}\ll 0.1$, the gluon density is overwhelmingly large and the leading
contribution to the $\Sigma(Q^{2},x_{B},\theta)$ is coming from

$\displaystyle\Sigma(Q^{2},x_{B},\theta)=\sum_{q}\frac{4\pi\alpha^{2}e_{q}^{2}}{Q^{4}}\,f_{q,{\rm
EEC}}(x_{B},\theta)\,,$ (8)

with

$\displaystyle f_{q,{\rm EEC}}(x,\theta)$ (9) $\displaystyle=$
$\displaystyle\frac{\alpha_{s}T_{R}}{2\pi\theta^{2}}\int_{x}^{1}\frac{d\xi}{\xi}(1-\xi)(\xi^{2}+(1-\xi)^{2})\,\left[\frac{x}{\xi}f_{g}\left(\frac{x}{\xi}\right)\right]\,.$

The collinear factorization predicts a $\frac{1}{\theta^{2}}$-scaling behavior
at ${\cal O}(\alpha_{s})$. For very small $\theta$, the scaling rule could
receive corrections from both the evolution of the $f_{\rm EEC}$ in Eq. (5)
and non-perturbative effects. But for generic small $\theta$, these effects
are mild and therefore $\theta^{2}\Sigma$ will be insensitive to the values of
$\theta$, up to ${\cal O}(\theta Q)$ power corrections. Furthermore, since the
energy weight kills the soft contribution, to all orders there will be no
perturbative Sudakov suppression in the small $\theta$ region in the collinear
factorization [57], as is clear from Eq. (5). Such a feature will be modified
by the small-$x$ dynamics as we will show.

Figure 4: The leading contribution to $f_{q,{\rm EEC}}(x,\theta)$ in the
small-$x$ region, where the double line represents the gauge link and the
gluon requires momentum $g^{+}=x_{g}P$ and $g_{t}\sim Q_{s}\sim\theta Q$.

_The nucleon EEC in the small- $x$ regime._ In the small-$x$ region, the gluon
density grows as $\frac{1}{x}$ and becomes overwhelmingly important and has to
be resummed to all orders. To realize such resummation in $f_{\rm EEC}$, we
invoke the CGC effective theory framework and follow the strategy in [94, 95,
96] to write the nucleon EEC in terms of the CGC dipole distribution 111The
complete calculation using the full dipole amplitude
$\psi_{T,L}^{\gamma^{\ast}\to q{\bar{q}}}$ for $\gamma^{\ast}\to q{\bar{q}}$
is presented in the Supplemental Material. Both approaches agree in the small
$\theta$ limit.. By evaluating the diagrams in Fig. 4, we find in the leading
logarithmic (LL) approximation

$\displaystyle f_{q,{\rm EEC}}(x_{B},\theta)$ $\displaystyle=$
$\displaystyle\frac{N_{C}S_{\perp}}{8\pi^{4}}\int d^{2}\vec{g}_{t}$ (10)
$\displaystyle\times\int_{\xi_{\rm cut}}^{1}\frac{d\xi}{\xi}{\cal
A}_{qg}\left(\xi,\theta,\vec{g}_{t}\right)\,F_{g,x_{B}}(\vec{g_{t}})\,,$

where $S_{\perp}$ is the averaged transverse area of the target nucleus and
$g_{t}\sim Q_{s}\sim\theta Q$ is the transverse momentum transfer.
$F_{g,x_{F}}=\int\frac{d^{2}\vec{r}}{4\pi^{2}}e^{-i\vec{g}_{t}\cdot\vec{r}_{t}}S_{x_{F}}^{(2)}(\vec{r}_{t})$
is the CGC dipole distribution evaluated at the scale $x_{F}$, where
$S^{(2)}_{x_{F}}(\vec{r}_{t})=\frac{1}{N_{C}}\langle{\rm
Tr}[W(\vec{r}_{t})W^{\dagger}(\vec{0})]\rangle_{x_{F}}$.
$x_{F}\frac{Q}{x_{B}}$ is the rapidity scale/boundary that separates the fast
moving modes being integrated out and the active slow moving partons in the
CGC effective frame work. In this work, we default to the natural choice
$x_{F}=x_{B}$. $\frac{1-\xi}{\xi}Q$ is the momentum “$+$”-component that
enters the detector. $\xi_{\rm cut}$ is determined by requiring the momentum
of the active quark does not exceed the rapidity boundary. Here the
coefficient ${\cal A}_{qg}$ is given by

$\displaystyle{\cal A}_{qg}(\xi,\theta,\vec{g}_{t})$ $\displaystyle=$
$\displaystyle\frac{1}{\theta^{2}}(1-\xi)\vec{k}_{t}^{2}(\vec{k}_{t}-\vec{g}_{t})^{2}$
(11)
$\displaystyle\times\,\left|\frac{\vec{k}_{t}}{\xi\vec{k}_{t}^{2}+(1-\xi)(\vec{k}_{t}-\vec{g}_{t})^{2}}-\frac{\vec{k}_{t}-\vec{g}_{t}}{(\vec{k}_{t}-\vec{g}_{t})^{2}}\right|^{2}\,,\quad$

with $k_{t}$ defined as $k_{t}=\frac{1-\xi}{\xi}\frac{Q}{2}\theta$, should be
of order $Q_{s}$.

It is easy to show that if $g_{t}\sim Q_{s}\ll Q\theta$, Eq. (10) reduces to
the $\frac{1}{\theta^{2}}$-scaling behavior of the collinear factorization in
Eq. (9). On the other hand, if $\theta Q\ll Q_{s}$, Eq. (10) scales as
$\theta^{0}$. We thus expect that in CGC, the $\theta^{2}\Sigma$ will be
independent of the $\theta$ for $\theta Q\gg Q_{s}$, however, contrary to the
collinear factorization, suppressed when $\theta Q\ll Q_{s}$. Meanwhile the
$\theta$ region between these two limits provide the opportunity to estimate
the saturation scale $Q_{s}$.

_Numerics._ Now we study the numerical impacts of the small-$x$ dynamics on
the shape of the $\theta^{2}\Sigma(Q^{2},x_{B},\theta)$ distribution from Eq.
(10), compared with the collinear prediction. We are particularly interested
in the region $\theta\ll 1$ where the $\theta$ distribution probes direcly the
$f_{\rm EEC}(x,\theta)$, see Eq. (4). For the small-$x$ dipole distribution
$S_{x_{F}}^{(2)}(\vec{r}_{t})$, we use both the MV model with rcBK running
[12, 97, 13, 98, 18, 99, 100, 101, 102, 103, 104] and the GBW model [105].

As for the MV model with rcBK running, we adopt the MV-like model [106] as the
initial condition, whose form is
$S^{(2)}_{x_{0}}(\vec{r}_{t})=\exp\left[-\frac{(r_{t}^{2}Q_{s0}^{2})^{\gamma}}{4}\ln\left(\frac{1}{\Lambda
r_{t}}+e\right)\right]$, where we choose $x_{0}=0.01$,
$\gamma=1.119,\Lambda=0.241~{}\text{GeV},Q_{s0}^{2}=A^{1/3}0.168~{}\text{GeV}^{2}$
with $A$ the atomic number. We use the solution to the LL BK evolution with
$\alpha_{s}$ running [98, 101, 106] of the dipole distribution to evolve the
dipole distribution from $x_{0}$ to $x_{F}$. In our calculation, we use the
result fitted from the HERA data for the transverse area of the nucleus
$S_{\perp}$ [19]. The GBW model is implemented using
$S^{(2)}_{x_{F}}(\vec{r}_{t})=\exp\left[-\frac{1}{4}r_{t}^{2}Q_{s}^{2}(x_{F})\right]$,
where $Q_{s}^{2}(x_{F})=A_{N}(x_{0}/x_{F})^{\lambda}\>{\rm GeV}^{2}$ and we
use $x_{0}=2.24\times 10^{-4}$, $\lambda=0.27$ and $A_{N}=1$ for the proton
while $A_{N}=5$ for the Au [107].

Figure 5:  Normalized $\theta^{2}\Sigma(Q^{2},x_{B},\theta)$ distribution for
proton and Au, with $x_{B}=0.003$, for both $Q^{2}=25\,{\rm
GeV}^{2}\,,\sqrt{s}=105\,{\rm GeV}$ (left panel) and $Q^{2}=100\,{\rm
GeV}^{2}\,,\sqrt{s}=318\,{\rm GeV}$ (right panel). Green dots are from
Pythia82 simulation, in which we demand $0.00295<x_{B}<0.00305$, and $25{\rm
GeV}^{2}\leq Q^{2}<35{\rm GeV}^{2}$ for the left panel while $100{\rm
GeV}^{2}\leq Q^{2}<110{\rm GeV}^{2}$ for the right.

In Fig. 5, we show the CGC predictions for
$\theta^{2}\Sigma(Q^{2},x_{B},\theta)$ as a function of $\theta$. Since we are
only interested in the shape, we normalized the distribution by
$\int_{\theta_{\rm min}}^{\theta_{\rm max}}d\theta\theta^{2}\Sigma$. We fixed
$x_{B}=3\times 10^{-3}$ and choose $Q^{2}=25\,{\rm GeV}^{2}$,
$\sqrt{s}=105\,\rm{GeV}$ (left panel) and $Q^{2}=100\,{\rm GeV}^{2}$,
$\sqrt{s}=318\,\rm{GeV}$ (right panel). We present predictions from CGC for
both proton (in purple lines) and Au (in orange), by the MV model with rcBK
running and GBW model. We see both models predict similar shapes in the
$\theta$ spectrum, in which the small-$\theta$ region is suppressed. They are
impressively different from the collinear expectations (in red lines and green
dots). In the figure, the collinear predictions (red lines) are made out of
the complete fixed order $\alpha_{s}$ calculation, without the $Q\theta\ll Q$
approximation, using CT18A [108] and EPPS21 [109] PDF sets for proton and Au,
respectively. To validate our collinear calculation and to estimate the size
of the evolution effect in Eq. (5), we also run a Pythia82 simulation [110]
for the proton case, where the LL resummation is performed. We see that for
large $\theta$, the fixed order calculations agree well with the the Pythia
simulation, while for small $\theta$ values, the resummation effects could be
sizable but does not suppress the small-$\theta$ region due to the absence of
the perturbative Sudakov factor in $f_{\rm EEC}$ in the collinear
factorization. The collinear prediction for the Au follows closely the
proton’s. The notable difference demonstrates that the $f_{\rm EEC}(x,\theta)$
could serve as a clean probe of the small-$x$ phenomenon. For comparison, we
also show the predictions from the full CGC calculation derived in the
Supplemental Material using the GBW model in purple circles.

In Fig. 5, the proton spectrum turns into a plateau for large values of
$\theta$, which is expected from Eq. (10) when $Q\theta\gg Q_{s}$. We can
define a turning point around which the slope of the distribution starts to
switch its monotonicity. The turning point allows us to estimate the size of
the saturation scale $Q_{s}$. For instance, from the left panel of Fig. 5, the
turning point for the proton is roughly around $\theta\sim 0.15-0.2$ and thus
$Q_{s}\sim\theta Q\sim 0.75-1.0\,{\rm GeV}$ which is consistent with the
values of $\Lambda_{\rm QCD}$. And we estimate the saturation scale for the Au
will be around $\theta\sim 0.4-0.5$ and thus $Q_{s}\sim\theta Q\sim
2-2.5\,{\rm GeV}$. The right panel of Fig. 5 is similar to the left, but with
$Q^{2}=100\,{\rm GeV}^{2}$. Since $Q^{2}$ is larger, the distribution enters
the plateau earlier as expected. Now the turning point for the Au is around
$\theta\sim 0.2-0.3$ which again indicates that $Q_{s}\sim 2-3\,{\rm GeV}$,
consistent with the $Q^{2}=25\,{\rm GeV}^{2}$ case.

We can further introduce the nuclear modification factor
$R_{pA}=\frac{A^{-1}\Sigma_{A}(Q^{2},x_{B},\theta)}{\Sigma_{p}(Q^{2},x_{B},\theta)}$,
which helps to reduce the systematics. In the collinear factorization, for
$\theta Q\gg\Lambda_{\rm QCD}$, the $\theta$ distribution is determined by the
matching coefficient $I_{ij}$ as predicted by Eq. (8), which is independent of
the incoming nucleus species. Thus taking the ratio $R_{Ap}$ reduces the
impacts from perturbative higher order corrections as well as possible non-
perturbative hadronization effects, and the collinear factorization predicts
the $R_{Ap}$ insensitive to the $\theta$ values, as showed explicitly as red
lines in Fig. 6.

Figure 6: $R_{pA}$ as a function of $\theta$, with $x_{B}=3\times 10^{-3}$
using the MV model with rcBK running and collinear factorization.

.

Once again, the small-$x$ formalism changes the pattern as we observed in Fig.
6, where the modification factor $R_{pA}$ is suppressed in the small $\theta$
region, while converges toward around unity as $\theta$ becomes large and
$Q\theta\gg Q_{s}$.

_Conclusions._ In this manuscript, we have proposed the nucleon energy-energy
correlator (nucleon EEC) as a new probe of the gluon saturation phenomenon in
DIS at the future electron-ion colliders. In particular, we have shown that
the $\theta$-shape of the nucleon EEC $f_{\rm EEC}(x,\theta)$ behaves
differently in the collinear factorization theorem and the CGC formalism. The
drastic difference is due to the intrinsic transverse momentum of order
$Q_{s}$ induced by the non-linear small-$x$ dynamics. We thus expect the
$f_{\rm EEC}$ to complement the other standard small-$x$ processes and offer a
great opportunity to pin down the onset of the gluon saturation phenomenon in
$eA$-collisions.

The advantage of the nucleon EEC probe, as compared to other standard
small-$x$ processes, is that it is fully inclusive, involving no fragmentation
functions and jet clustering. Therefore the observable is expected to be clean
both theoretically and experimentally. Extensions to other observables that
are induced by the intrinsic transverse dynamics of the nucleon/nucleus shall
follow. For polarized hadron beam, we can study the spin asymmetry by adding
the azimuthal dependence to the energy operator ${\cal E}$ and we expect
different asymmetries to be predicted in the collinear and CGC frameworks. We
hope that the results presented in this manuscript will motivate carrying out
the proposed measurement at the current and future electron-ion facilitieis,
meanwhile stimulate further applications of the nucleon EEC in nuclear
structure studies.

###### Acknowledgements.

_Acknowledgement._ We are grateful to Farid Salazar, Hongxi Xing, Jian Zhou
for useful discussions. This work is supported by the Natural Science
Foundation of China under contract No. 12175016 (X. L.), No. 11975200 (H. X.
Z.), and the Office of Science of the U.S. Department of Energy under Contract
No. DE-AC02-05CH11231 (F.Y.).

## References

  * [1] L. V. Gribov, E. M. Levin, and M. G. Ryskin,  Phys. Rept. 100, 1 (1983).
  * [2] A. H. Mueller and J.-w. Qiu,  Nucl. Phys. B 268, 427 (1986).
  * [3] A. H. Mueller,  Nucl. Phys. B 335, 115 (1990).
  * [4] L. D. McLerran and R. Venugopalan,  Phys. Rev. D 49, 2233 (1994), hep-ph/9309289.
  * [5] L. D. McLerran and R. Venugopalan,  Phys. Rev. D 49, 3352 (1994), hep-ph/9311205.
  * [6] L. D. McLerran and R. Venugopalan,  Phys. Rev. D 50, 2225 (1994), hep-ph/9402335.
  * [7] A. Accardi et al.,  Eur. Phys. J. A 52, 268 (2016), 1212.1701.
  * [8] R. Abdul Khalek et al.,  (2021), 2103.05419.
  * [9] Proceedings, Probing Nucleons and Nuclei in High Energy Collisions: Dedicated to the Physics of the Electron Ion Collider: Seattle (WA), United States, October 1 - November 16, 2018, WSP, 2020, 2002.12333.
  * [10] F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan,  Ann. Rev. Nucl. Part. Sci. 60, 463 (2010), 1002.0333.
  * [11] E. Iancu and R. Venugopalan,  The Color glass condensate and high-energy scattering in QCD (, 2003), pp. 249–3363, hep-ph/0303204.
  * [12] I. Balitsky,  Nucl. Phys. B 463, 99 (1996), hep-ph/9509348.
  * [13] Y. V. Kovchegov,  Phys. Rev. D 60, 034008 (1999), hep-ph/9901281.
  * [14] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert,  Nucl. Phys. B 504, 415 (1997), hep-ph/9701284.
  * [15] J. Jalilian-Marian, A. Kovner, A. Leonidov, and H. Weigert,  Phys. Rev. D 59, 014014 (1998), hep-ph/9706377.
  * [16] E. Iancu, A. Leonidov, and L. D. McLerran,  Nucl. Phys. A 692, 583 (2001), hep-ph/0011241.
  * [17] E. Ferreiro, E. Iancu, A. Leonidov, and L. McLerran,  Nucl. Phys. A 703, 489 (2002), hep-ph/0109115.
  * [18] J. L. Albacete, N. Armesto, J. G. Milhano, P. Quiroga-Arias, and C. A. Salgado,  Eur. Phys. J. C 71, 1705 (2011), 1012.4408.
  * [19] T. Lappi and H. Mäntysaari,  Phys. Rev. D 88, 114020 (2013), 1309.6963.
  * [20] B. Ducloué, H. Hänninen, T. Lappi, and Y. Zhu,  Phys. Rev. D 96, 094017 (2017), 1708.07328.
  * [21] F. Dominguez, B.-W. Xiao, and F. Yuan,  Phys. Rev. Lett. 106, 022301 (2011), 1009.2141.
  * [22] F. Dominguez, C. Marquet, B.-W. Xiao, and F. Yuan,  Phys. Rev. D 83, 105005 (2011), 1101.0715.
  * [23] A. H. Mueller, B.-W. Xiao, and F. Yuan,  Phys. Rev. D 88, 114010 (2013), 1308.2993.
  * [24] A. Metz and J. Zhou,  Phys. Rev. D 84, 051503 (2011), 1105.1991.
  * [25] F. Dominguez, J.-W. Qiu, B.-W. Xiao, and F. Yuan,  Phys. Rev. D 85, 045003 (2012), 1109.6293.
  * [26] A. Dumitru, T. Lappi, and V. Skokov,  Phys. Rev. Lett. 115, 252301 (2015), 1508.04438.
  * [27] A. Dumitru and V. Skokov,  Phys. Rev. D 94, 014030 (2016), 1605.02739.
  * [28] D. Boer, P. J. Mulders, C. Pisano, and J. Zhou,  JHEP 08, 001 (2016), 1605.07934.
  * [29] G. Beuf,  Phys. Rev. D 96, 074033 (2017), 1708.06557.
  * [30] A. Dumitru, V. Skokov, and T. Ullrich,  Phys. Rev. C 99, 015204 (2019), 1809.02615.
  * [31] H. Mäntysaari, N. Mueller, F. Salazar, and B. Schenke,  Phys. Rev. Lett. 124, 112301 (2020), 1912.05586.
  * [32] Y.-Y. Zhao, M.-M. Xu, L.-Z. Chen, D.-H. Zhang, and Y.-F. Wu,  Phys. Rev. D 104, 114032 (2021), 2105.08818.
  * [33] R. Boussarie, H. Mäntysaari, F. Salazar, and B. Schenke,  JHEP 09, 178 (2021), 2106.11301.
  * [34] P. Caucal, F. Salazar, and R. Venugopalan,  JHEP 11, 222 (2021), 2108.06347.
  * [35] Y.-Y. Zhang and X.-N. Wang,  Phys. Rev. D 105, 034015 (2022), 2104.04520.
  * [36] P. Taels, T. Altinoluk, G. Beuf, and C. Marquet,  JHEP 10, 184 (2022), 2204.11650.
  * [37] P. Caucal, F. Salazar, B. Schenke, and R. Venugopalan,  JHEP 11, 169 (2022), 2208.13872.
  * [38] R. Boussarie, A. V. Grabovsky, L. Szymanowski, and S. Wallon,  JHEP 09, 026 (2014), 1405.7676.
  * [39] R. Boussarie, A. V. Grabovsky, L. Szymanowski, and S. Wallon,  JHEP 11, 149 (2016), 1606.00419.
  * [40] F. Salazar and B. Schenke,  Phys. Rev. D 100, 034007 (2019), 1905.03763.
  * [41] R. Boussarie, A. V. Grabovsky, L. Szymanowski, and S. Wallon,  Phys. Rev. D 100, 074020 (2019), 1905.07371.
  * [42] D. Boer and C. Setyadi,  Phys. Rev. D 104, 074006 (2021), 2106.15148.
  * [43] E. Iancu, A. H. Mueller, and D. N. Triantafyllopoulos,  Phys. Rev. Lett. 128, 202001 (2022), 2112.06353.
  * [44] E. Iancu, A. H. Mueller, D. N. Triantafyllopoulos, and S. Y. Wei,  JHEP 10, 103 (2022), 2207.06268.
  * [45] Y. Hatta, B.-W. Xiao, and F. Yuan,  Phys. Rev. Lett. 116, 202301 (2016), 1601.01585.
  * [46] T. Altinoluk, N. Armesto, G. Beuf, and A. H. Rezaeian,  Phys. Lett. B 758, 373 (2016), 1511.07452.
  * [47] H. Mäntysaari, N. Mueller, and B. Schenke,  Phys. Rev. D 99, 074004 (2019), 1902.05087.
  * [48] Y. Hagiwara, C. Zhang, J. Zhou, and Y.-j. Zhou,  Phys. Rev. D 104, 094021 (2021), 2106.13466.
  * [49] L. Zheng, E. C. Aschenauer, J. H. Lee, and B.-W. Xiao,  Phys. Rev. D 89, 074037 (2014), 1403.2413.
  * [50] F. Bergabo and J. Jalilian-Marian,  Nucl. Phys. A 1018, 122358 (2022), 2108.10428.
  * [51] F. Bergabo and J. Jalilian-Marian,  Phys. Rev. D 106, 054035 (2022), 2207.03606.
  * [52] E. Iancu and Y. Mulian,  (2022), 2211.04837.
  * [53] M. Fucilla, A. V. Grabovsky, E. Li, L. Szymanowski, and S. Wallon,  (2022), 2211.05774.
  * [54] I. Kolbé, K. Roy, F. Salazar, B. Schenke, and R. Venugopalan,  JHEP 01, 052 (2021), 2008.04372.
  * [55] X.-B. Tong, B.-W. Xiao, and Y.-Y. Zhang,  (2022), 2211.01647.
  * [56] T. Altinoluk, G. Beuf, A. Czajka, and A. Tymowska,  (2022), 2212.10484.
  * [57] X. Liu and H. X. Zhu,  (2022), 2209.02080.
  * [58] C. Basham, L. S. Brown, S. D. Ellis, and S. T. Love,  Phys. Rev. Lett. 41, 1585 (1978).
  * [59] C. Basham, L. Brown, S. Ellis, and S. Love,  Phys. Rev. D 19, 2018 (1979).
  * [60] H. Chen, I. Moult, X. Zhang, and H. X. Zhu,  Phys. Rev. D 102, 054012 (2020), 2004.11381.
  * [61] D. M. Hofman and J. Maldacena,  JHEP 05, 012 (2008), 0803.1467.
  * [62] A. Belitsky, S. Hohenegger, G. Korchemsky, E. Sokatchev, and A. Zhiboedov,  Phys. Rev. Lett. 112, 071601 (2014), 1311.6800.
  * [63] A. Belitsky, S. Hohenegger, G. Korchemsky, E. Sokatchev, and A. Zhiboedov,  Nucl. Phys. B 884, 305 (2014), 1309.0769.
  * [64] M. Kologlu, P. Kravchuk, D. Simmons-Duffin, and A. Zhiboedov,  JHEP 01, 128 (2021), 1905.01311.
  * [65] G. Korchemsky,  JHEP 01, 008 (2020), 1905.01444.
  * [66] L. J. Dixon, I. Moult, and H. X. Zhu,  Phys. Rev. D 100, 014009 (2019), 1905.01310.
  * [67] H. Chen et al.,  JHEP 08, 028 (2020), 1912.11050.
  * [68] H. Chen, I. Moult, and H. X. Zhu,  Phys. Rev. Lett. 126, 112003 (2021), 2011.02492.
  * [69] C.-H. Chang, M. Kologlu, P. Kravchuk, D. Simmons-Duffin, and A. Zhiboedov,  JHEP 05, 059 (2022), 2010.04726.
  * [70] Y. Li, I. Moult, S. S. van Velzen, W. J. Waalewijn, and H. X. Zhu,  Phys. Rev. Lett. 128, 182001 (2022), 2108.01674.
  * [71] M. Jaarsma, Y. Li, I. Moult, W. Waalewijn, and H. X. Zhu,  JHEP 06, 139 (2022), 2201.05166.
  * [72] P. T. Komiske, I. Moult, J. Thaler, and H. X. Zhu,  (2022), 2201.07800.
  * [73] J. Holguin, I. Moult, A. Pathak, and M. Procura,  (2022), 2201.08393.
  * [74] K. Yan and X. Zhang,  Phys. Rev. Lett. 129, 021602 (2022), 2203.04349.
  * [75] H. Chen, I. Moult, J. Sandor, and H. X. Zhu,  (2022), 2202.04085.
  * [76] C.-H. Chang and D. Simmons-Duffin,  (2022), 2202.04090.
  * [77] H. Chen, I. Moult, J. Thaler, and H. X. Zhu,  JHEP 07, 146 (2022), 2205.02857.
  * [78] K. Lee, B. Meçaj, and I. Moult,  (2022), 2205.03414.
  * [79] A. J. Larkoski,  (2022), 2205.12375.
  * [80] L. Ricci and M. Riembau,  (2022), 2207.03511.
  * [81] T.-Z. Yang and X. Zhang,  (2022), 2208.01051.
  * [82] C. Andres et al.,  (2022), 2209.11236.
  * [83] H. Chen et al.,  (2022), 2210.10058.
  * [84] E. Craft, K. Lee, B. Meçaj, and I. Moult,  (2022), 2210.09311.
  * [85] C. W. Bauer, S. Fleming, D. Pirjol, and I. W. Stewart,  Phys. Rev. D 63, 114020 (2001), hep-ph/0011336.
  * [86] C. W. Bauer, D. Pirjol, and I. W. Stewart,  Phys. Rev. D 65, 054022 (2002), hep-ph/0109045.
  * [87] C. W. Bauer and I. W. Stewart,  Phys. Lett. B 516, 134 (2001), hep-ph/0107001.
  * [88] M. Beneke, A. P. Chapovsky, M. Diehl, and T. Feldmann,  Nucl. Phys. B 643, 431 (2002), hep-ph/0206152.
  * [89] C. W. Bauer, S. Fleming, D. Pirjol, I. Z. Rothstein, and I. W. Stewart,  Phys. Rev. D 66, 014017 (2002), hep-ph/0202088.
  * [90] N. A. Sveshnikov and F. V. Tkachov,  Phys. Lett. B 382, 403 (1996), hep-ph/9512370.
  * [91] F. V. Tkachov,  Int. J. Mod. Phys. A 12, 5411 (1997), hep-ph/9601308.
  * [92] G. P. Korchemsky and G. F. Sterman,  Nucl. Phys. B 555, 335 (1999), hep-ph/9902341.
  * [93] C. W. Bauer, S. P. Fleming, C. Lee, and G. F. Sterman,  Phys. Rev. D 78, 034027 (2008), 0801.4569.
  * [94] C. Marquet, B.-W. Xiao, and F. Yuan,  Phys. Lett. B 682, 207 (2009), 0906.1454.
  * [95] B.-W. Xiao, F. Yuan, and J. Zhou,  Nucl. Phys. B 921, 104 (2017), 1703.06163.
  * [96] J. Zhou,  Phys. Rev. D 99, 054026 (2019), 1807.00506.
  * [97] Y. V. Kovchegov and H. Weigert,  Nucl. Phys. A 789, 260 (2007), hep-ph/0612071.
  * [98] Y. V. Kovchegov and H. Weigert,  Nucl. Phys. A 784, 188 (2007), hep-ph/0609090.
  * [99] K. J. Golec-Biernat, L. Motyka, and A. M. Stasto,  Phys. Rev. D 65, 074037 (2002), hep-ph/0110325.
  * [100] J. L. Albacete and Y. V. Kovchegov,  Phys. Rev. D 75, 125021 (2007), 0704.0612.
  * [101] I. Balitsky,  Phys. Rev. D 75, 014001 (2007), hep-ph/0609105.
  * [102] E. Gardi, J. Kuokkanen, K. Rummukainen, and H. Weigert,  Nucl. Phys. A 784, 282 (2007), hep-ph/0609087.
  * [103] I. Balitsky and G. A. Chirilli,  Phys. Rev. D 77, 014019 (2008), 0710.4330.
  * [104] J. Berger and A. Stasto,  Phys. Rev. D 83, 034015 (2011), 1010.0671.
  * [105] K. J. Golec-Biernat and M. Wusthoff,  Phys. Rev. D 59, 014017 (1998), hep-ph/9807513.
  * [106] H. Fujii and K. Watanabe,  Nucl. Phys. A 915, 1 (2013), 1304.2221.
  * [107] K. Golec-Biernat and S. Sapeta,  JHEP 03, 102 (2018), 1711.11360.
  * [108] T.-J. Hou et al.,  Phys. Rev. D 103, 014013 (2021), 1912.10053.
  * [109] K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado,  Eur. Phys. J. C 82, 413 (2022), 2112.12462.
  * [110] T. Sjöstrand et al.,  Comput. Phys. Commun. 191, 159 (2015), 1410.3012.

See pages 1 of sm.pdf See pages 2 of sm.pdf