# Harmonics of Parton Saturation in Lepton-Jet Correlations at the EIC

Xuan-Bo Tong<EMAIL_ADDRESS>School of Science and Engineering, The
Chinese University of Hong Kong, Shenzhen, Shenzhen, Guangdong, 518172, P.R.
China  University of Science and Technology of China, Hefei, Anhui, 230026,
P.R.China  Bo-Wen Xiao<EMAIL_ADDRESS>School of Science and
Engineering, The Chinese University of Hong Kong, Shenzhen, Shenzhen,
Guangdong, 518172, P.R. China  Yuan-Yuan Zhang<EMAIL_ADDRESS>School of Science and Engineering, The Chinese University of Hong Kong,
Shenzhen, Shenzhen, Guangdong, 518172, P.R. China  University of Science and
Technology of China, Hefei, Anhui, 230026, P.R.China

###### Abstract

Parton saturation is one of the most intriguing phenomena in the high energy
nuclear physics research frontier, especially in the upcoming era of the
Electron-Ion Collider (EIC). The lepton-jet correlation in deep inelastic
scattering provides us with a new gateway to the parton saturation at the EIC.
In particular, we demonstrate that azimuthal angle anisotropies of the lepton-
jet correlation are sensitive to the strength of the saturation momentum in
the EIC kinematic region. In contrast to the predictions based on the
collinear framework calculation, significant nuclear modification of the
anisotropies is observed when we compare the saturation physics results in
$e+p$ and $e+\text{Au}$ scatterings. By measuring these harmonic coefficients
at the EIC, one can conduct quantitative analysis in different collisional
systems and unveil compelling evidence for saturation effects.

Introduction. The parton saturation phenomenon Gribov:1983ivg ; Mueller:1985wy
; Mueller:1989st ; McLerran:1993ni ; McLerran:1993ka ; McLerran:1994vd ,
predicted by the color glass condensate (CGC) theory Gelis:2010nm ;
Iancu:2003xm , is one of the cutting-edge topics that many physicists strive
to explore, and it is also one of the profound questions that the electron-ion
collider (EIC) Boer:2011fh ; Accardi:2012qut ; AbdulKhalek:2021gbh ;
AbdulKhalek:2022hcn  intends to address. Various observables have been
proposed to search for the signals of the parton saturation at the future EIC
(see e.g., a recent review Morreale:2021pnn ). Especially, two-particle
correlations Dominguez:2010xd ; Dominguez:2011wm ; Mueller:2013wwa ;
Metz:2011wb ; Dominguez:2011br ; Dumitru:2015gaa ; Dumitru:2016jku ;
Boer:2016fqd ; Dumitru:2018kuw ; Mantysaari:2019hkq ; Zhao:2021kae ;
Boussarie:2021ybe ; Caucal:2021ent ; Zhang:2021tcc ; Taels:2022tza ;
Caucal:2022ulg ; Boussarie:2014lxa ; Boussarie:2016ogo ; Salazar:2019ncp ;
Boussarie:2019ero ; Boer:2021upt ; Iancu:2021rup ; Iancu:2022lcw ;
Hatta:2016dxp ; Altinoluk:2015dpi ; Mantysaari:2019csc ; Hagiwara:2021xkf ;
Zheng:2014vka ; Bergabo:2021woe ; Bergabo:2022tcu ; Iancu:2022gpw ;
Fucilla:2022wcg ; Kolbe:2020tlq , including the azimuthal correlations of di-
jet Dominguez:2010xd ; Dominguez:2011wm ; Mueller:2013wwa ; Metz:2011wb ;
Dominguez:2011br ; Dumitru:2015gaa ; Dumitru:2016jku ; Boer:2016fqd ;
Dumitru:2018kuw ; Mantysaari:2019hkq ; Zhao:2021kae ; Boussarie:2021ybe ;
Caucal:2021ent ; Zhang:2021tcc ; Taels:2022tza ; Caucal:2022ulg ;
Boussarie:2014lxa ; Boussarie:2016ogo ; Salazar:2019ncp ; Boussarie:2019ero ;
Boer:2021upt ; Iancu:2021rup ; Iancu:2022lcw ; Hatta:2016dxp ;
Altinoluk:2015dpi ; Mantysaari:2019csc  or di-hadron Zheng:2014vka ;
Bergabo:2021woe ; Bergabo:2022tcu ; Iancu:2022gpw ; Fucilla:2022wcg  in the
inclusive or diffractive processes, are one of the promising observables under
intensive study in the last decade, see also the tremendous applications in
the LHC and RHIC Jalilian-Marian:2004vhw ; Kharzeev:2004bw ; Marquet:2007vb ;
Tuchin:2009nf ; Dumitru:2010ak ; Dominguez:2010xd ; Dominguez:2011wm ;
Mueller:2013wwa ; Kutak:2012rf ; vanHameren:2014ala ; Kotko:2015ura ;
vanHameren:2016ftb ; vanHameren:2019ysa ; vanHameren:2020rqt ; Marquet:2016cgx
; Klein:2019qfb ; Iancu:2020mos ; Bolognino:2021mrc ; Al-Mashad:2022zbq ;
Hagiwara:2017fye ; Kotko:2017oxg ; Bhattacharya:2018lgm ; Boussarie:2018zwg ;
Albacete:2010pg ; Stasto:2011ru ; Lappi:2012nh ; Iancu:2013dta ;
Albacete:2018ruq ; Stasto:2018rci ; Jalilian-Marian:2005qbq ; Jalilian-
Marian:2012wwi ; Stasto:2012ru ; Rezaeian:2012wa ; Basso:2015pba ;
Rezaeian:2016szi ; Basso:2016ulb ; Boer:2017xpy ; Benic:2017znu ;
Goncalves:2020tvh ; Benic:2022ixp ; Boer:2022njw ; Gelis:2002fw ;
Kovner:2014qea ; Kovner:2015rna ; Marquet:2019ltn ; Akcakaya:2012si ;
Marquet:2017xwy . The azimuthal angle anisotropies of two-particle
correlations can be obtained through harmonic analysis, and they offer unique
perspectives and channels to identify and probe the parton saturation at the
EIC Metz:2011wb ; Dominguez:2011br ; Dumitru:2015gaa ; Dumitru:2016jku ;
Dumitru:2018kuw ; Boer:2016fqd ; Mantysaari:2019hkq ; Zhao:2021kae ;
Boussarie:2021ybe ; Hatta:2016dxp ; Altinoluk:2015dpi ; Mantysaari:2019csc ;
Hagiwara:2021xkf .

Recently, the lepton-jet correlation (LJC), a novel type of two-particle
correlation, was introduced in Refs. Liu:2020dct ; Liu:2018trl  to study the
quark transverse-momentum-dependent (TMD) distributions in deep inelastic
scatterings. It is interesting to note that this observable is directly
defined in the lab frame (or the center-of-mass frame of the incoming lepton
and hadron) and thus can be conveniently measured by the detector. This
process allows one to use the scattered lepton as the tagging reference, which
can be identified and measured rather precisely, and the recoiled jet as the
probe Arratia:2019vju  to study the soft dynamics of the process in the back-
to-back correlation region, where the system’s imbalance transverse momentum
$|\vec{q}_{\perp}|$ is much softer than its relative transverse momentum
$|\vec{P}_{\perp}|$. Since the Bremsstrahlung emission of gluon (and photon)
favors the collinear region, it was shown that large radiation-induced
anisotropy in LJC Hatta:2021jcd  can be generated within the TMD and collinear
factorization framework Hatta:2021jcd ; Hatta:2020bgy ; Catani:2014qha ;
Catani:2017tuc ; Ju:2022wia .

The objective of this paper is to use the LJC anisotropy as a novel probe to
study parton saturation in detail. As an early exploration, we present
quantitative evaluations of several LJC anisotropies with various inputs in
the CGC and collinear formalisms. In particular, we demonstrate that the
$\langle\cos n\phi\rangle$ coefficients, namely, the harmonics of LJC, can
shed light on distinct behaviors of the saturation phenomenon in
small-$q_{\perp}$ region. Here $\phi$ is the azimuthal angle difference
between $\vec{q}_{\perp}$ and $\vec{P}_{\perp}$, and $n$ is the harmonic
number.

First, due to the collinear enhancement of final-state jet gluon emission, the
near-cone gluon emission has the largest probability in general, and thus
produces significant initial anisotropies. Second, through numerical and
analytic analysis, we find that the saturation effects can significantly
suppress the anisotropy in the small $q_{\perp}$ region especially in
$e+\textrm{Au}$ collisions. In the CGC framework, the typical transverse
momentum of incoming quarks is related to and determined by the saturation
momentum $Q_{s}$, and it is due to random multiple scatterings and small-$x$
gluon emissions, which have no angular preference. Thus, one expects that
saturation effects tend to reduce the final measured LJC anisotropies. For the
cases with appreciable saturation momentum $Q_{s}$, the anisotropy can be
largely washed out in the region $|\vec{q}_{\perp}|\lesssim Q_{s}$. In
addition, due to the nuclear enhancement of the saturation momentum
$Q_{s}^{2}\propto A^{1/3}$ with $A$ the nucleon number, a strong nuclear
modification of the final anisotropies is then predicted by changing the
target hadron from the proton to heavy nuclei, such as Au. As a result, this
new observable at the upcoming EIC may yield additional compelling evidence
for the parton saturation.

Lepton-jet correlation at small $x$. Let us briefly describe the theoretical
framework as follows. We investigate the production of the lepton and jet in
deep inelastic scatterings as follows

$\displaystyle\ell(k)+\textrm{A}(P)\rightarrow\ell^{\prime}\left(k_{\ell}\right)+\operatorname{Jet}\left(k\textbf{}_{J}\right)+X,$
(1)

and measure their azimuthal angular correlations from their momenta defined
above. We are interested in the correlation limit, where
$|\vec{q}_{\perp}|=|\vec{k}_{\ell\perp}+\vec{k}_{J\perp}|\ll|\vec{P}_{\perp}|=|(\vec{k}_{\ell\perp}-\vec{k}_{J\perp})/2|\sim|\vec{k}_{\ell\perp}|\sim|\vec{k}_{J\perp}|$.
Here the center-of-mass energy is $s_{eN}=(k+P)^{2}$ and the momentum transfer
is $Q^{2}=(k-k_{\ell})^{2}$. Based on the TMD factorization and the small-$x$
formalism Marquet:2009ca ; Dominguez:2010xd ; Xiao:2010sa ; Xiao:2010sp ;
Dominguez:2011wm , we derive the azimuthal angle dependent cross section for
LJC in the correlation limit. To include the parton showers effects, the
Sudakov logarithms are resumed to all orders using the technique developed in
Refs. Mueller:2012uf ; Mueller:2013wwa . In terms of the harmonic expansion,
the cross section can be cast into

$\displaystyle\frac{d^{5}\sigma(\ell
P\rightarrow\ell^{\prime}J)}{dy_{\ell}d^{2}P_{\perp}d^{2}q_{\perp}}=\sigma_{0}\int\frac{b_{\perp}db_{\perp}}{2\pi}W(x,b_{\perp})\Big{[}J_{0}(q_{\perp}b_{\perp})$
$\displaystyle\qquad+\sum_{n=1}^{\infty}2\cos(n\phi)\alpha_{s}(\mu_{b})\frac{C_{F}c_{n}(R)}{n\pi}J_{n}(q_{\perp}b_{\perp})\Big{]}~{}.$
(2)

At leading order (LO), the correlation is produced by a quark (or antiquark)
jet and the away side lepton. In Eq.(2), $x$ is the longitudinal momentum
fraction of the incoming quark w.r.t. the target proton or nucleus. From
kinematics, one finds
$x=k_{\ell\perp}/\sqrt{s_{eN}}\left(e^{-y_{\ell}}+e^{-y_{J}}\right)$ and
$1=k_{\ell\perp}/\sqrt{s_{eN}}\left(e^{y_{\ell}}+e^{y_{J}}\right)$, where
$y_{\ell},\,y_{J}$ are the rapidities of the produced lepton and jet,
respectively. $\sigma_{0}$ stands for the LO hard cross section, and $W$ is
defined as

$W(x,b_{\perp})=\sum_{q}e_{q}^{2}xf_{q}(x,b_{\perp})e^{-\text{Sud}(b_{\perp})}~{},$
(3)

where $f_{q}(x,b_{\perp})$ denotes the unintegrated quark distribution in the
coordinate space. In small-$x$ limit, it is related to the dipole scattering
amplitude McLerran:1998nk ; Venugopalan:1999wu ; Mueller:1999wm ;
Marquet:2009ca ; Xiao:2017yya  as follows

$\displaystyle xf_{q}($ $\displaystyle
x,b_{\perp})=\frac{N_{c}S_{\perp}}{8\pi^{4}}\int
dyd^{2}r_{\perp}\frac{(\vec{b}_{\perp}+\vec{r}_{\perp})\cdot\vec{r}_{\perp}}{|\vec{b}_{\perp}+\vec{r}_{\perp}||\vec{r}_{\perp}|}$
$\displaystyle\times
yK_{1}(\sqrt{y}|\vec{b}_{\perp}+\vec{r}_{\perp}|)K_{1}\left(\sqrt{y}|\vec{r}_{\perp}|\right)$
$\displaystyle\times\Big{[}1+{\cal S}_{x}(b_{\perp})-{\cal
S}_{x}(b_{\perp}+r_{\perp})-{\cal S}_{x}(r_{\perp})\Big{]}~{},$ (4)

where $S_{\perp}$ is the averaged transverse area of the target hadron and
${\cal S}_{x}(r_{\perp})$ represents the dipole scattering matrix with
$r_{\perp}$ the dipole transverse size.

The Sudakov form factor $\text{Sud}(b_{\perp})$ in Eq.(3) can be separated
into the perturbative and non-perturbative parts,
$\text{Sud}(b_{\perp})=\text{Sud}_{\text{P}}(b_{\perp})+\text{Sud}_{\text{NP}}(b_{\perp})$.
The perturbative Sudakov factor reads

$\text{Sud}_{\text{P}}=\int_{\mu_{b}}^{Q}\frac{d\mu}{\mu}\frac{\alpha_{s}(\mu)C_{F}}{\pi}\Big{[}\ln\frac{Q^{2}}{\mu^{2}}+\ln\frac{Q^{2}}{P_{\perp}^{2}}+c_{0}(R)\Big{]}~{},$
(5)

where $\mu_{b}=b_{0}/b_{\perp}^{*}$ with $b_{0}\equiv 2e^{-\gamma_{E}}$ and
$\gamma_{E}$ the Euler constant;
$b^{*}_{\perp}=b_{\perp}/\sqrt{1+b_{\perp}^{2}/b_{\text{max}}^{2}}$ with
$b_{\text{max}}=1.5~{}$GeV-1. Since the small-$x$ quark distributions have
already contained the non-perturbative information on the target,
$\text{Sud}_{\text{NP}}$ is usually not included in the CGC formalism for
simplicity. For comparison purposes, we also compute the LJC in the TMD
framework Hatta:2021jcd  with the collinear quark distribution and the
standard non-perturbative $\text{Sud}_{\text{NP}}$ Sun:2014dqm ;
Prokudin:2015ysa .

The Fourier coefficient $c_{n}(R)$ is extracted by analyzing the soft gluon
emission from the initial quark and final jet, following the procedure in Ref.
Hatta:2021jcd . After subtracting soft gluon emissions inside jets, one
obtains the residual dependence on the jet cone size $R$ in $c_{n}(R)$. The
exact expression of $c_{n}(R)$ can be found in Ref. Hatta:2021jcd . In our
calculation, we choose $R=0.4$ , where
$c_{0}=1.79,~{}c_{1}=2.60,~{}c_{2}=0.95,~{}c_{3}=0.25$.

With the harmonic expansion of the differential cross section in Eq.(2), the
azimuthal anisotropy of the LJC reads

$\langle\cos n\phi\rangle=\frac{\int
b_{\perp}db_{\perp}J_{n}\left(q_{\perp}b_{\perp}\right)W(x,b_{\perp})\alpha_{s}(\mu_{b})\frac{C_{F}c_{n}(R)}{n\pi}}{\int
b_{\perp}db_{\perp}J_{0}\left(q_{\perp}b_{\perp}\right)W(x,b_{\perp})}~{}.$
(6)

Using various models for the dipole S-matrix, we can predict the LJC azimuthal
anisotropy $\langle\cos n\phi\rangle$ at small $x$, and show the dependence on
the momentum imbalance $q_{\perp}$ and the sensitivity to the saturation
momentum $Q_{s}$.

At last, let us discuss the QED modification to the anisotropy $\langle\cos
n\phi\rangle$. As noted by Ref. Hatta:2021jcd , this modification can occur
when we turn on the final-state photon emissions from the electron side. In
this case, we find that the formula of $\langle\cos n\phi\rangle$ in Eq.(6)
can be modified by the following two corrections. First, the QED Sudakov form
factor should be added, which reads

$\displaystyle\text{Sud}^{\gamma}=\int^{Q}_{b_{0}/b_{\perp}}\frac{d\mu}{\mu}\frac{\alpha_{e}}{\pi}\Big{[}\ln\frac{Q^{2}}{\mu^{2}}+\ln\frac{Q^{2}}{P_{\perp}^{2}}-\frac{3}{2}+\ln\frac{P_{\perp}^{2}}{m_{e}^{2}}\Big{]}$
(7)

with $m_{e}$ the electron mass and $\alpha_{e}$ the QED coupling. Second, the
Fourier coefficient $\alpha_{s}C_{F}c_{n}$ is replaced by
$\alpha_{s}C_{F}c_{n}+\alpha_{e}c_{n}^{\gamma}$, where

$c_{n}^{\gamma}=(-1)^{n}\Big{[}\ln\frac{P_{\perp}^{2}}{m_{e}^{2}}+\frac{2}{\pi}\int^{\pi}_{0}d\phi(\pi-\phi)\frac{\cos\phi}{\sin\phi}(\cos
n\phi-1)\Big{]}~{}.$ (8)

We note that $c_{n}^{\gamma}$ is in line with $c_{n}(R)$ in the small-$R$
limit when $R$ is substituted with $m_{e}^{2}/P_{\perp}^{2}$, except for the
overall factor $(-1)^{n}$. Moreover, despite the fact that the QED coupling
$\alpha_{e}\approx 1/137$ is much smaller than that in QCD, the soft photon
contribution to the anisotropy is enhanced by the large logarithm
$\ln(P_{\perp}^{2}/m_{e}^{2})$ in the correlation region. The numerical impact
of the QED modifications is evaluated in the next section.

Figure 1:  (a) Harmonics in $e+p$ collisions with the inputs from the rcBK
solution, GBW model, and CT18A PDFs; (b) Harmonics predicted for $e+\text{Au}$
collisions with the inputs: the rcBK solution, GBW model, and EPPS21 PDFs; (c)
Comparison of the exact results and the small-$q_{\perp}$ asymptotic
behaviors. Solid lines: the exact results. Dashdotted lines: the
small-$q_{\perp}$ asymptotic expansions in Eq.(10).

Numeric analysis  To numerically study the harmonic coefficients $\langle\cos
n\phi\rangle$ in LJC, we employ the GBW model Golec-Biernat:1998zce  and the
solution of the running-coupling Balitsky-Kovchegov (rcBK) equation
Balitsky:1995ub ; Kovchegov:2006wf ; Kovchegov:1999yj ; Kovchegov:2006vj ;
Albacete:2010sy ; Golec-Biernat:2001dqn ; Albacete:2007yr ; Balitsky:2006wa ;
Gardi:2006rp ; Albacete:2007yr ; Balitsky:2007feb ; Berger:2010sh  for the
dipole S-matrix ${\cal S}_{x}(r_{\perp})$. For the gold nucleus target, the
corresponding saturation momentum is chosen as $Q_{s,A}^{2}\approx
5Q_{s,p}^{2}$, where $Q_{s,p}$ is the proton saturation momentum. In the GBW
model, ${\cal S}_{x}(r_{\perp})=e^{-r_{\perp}^{2}Q_{s}^{2}(x)/4}$,
$Q_{s,p}^{2}(x)=(x_{0}/x)^{0.28}~{}\text{GeV}^{2}$ with $x_{0}=3\times
10^{-4}$. The modefied MV model Fujii:2013gxa ; Albacete:2010sy , which gives
${\cal
S}_{x_{0}}(r_{\perp})=\exp\left[-(r_{\perp}^{2}Q_{s0}^{2})^{1.118}\ln(1/(0.24r_{\perp})+e)/4\right]$
with $Q_{s0,p}^{2}=0.16$ GeV2 at $x_{0}=0.01$, sets the initial conditions for
the rcBK evolution. As a comparison, we also compute $\langle\cos
n\phi\rangle$ in the TMD framework following Ref. Hatta:2021jcd  with
collinear parton distribution functions (PDFs) as inputs. Here we adopt the
NLO PDF sets of CT18A Hou:2019efy  and EPPS21 Eskola:2021nhw  parametrizations
for proton and the gold nucleus, respectively.

With the above inputs, we predict the $q_{\perp}$-distribution of the
$\langle\cos n\phi\rangle$ anisotropies for the following kinematics bin:
$\sqrt{s_{eN}}=89$ GeV, $y_{\ell}=2.41$, $0.008\leq x\leq 0.01,\ 4$ GeV$\leq
P_{\perp}\leq 4.4$ GeV, 5.6 GeV$\leq Q\leq$ 5.9 GeV. This choice of kinematics
coincides with the event simulation study of the LJC in Ref. Arratia:2019vju
for the EIC. Due to limited collisional energies, it is difficult to access
lower $x$ regions ($x\leq 1\times 10^{-3}$). Nevertheless, we later show that
the saturation effects can be amplified with the gold nuclear target. Here the
lower cut for $x$ is chosen as the minimum value allowed by the kinematic
constraints at given $s_{\text{eN}},P_{\perp}$, namely $x_{\text{min}}\approx
4P_{\perp}^{2}/s_{\textrm{eN}}$.

Figs. 1 (a) and (b) show $\langle\cos n\phi\rangle$ for $n=1,2,3$ as a
function of $q_{\perp}$ in different models with the proton and gold nucleus
targets. The results based on the rcBK solution and the GBW model are
presented in the solid and dotted lines, respectively. Here the predictions
are made at the jet cone size $R=0.4$.

Let us comment on the features of the harmonics $\langle\cos n\phi\rangle$
shown in Figs. 1 (a,b). First, the figures show that the anisotropies start
from zero with sharp rises in the small-$q_{\perp}$ region, and then their
increase becomes gradual in the larger-$q_{\perp}$ region. More specifically,
we find that the anisotropies follows a power law at small-$q_{\perp}$,
namely, $\langle\cos n\phi\rangle\sim q_{\perp}^{n}$ , see e.g., an
illustration in Fig. 1(c). As we show later, the small-$q_{\perp}$ behaviour
of the harmonics is an interesting tool for us to study the saturation effects
in LJC. Moreover, the anisotropy decreases quickly with the harmonic number
$n$. Since the dominant asymmetry arises due to the gluon emission from the
quark jet side, the first harmonic ($n=1$) is the largest coefficient. Due to
this hierarchy, we only plot the first three $n$. Nevertheless, our results
show that sizable anisotropies can be measured at the EIC at least for
harmonics with $n=1,2$. The above features are consistent with the collinear
framework results, see e.g., the dashed lines in Figs. 1 (a,b) and a previous
calculation for proton in Ref. Hatta:2021jcd .

Figure 2: QED modifications to the LJC harmonics with the input from the rcBK
solution for $e+\text{Au}$ collisions. Solid and dashed lines represent the
harmonics without and with QED modifications, respectively.

Furthermore, in Fig. 2, we show the effect of the QED correction to
$\langle\cos n\phi\rangle$ in dashed lines. Using the rcBK solution as the
input, one can find that the QED corrections are fairly apparent. Especially,
we find that the QED correction reduces the odd harmonics while it increases
the even harmonics. This is expected since the QED Bremsstrahlung radiation
favors the collinear region along the lepton direction, which is the away side
of the gluon radiation. Analytically, this feature is reflected by the sign
difference between the Fourier coefficient $c_{n}^{\gamma}$ and $c_{n}(R)$ for
the odd harmonics, as discussed in Eq. (8).

In fact, by comparing the $q_{\perp}$-distributions of the harmonics between
$e+p$ and $e+\text{Au}$ collisions in Figs. 1 (a) and (b), we observe a
sizable decrease of the anisotropy $\langle\cos n\phi\rangle$ in the small-$x$
formalism when switching from the proton target to the gold nucleus target.
Since the saturation effects are isotropic in this process, the anisotropies
are suppressed when the saturation momentum increases. To make this
suppression more transparent, we define the nuclear modification factor as
follows

$\displaystyle{R}_{eA}^{(n)}=\frac{\langle\cos n\phi\rangle_{eA}}{\langle\cos
n\phi\rangle_{ep}}~{}.$ (9)

As shown in Fig. 3, significant nuclear suppression is predicted for all three
harmonics in the low-$q_{\perp}$ region. On top of that, the saturation
formalism indicates a numerical hierarchy of suppression as the harmonic
number $n$ increases. This implies that the higher harmonics are more
sensitive to the parton saturation effects.

Figure 3: Nuclear modification factors of the harmonics. Solid and dashed
lines represent ${R}_{eA}^{(n)}$ from the rcBK solution without and with the
QED corrections, respectively; Bands represent ${R}_{eA}^{(n)}$ from EPPS21
nuclear PDFs without the QED corrections.

Also, we find that the QED modifications to the ratio ${R}_{eA}^{(n)}$ in Eq.
(9) are negligible, as we observe little differences between the solid and
dashed lines in Fig. 3. This observation can be understood as follows. For one
thing, although the QED correction to the Fourier coefficient,
$\alpha_{e}c_{n}^{\gamma}$, can have apparent contributions to the harmonics
as shown in Fig. 2, it cancels between $\langle\cos n\phi\rangle_{eA}$ and
$\langle\cos n\phi\rangle_{ep}$ in the ratio ${R}_{eA}^{(n)}$; For another,
the Sudakov form factor is dominated by the QCD contributions, which can be
checked in Eqs. (5,7).

In addition, one can conduct a comparative study by computing the above
nuclear modification factor ${R}_{eA}^{(n)}$ based on the inputs of the
nuclear PDFs. Within the collinear framework, one can encode various nuclear
effects in the nuclear PDFs by fitting them with existing experimental data.
The ratios $R^{(n)}_{eA}$ calculated with the EPPS21 gold nuclear PDFs
Eskola:2021nhw  and the CT18A proton PDFs Hou:2019efy  are presented in Fig. 3
with the error band at $90\%$ confidence-level. Here we neglect the
uncertainty from the baseline proton PDFs, since they are small.
Interestingly, we find that the predictions of $R^{(n)}_{eA}$ made with the
rcBK solution and EPPS21 PDFs are substantially different at small
$q_{\perp}$, although the predictions at large $q_{\perp}$ all start to
converge toward unity.

The nuclear modification factor plotted in Fig. 3 is the main result of this
work, since it illustrates the striking difference between the saturation and
non-saturation frameworks. Therefore, by measuring the LJC anisotropy in $e+p$
to $e+\text{Au}$ collisions, one can quantitatively study the parton
saturation and pinpoint its signature features at the upcoming EIC.

Analytic explanations of the numerical results.  Given the anisotropies
derived in correlation limit in Eq. (6), we further expand the expression in
the small-$q_{\perp}$ limit. Since saturation effects are expected to be
dominant in the small-$q_{\perp}$ region, this asymptotic expansion can
provide an analytic understanding of the nuclear suppression of the harmonics
due to parton saturation.

By approximating the Bessel functions
$J_{n}\left(q_{\perp}b_{\perp}\right)\sim{(q_{\perp}b_{\perp}/2)^{n}}/{\Gamma(n+1)}$
in Eq. (6), we can obtain the simple power law behavior for the harmonics,
which reads $\langle\cos n\phi\rangle\approx{\cal C}_{n}q_{\perp}^{n}$ with
${\cal C}_{n}$ defined as the ratio of two integrals over $b_{\perp}$. In the
correlation limit which gives $Q\gtrsim P_{\perp}\gg q_{\perp}$, one can
justify the use of the saddle-point approximation (see e.g., the applications
in Refs. Parisi:1979se ; Collins:1981va ; Collins:1984kg ; Shi:2021hwx ) in
the evaluation of these two integrals, and thus obtain the following
expression of the harmonics at small $q_{\perp}$

$\displaystyle\langle\cos n\phi\rangle$ $\displaystyle\approx{\cal
C}_{n}q_{\perp}^{n}$
$\displaystyle\approx\left(\frac{q_{\perp}b_{0}}{2\Lambda_{\text{QCD}}}\right)^{n}\frac{\alpha_{s}(\mu^{\text{sp}}_{n})C_{F}c_{n}(R)}{\pi
n\Gamma(n+1)}\frac{f_{q}(x,b^{\text{sp}}_{\perp
n})}{f_{q}(x,b^{\text{sp}}_{\perp 0})}~{}$
$\displaystyle\times\left[\frac{2\beta_{1}+C_{F}}{(n+2)\beta_{1}+C_{F}}\right]^{1+\frac{C_{F}}{2\beta_{1}}\ln\frac{e^{c_{0}(R)}Q^{4}}{\Lambda_{\text{QCD}}^{2}P_{\perp}^{2}}}~{},$
(10)

where $\beta_{1}=(33-2n_{f})/12$,
$\mu^{\text{sp}}_{n}=b_{0}/b^{\text{sp}}_{\perp n}$, $b_{0}\equiv
2e^{-\gamma_{E}}$. $b^{\text{sp}}_{\perp n}$ denotes the position of the
saddle point:

$\displaystyle b^{\text{sp}}_{\perp
n}=\frac{b_{0}}{\Lambda_{\text{QCD}}}\Bigg{[}\frac{e^{c_{0}(R)}Q^{4}}{\Lambda_{\text{QCD}}^{2}P_{\perp}^{2}}\Bigg{]}^{-\frac{C_{F}}{2(2+n)\beta_{1}+2C_{F}}}~{},$
(11)

which is mainly determined by the Sudakov factor in Eq.(5). In particular,
Eq.(11) in the $n=0$ case is akin to the saddle point originally found in Ref.
Collins:1984kg .

Fig. 1(c) presents the numerical comparison between the small-$q_{\perp}$
expansion of $\langle\cos n\phi\rangle$ in Eq. (10) and the exact expression
in Eq. (6) with the rcBK solution as the input for the gold target. One can
see that the two results coincide in the small-$q_{\perp}$ region.

Let us discuss the dependence of the saturation momentum $Q_{s}$ in the
harmonics, which has been entirely contained in the ratio of unintegrated
quark distribution $f_{q}(x,b_{\perp}^{\text{sp}})$ in Eq. (10). One finds
that $b^{\text{sp}}_{\perp n}\rightarrow 0$ in the limit $Q\gtrsim
P_{\perp}\rightarrow\infty$ from Eq. (11). Thus we can apply small-$b_{\perp}$
approximation on $f_{q}(x,b_{\perp}^{\text{sp}})$, which yields
$f_{q}(x,b_{\perp}^{\text{sp}})\propto
Q_{s}^{2}\ln[1/(Q_{s}b_{\perp}^{\text{sp}})]$ Mueller:1999wm ; Xiao:2017yya .
As a result, the anisotropic harmonics has the following asymptotic form,

$\displaystyle\langle\cos
n\phi\rangle\propto\frac{f_{q}(x,b^{\text{sp}}_{\perp
n})}{f_{q}(x,b^{\text{sp}}_{\perp
0})}\approx\frac{\ln(Q_{s}b^{\text{sp}}_{\perp
n})}{\ln(Q_{s}b^{\text{sp}}_{\perp 0})}~{}.$ (12)

Since Eq. (11) implies $b^{\text{sp}}_{\perp n}>b^{\text{sp}}_{\perp 0}$ for
positive $n$, one can check the derivative of $\langle\cos n\phi\rangle$ with
respective to $Q_{s}$ is negative. That means as the saturation momentum
$Q_{s}$ increases, the anisotropic harmonics $\langle\cos n\phi\rangle$
decrease. Therefore, we can qualitatively explain the suppression of the
anisotropy observed in the numerical study. As illustrated in Fig. 3, there is
a numerical hierarchy of the nuclear modifications, which the nuclear
suppression becomes stronger as the harmonic number $n$ increases. This can be
qualitatively explained by Eq. (12) and the fact that $b^{\text{sp}}_{\perp
n}$ increases with $n$ according to Eq. (11).

Conclusion  In summary, by investigating the azimuthal angle anisotropy of the
LJC in DIS in the correlation limit, we find that the LJC anisotropy is
sensitive to saturation effects and thus provides us a novel channel to study
the onset of the saturation phenomenon at the EIC. Through a comprehensive
numerical study, we demonstrate that there is a considerable suppression in
the anisotropy when the target changes from proton to the gold nucleus. This
nuclear suppression is significant in the small-$q_{\perp}$ region, where the
saturation effects are dominant. In contrast, we find much smaller nuclear
suppression in the collinear factorization framework. Therefore, the nuclear
modification factor $R_{eA}^{(n)}$ can be used as a benchmark to distinguish
the saturation and non-saturation frameworks.

Acknowledgments: We thank Feng Yuan for useful comments. This work is partly
supported by the CUHK-Shenzhen university development fund under grant No.
UDF01001859.

## References

  * (1) L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rept. 100, 1-150 (1983)
  * (2) A. H. Mueller and J. w. Qiu, Nucl. Phys. B 268, 427-452 (1986)
  * (3) A. H. Mueller, Nucl. Phys. B 335, 115-137 (1990)
  * (4) L. D. McLerran and R. Venugopalan, Phys. Rev. D 49, 2233-2241 (1994) [arXiv:hep-ph/9309289 [hep-ph]].
  * (5) L. D. McLerran and R. Venugopalan, Phys. Rev. D 49, 3352-3355 (1994) [arXiv:hep-ph/9311205 [hep-ph]].
  * (6) L. D. McLerran and R. Venugopalan, Phys. Rev. D 50, 2225-2233 (1994) [arXiv:hep-ph/9402335 [hep-ph]].
  * (7) F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan, Ann. Rev. Nucl. Part. Sci. 60, 463-489 (2010) [arXiv:1002.0333 [hep-ph]].
  * (8) E. Iancu and R. Venugopalan, [arXiv:hep-ph/0303204 [hep-ph]].
  * (9) D. Boer, M. Diehl, R. Milner, R. Venugopalan, W. Vogelsang, D. Kaplan, H. Montgomery, S. Vigdor, A. Accardi and E. C. Aschenauer, et al. [arXiv:1108.1713 [nucl-th]].
  * (10) A. Accardi, J. L. Albacete, M. Anselmino, N. Armesto, E. C. Aschenauer, A. Bacchetta, D. Boer, W. K. Brooks, T. Burton and N. B. Chang, et al. Eur. Phys. J. A 52, no.9, 268 (2016) [arXiv:1212.1701 [nucl-ex]].
  * (11) R. Abdul Khalek, A. Accardi, J. Adam, D. Adamiak, W. Akers, M. Albaladejo, A. Al-bataineh, M. G. Alexeev, F. Ameli and P. Antonioli, et al. Nucl. Phys. A 1026, 122447 (2022) [arXiv:2103.05419 [physics.ins-det]].
  * (12) R. Abdul Khalek, U. D’Alesio, M. Arratia, A. Bacchetta, M. Battaglieri, M. Begel, M. Boglione, R. Boughezal, R. Boussarie and G. Bozzi, et al. [arXiv:2203.13199 [hep-ph]].
  * (13) A. Morreale and F. Salazar, Universe 7, no.8, 312 (2021) [arXiv:2108.08254 [hep-ph]].
  * (14) F. Dominguez, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 106, 022301 (2011) [arXiv:1009.2141 [hep-ph]].
  * (15) F. Dominguez, C. Marquet, B. W. Xiao and F. Yuan, Phys. Rev. D 83, 105005 (2011) [arXiv:1101.0715 [hep-ph]].
  * (16) A. H. Mueller, B. W. Xiao and F. Yuan, Phys. Rev. D 88, no.11, 114010 (2013) [arXiv:1308.2993 [hep-ph]].
  * (17) A. Metz and J. Zhou, Phys. Rev. D 84, 051503 (2011) [arXiv:1105.1991 [hep-ph]].
  * (18) F. Dominguez, J. W. Qiu, B. W. Xiao and F. Yuan, Phys. Rev. D 85, 045003 (2012) [arXiv:1109.6293 [hep-ph]].
  * (19) A. Dumitru, T. Lappi and V. Skokov, Phys. Rev. Lett. 115, no.25, 252301 (2015) [arXiv:1508.04438 [hep-ph]].
  * (20) A. Dumitru and V. Skokov, Phys. Rev. D 94, no.1, 014030 (2016) doi:10.1103/PhysRevD.94.014030 [arXiv:1605.02739 [hep-ph]].
  * (21) D. Boer, P. J. Mulders, C. Pisano and J. Zhou, JHEP 08, 001 (2016) doi:10.1007/JHEP08(2016)001 [arXiv:1605.07934 [hep-ph]].
  * (22) A. Dumitru, V. Skokov and T. Ullrich, Phys. Rev. C 99, no.1, 015204 (2019) [arXiv:1809.02615 [hep-ph]].
  * (23) H. Mäntysaari, N. Mueller, F. Salazar and B. Schenke, Phys. Rev. Lett. 124, no.11, 112301 (2020) [arXiv:1912.05586 [nucl-th]].
  * (24) Y. Y. Zhao, M. M. Xu, L. Z. Chen, D. H. Zhang and Y. F. Wu, Phys. Rev. D 104, no.11, 114032 (2021) [arXiv:2105.08818 [hep-ph]].
  * (25) R. Boussarie, H. Mäntysaari, F. Salazar and B. Schenke, JHEP 09, 178 (2021) [arXiv:2106.11301 [hep-ph]].
  * (26) P. Caucal, F. Salazar and R. Venugopalan, JHEP 11, 222 (2021) [arXiv:2108.06347 [hep-ph]].
  * (27) Y. Y. Zhang and X. N. Wang, Phys. Rev. D 105, no.3, 034015 (2022) [arXiv:2104.04520 [hep-ph]].
  * (28) P. Taels, T. Altinoluk, G. Beuf and C. Marquet, JHEP 10, 184 (2022) [arXiv:2204.11650 [hep-ph]].
  * (29) P. Caucal, F. Salazar, B. Schenke and R. Venugopalan, [arXiv:2208.13872 [hep-ph]].
  * (30) R. Boussarie, A. V. Grabovsky, L. Szymanowski and S. Wallon, JHEP 09, 026 (2014) [arXiv:1405.7676 [hep-ph]].
  * (31) R. Boussarie, A. V. Grabovsky, L. Szymanowski and S. Wallon, JHEP 11, 149 (2016) [arXiv:1606.00419 [hep-ph]].
  * (32) F. Salazar and B. Schenke, Phys. Rev. D 100, no.3, 034007 (2019) [arXiv:1905.03763 [hep-ph]].
  * (33) R. Boussarie, A. V. Grabovsky, L. Szymanowski and S. Wallon, Phys. Rev. D 100, no.7, 074020 (2019) [arXiv:1905.07371 [hep-ph]].
  * (34) D. Boer and C. Setyadi, Phys. Rev. D 104, no.7, 074006 (2021) [arXiv:2106.15148 [hep-ph]].
  * (35) E. Iancu, A. H. Mueller and D. N. Triantafyllopoulos, Phys. Rev. Lett. 128, no.20, 202001 (2022) [arXiv:2112.06353 [hep-ph]].
  * (36) E. Iancu, A. H. Mueller, D. N. Triantafyllopoulos and S. Y. Wei, JHEP 10, 103 (2022) [arXiv:2207.06268 [hep-ph]].
  * (37) Y. Hatta, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 116, no.20, 202301 (2016) [arXiv:1601.01585 [hep-ph]].
  * (38) T. Altinoluk, N. Armesto, G. Beuf and A. H. Rezaeian, Phys. Lett. B 758, 373-383 (2016) [arXiv:1511.07452 [hep-ph]].
  * (39) H. Mäntysaari, N. Mueller and B. Schenke, Phys. Rev. D 99, no.7, 074004 (2019) [arXiv:1902.05087 [hep-ph]].
  * (40) Y. Hagiwara, C. Zhang, J. Zhou and Y. j. Zhou, Phys. Rev. D 104, no.9, 094021 (2021) [arXiv:2106.13466 [hep-ph]].
  * (41) L. Zheng, E. C. Aschenauer, J. H. Lee and B. W. Xiao, Phys. Rev. D 89, no.7, 074037 (2014) [arXiv:1403.2413 [hep-ph]].
  * (42) F. Bergabo and J. Jalilian-Marian, Nucl. Phys. A 1018, 122358 (2022) [arXiv:2108.10428 [hep-ph]].
  * (43) F. Bergabo and J. Jalilian-Marian, Phys. Rev. D 106, no.5, 054035 (2022) [arXiv:2207.03606 [hep-ph]].
  * (44) E. Iancu and Y. Mulian, [arXiv:2211.04837 [hep-ph]].
  * (45) M. Fucilla, A. V. Grabovsky, E. Li, L. Szymanowski and S. Wallon, [arXiv:2211.05774 [hep-ph]].
  * (46) I. Kolbé, K. Roy, F. Salazar, B. Schenke and R. Venugopalan, JHEP 01, 052 (2021) [arXiv:2008.04372 [hep-ph]].
  * (47) J. Jalilian-Marian and Y. V. Kovchegov, Phys. Rev. D 70, 114017 (2004) [erratum: Phys. Rev. D 71, 079901 (2005)] [arXiv:hep-ph/0405266 [hep-ph]].
  * (48) D. Kharzeev, E. Levin and L. McLerran, Nucl. Phys. A 748, 627-640 (2005) [arXiv:hep-ph/0403271 [hep-ph]].
  * (49) C. Marquet, Nucl. Phys. A 796, 41-60 (2007) [arXiv:0708.0231 [hep-ph]].
  * (50) K. Tuchin, Nucl. Phys. A 846, 83-94 (2010) [arXiv:0912.5479 [hep-ph]].
  * (51) A. Dumitru and J. Jalilian-Marian, Phys. Rev. D 82, 074023 (2010) [arXiv:1008.0480 [hep-ph]].
  * (52) K. Kutak and S. Sapeta, Phys. Rev. D 86, 094043 (2012) [arXiv:1205.5035 [hep-ph]].
  * (53) A. van Hameren, P. Kotko, K. Kutak and S. Sapeta, Phys. Lett. B 737, 335-340 (2014) [arXiv:1404.6204 [hep-ph]].
  * (54) P. Kotko, K. Kutak, C. Marquet, E. Petreska, S. Sapeta and A. van Hameren, JHEP 09, 106 (2015) [arXiv:1503.03421 [hep-ph]].
  * (55) A. van Hameren, P. Kotko, K. Kutak, C. Marquet, E. Petreska and S. Sapeta, JHEP 12, 034 (2016) [erratum: JHEP 02, 158 (2019)] [arXiv:1607.03121 [hep-ph]].
  * (56) A. van Hameren, P. Kotko, K. Kutak and S. Sapeta, Phys. Lett. B 795, 511-515 (2019) [arXiv:1903.01361 [hep-ph]].
  * (57) A. van Hameren, P. Kotko, K. Kutak and S. Sapeta, Phys. Lett. B 814, 136078 (2021) [arXiv:2010.13066 [hep-ph]].
  * (58) C. Marquet, E. Petreska and C. Roiesnel, JHEP 10, 065 (2016) [arXiv:1608.02577 [hep-ph]].
  * (59) S. R. Klein and H. Mäntysaari, Nature Rev. Phys. 1, no.11, 662-674 (2019) doi:10.1038/s42254-019-0107-6 [arXiv:1910.10858 [hep-ex]].
  * (60) E. Iancu and Y. Mulian, JHEP 03, 005 (2021) [arXiv:2009.11930 [hep-ph]].
  * (61) A. D. Bolognino, F. G. Celiberto, M. Fucilla, D. Y. Ivanov and A. Papa, Phys. Rev. D 103, no.9, 094004 (2021) [arXiv:2103.07396 [hep-ph]].
  * (62) M. A. Al-Mashad, A. van Hameren, H. Kakkad, P. Kotko, K. Kutak, P. van Mechelen and S. Sapeta, [arXiv:2210.06613 [hep-ph]].
  * (63) Y. Hagiwara, Y. Hatta, R. Pasechnik, M. Tasevsky and O. Teryaev, Phys. Rev. D 96, no.3, 034009 (2017) [arXiv:1706.01765 [hep-ph]].
  * (64) P. Kotko, K. Kutak, S. Sapeta, A. M. Stasto and M. Strikman, Eur. Phys. J. C 77, no.5, 353 (2017) [arXiv:1702.03063 [hep-ph]].
  * (65) S. Bhattacharya, A. Metz, V. K. Ojha, J. Y. Tsai and J. Zhou, Phys. Lett. B 833, 137383 (2022) [arXiv:1802.10550 [hep-ph]].
  * (66) R. Boussarie, Y. Hatta, B. W. Xiao and F. Yuan, Phys. Rev. D 98, no.7, 074015 (2018) [arXiv:1807.08697 [hep-ph]].
  * (67) J. L. Albacete and C. Marquet, Phys. Rev. Lett. 105, 162301 (2010) [arXiv:1005.4065 [hep-ph]].
  * (68) A. Stasto, B. W. Xiao and F. Yuan, Phys. Lett. B 716, 430-434 (2012) [arXiv:1109.1817 [hep-ph]].
  * (69) T. Lappi and H. Mantysaari, Nucl. Phys. A 908, 51-72 (2013) [arXiv:1209.2853 [hep-ph]].
  * (70) E. Iancu and J. Laidet, Nucl. Phys. A 916, 48-78 (2013) [arXiv:1305.5926 [hep-ph]].
  * (71) J. L. Albacete, G. Giacalone, C. Marquet and M. Matas, Phys. Rev. D 99, no.1, 014002 (2019) [arXiv:1805.05711 [hep-ph]].
  * (72) A. Stasto, S. Y. Wei, B. W. Xiao and F. Yuan, Phys. Lett. B 784, 301-306 (2018) [arXiv:1805.05712 [hep-ph]].
  * (73) J. Jalilian-Marian, Nucl. Phys. A 770, 210-220 (2006) [arXiv:hep-ph/0509338 [hep-ph]].
  * (74) J. Jalilian-Marian and A. H. Rezaeian, Phys. Rev. D 86, 034016 (2012) [arXiv:1204.1319 [hep-ph]].
  * (75) A. Stasto, B. W. Xiao and D. Zaslavsky, Phys. Rev. D 86, 014009 (2012) [arXiv:1204.4861 [hep-ph]].
  * (76) A. H. Rezaeian, Phys. Rev. D 86, 094016 (2012) [arXiv:1209.0478 [hep-ph]].
  * (77) E. Basso, V. P. Goncalves, J. Nemchik, R. Pasechnik and M. Sumbera, Phys. Rev. D 93, no.3, 034023 (2016) [arXiv:1510.00650 [hep-ph]].
  * (78) A. H. Rezaeian, Phys. Rev. D 93, no.9, 094030 (2016) [arXiv:1603.07354 [hep-ph]].
  * (79) E. Basso, V. P. Goncalves, M. Krelina, J. Nemchik and R. Pasechnik, Phys. Rev. D 93, no.9, 094027 (2016) [arXiv:1603.01893 [hep-ph]].
  * (80) D. Boer, P. J. Mulders, J. Zhou and Y. j. Zhou, JHEP 10, 196 (2017) [arXiv:1702.08195 [hep-ph]].
  * (81) S. Benić and A. Dumitru, Phys. Rev. D 97, no.1, 014012 (2018) [arXiv:1710.01991 [hep-ph]].
  * (82) V. P. Goncalves, Y. Lima, R. Pasechnik and M. Šumbera, Phys. Rev. D 101, no.9, 094019 (2020) [arXiv:2003.02555 [hep-ph]].
  * (83) S. Benić, O. Garcia-Montero and A. Perkov, Phys. Rev. D 105, no.11, 114052 (2022) [arXiv:2203.01685 [hep-ph]].
  * (84) D. Boer, Y. Hagiwara, J. Zhou and Y. j. Zhou, Phys. Rev. D 105, no.9, 096017 (2022) [arXiv:2203.00267 [hep-ph]].
  * (85) F. Gelis and J. Jalilian-Marian, Phys. Rev. D 66, 094014 (2002) [arXiv:hep-ph/0208141 [hep-ph]].
  * (86) A. Kovner and A. H. Rezaeian, Phys. Rev. D 90, no.1, 014031 (2014) [arXiv:1404.5632 [hep-ph]].
  * (87) A. Kovner and A. H. Rezaeian, Phys. Rev. D 92, no.7, 074045 (2015) [arXiv:1508.02412 [hep-ph]].
  * (88) C. Marquet, S. Y. Wei and B. W. Xiao, Phys. Lett. B 802, 135253 (2020) [arXiv:1909.08572 [hep-ph]].
  * (89) E. Akcakaya, A. Schäfer and J. Zhou, Phys. Rev. D 87, no.5, 054010 (2013) [arXiv:1208.4965 [hep-ph]].
  * (90) C. Marquet, C. Roiesnel and P. Taels, Phys. Rev. D 97, no.1, 014004 (2018) [arXiv:1710.05698 [hep-ph]].
  * (91) X. Liu, F. Ringer, W. Vogelsang and F. Yuan, Phys. Rev. D 102, no.9, 094022 (2020) [arXiv:2007.12866 [hep-ph]].
  * (92) X. Liu, F. Ringer, W. Vogelsang and F. Yuan, Phys. Rev. Lett. 122, no.19, 192003 (2019) [arXiv:1812.08077 [hep-ph]].
  * (93) M. Arratia, Y. Song, F. Ringer and B. V. Jacak, Phys. Rev. C 101, no.6, 065204 (2020) [arXiv:1912.05931 [nucl-ex]].
  * (94) Y. Hatta, B. W. Xiao, F. Yuan and J. Zhou, Phys. Rev. D 104, no.5, 054037 (2021) [arXiv:2106.05307 [hep-ph]].
  * (95) Y. Hatta, B. W. Xiao, F. Yuan and J. Zhou, Phys. Rev. Lett. 126 (2021) no.14, 142001 [arXiv:2010.10774 [hep-ph]].
  * (96) S. Catani, M. Grazzini and A. Torre, Nucl. Phys. B 890, 518-538 (2014) [arXiv:1408.4564 [hep-ph]].
  * (97) S. Catani, M. Grazzini and H. Sargsyan, JHEP 06, 017 (2017) [arXiv:1703.08468 [hep-ph]].
  * (98) W. L. Ju and M. Schönherr, [arXiv:2210.09272 [hep-ph]].
  * (99) C. Marquet, B. W. Xiao and F. Yuan, Phys. Lett. B 682, 207-211 (2009) [arXiv:0906.1454 [hep-ph]].
  * (100) B. W. Xiao and F. Yuan, Phys. Rev. D 82, 114009 (2010) [arXiv:1008.4432 [hep-ph]].
  * (101) B. W. Xiao and F. Yuan, Phys. Rev. Lett. 105, 062001 (2010) [arXiv:1003.0482 [hep-ph]].
  * (102) A. H. Mueller, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 110, no.8, 082301 (2013) [arXiv:1210.5792 [hep-ph]].
  * (103) L. D. McLerran and R. Venugopalan, Phys. Rev. D 59, 094002 (1999) [arXiv:hep-ph/9809427 [hep-ph]].
  * (104) R. Venugopalan, Acta Phys. Polon. B 30, 3731-3761 (1999) [arXiv:hep-ph/9911371 [hep-ph]].
  * (105) A. H. Mueller, Nucl. Phys. B 558, 285-303 (1999) [arXiv:hep-ph/9904404 [hep-ph]].
  * (106) B. W. Xiao, F. Yuan and J. Zhou, Nucl. Phys. B 921, 104-126 (2017) [arXiv:1703.06163 [hep-ph]].
  * (107) P. Sun, J. Isaacson, C. P. Yuan and F. Yuan, Int. J. Mod. Phys. A 33, no.11, 1841006 (2018) [arXiv:1406.3073 [hep-ph]].
  * (108) A. Prokudin, P. Sun and F. Yuan, Phys. Lett. B 750, 533-538 (2015) [arXiv:1505.05588 [hep-ph]].
  * (109) K. J. Golec-Biernat and M. Wusthoff, Phys. Rev. D 59, 014017 (1998) [arXiv:hep-ph/9807513 [hep-ph]].
  * (110) I. Balitsky, Nucl. Phys. B 463, 99-160 (1996) [arXiv:hep-ph/9509348 [hep-ph]].
  * (111) Y. V. Kovchegov and H. Weigert, Nucl. Phys. A 789, 260-284 (2007) [arXiv:hep-ph/0612071 [hep-ph]].
  * (112) Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999) [arXiv:hep-ph/9901281 [hep-ph]].
  * (113) Y. V. Kovchegov and H. Weigert, Nucl. Phys. A 784, 188-226 (2007) [arXiv:hep-ph/0609090 [hep-ph]].
  * (114) J. L. Albacete, N. Armesto, J. G. Milhano, P. Quiroga-Arias and C. A. Salgado, Eur. Phys. J. C 71, 1705 (2011) [arXiv:1012.4408 [hep-ph]].
  * (115) K. J. Golec-Biernat, L. Motyka and A. M. Stasto, Phys. Rev. D 65, 074037 (2002) [arXiv:hep-ph/0110325 [hep-ph]].
  * (116) J. L. Albacete and Y. V. Kovchegov, Phys. Rev. D 75, 125021 (2007) [arXiv:0704.0612 [hep-ph]].
  * (117) I. Balitsky, Phys. Rev. D 75, 014001 (2007) [arXiv:hep-ph/0609105 [hep-ph]].
  * (118) E. Gardi, J. Kuokkanen, K. Rummukainen and H. Weigert, Nucl. Phys. A 784, 282-340 (2007) [arXiv:hep-ph/0609087 [hep-ph]].
  * (119) I. Balitsky and G. A. Chirilli, Phys. Rev. D 77, 014019 (2008) [arXiv:0710.4330 [hep-ph]].
  * (120) J. Berger and A. Stasto, Phys. Rev. D 83, 034015 (2011) [arXiv:1010.0671 [hep-ph]].
  * (121) H. Fujii and K. Watanabe, Nucl. Phys. A 915, 1-23 (2013) [arXiv:1304.2221 [hep-ph]].
  * (122) T. J. Hou, J. Gao, T. J. Hobbs, K. Xie, S. Dulat, M. Guzzi, J. Huston, P. Nadolsky, J. Pumplin and C. Schmidt, et al. Phys. Rev. D 103, no.1, 014013 (2021) [arXiv:1912.10053 [hep-ph]].
  * (123) K. J. Eskola, P. Paakkinen, H. Paukkunen and C. A. Salgado, Eur. Phys. J. C 82, no.5, 413 (2022) [arXiv:2112.12462 [hep-ph]].
  * (124) G. Parisi and R. Petronzio, Nucl. Phys. B 154, 427-440 (1979)
  * (125) J. C. Collins and D. E. Soper, Nucl. Phys. B 197, 446-476 (1982)
  * (126) J. C. Collins, D. E. Soper and G. F. Sterman, Nucl. Phys. B 250, 199-224 (1985)
  * (127) Y. Shi, L. Wang, S. Y. Wei and B. W. Xiao, Phys. Rev. Lett. 128, no.20, 202302 (2022) [arXiv:2112.06975 [hep-ph]].


* Dupuis [1988] Dupuis, P.  Large deviations analysis of some recursive algorithms with state dependent noise.  _The Annals of Probability_ , 16(4):1509–1536, October 1988.
  * Dupuis & Kushner [1985] Dupuis, P. and Kushner, H. J.  Stochastic approximations via large deviations: Asymptotic properties.  _SIAM Journal on Control and Optimization_ , 23(5):675–696, September 1985.
  * Eberle [2016] Eberle, A.  Reflection couplings and contraction rates for diffusions.  _Probability Theory and Related Fields_ , 166(3-4):851–886, December 2016.
  * Erdogdu et al. [2018] Erdogdu, M. A., Mackey, L. W., and Shamir, O.  Global non-convex optimization with discretized diffusions.  In _NeurIPS ’18: Proceedings of the 32nd International Conference of Neural Information Processing Systems_ , 2018.
  * Feng et al. [2020] Feng, Y., Gao, T., Li, L., Liu, J.-G., and Lu, Y.  Uniform-in-time weak error analysis for stochastic gradient descent algorithms via diffusion approximation.  _Communications in Mathematical Sciences_ , 18(1):163–188, 2020.
  * Freidlin & Wentzell [1998] Freidlin, M. I. and Wentzell, A. D.  _Random Perturbations of Dynamical Systems_.  Springer, 1 edition, 1998.
  * Ge et al. [2015] Ge, R., Huang, F., Jin, C., and Yuan, Y.  Escaping from saddle points – Online stochastic gradient for tensor decomposition.  In _COLT ’15: Proceedings of the 28th Annual Conference on Learning Theory_ , 2015.
  * Gulinsky & Veretennikov [1993] Gulinsky, O. V. and Veretennikov, A. Y.  _Large Deviations for Discrete-Time Processes with Averaging_.  De Gruyter, 1993.
  * Gürbüzbalaban et al. [2021] Gürbüzbalaban, M., Şimşekli, U., and Zhu, L.  The heavy-tail phenomenon in SGD.  In _ICML ’21: Proceedings of the 38th International Conference on Machine Learning_ , 2021.
  * Hernández-Lerma & Lasserre [2003] Hernández-Lerma, O. and Lasserre, J.-B.  _Markov Chains and Invariant Probabilities_.  Birkhäuser, Basel, 2003.
  * Hodgkinson & Mahoney [2021] Hodgkinson, L. and Mahoney, M.  Multiplicative noise and heavy tails in stochastic optimization.  In _ICML ’21: Proceedings of the 38th International Conference on Machine Learning_ , 2021.
  * Hsieh et al. [2021] Hsieh, Y.-P., Mertikopoulos, P., and Cevher, V.  The limits of min-max optimization algorithms: Convergence to spurious non-critical sets.  In _ICML ’21: Proceedings of the 38th International Conference on Machine Learning_ , 2021.
  * Hsieh et al. [2023] Hsieh, Y.-P., Karimi, M. R., Krause, A., and Mertikopoulos, P.  Riemannian stochastic optimization methods avoid strict saddle points.  In _NeurIPS ’23: Proceedings of the 37th International Conference on Neural Information Processing Systems_ , 2023.
  * Hu et al. [2019] Hu, W., Li, C. J., Li, L., and Liu, J.-G.  On the diffusion approximation of nonconvex stochastic gradient descent.  _Annals of Mathematical Sciences and Applications_ , 4(1):3–32, 2019.
  * Ioffe et al. [1979] Ioffe, A. D., Tikhomirov, V. M., and Makowski, K.  _Theory of Extremal Problems_.  North Holland, Amsterdam, NL, 1979.
  * Jastrzębski et al. [2017] Jastrzębski, S., Kenton, Z., Arpit, D., Ballas, N., Fischer, A., Bengio, Y., and Storkey, A.  Three factors influencing minima in SGD.  https://arxiv.org/abs/1711.04623, 2017.
  * Jin et al. [2017] Jin, C., Ge, R., Netrapalli, P., Kakade, S. M., and Jordan, M. I.  How to escape saddle points efficiently.  In _ICML ’17: Proceedings of the 34th International Conference on Machine Learning_ , 2017.
  * Kallenberg [2021] Kallenberg, O.  _Foundations of Modern Probability_ , volume 99 of _Probability Theory and Stochastic Modelling_.  Springer, 2021.
  * Kawaguchi [2016] Kawaguchi, K.  Deep learning without poor local minima.  In _NIPS ’16: Proceedings of the 30th International Conference on Neural Information Processing Systems_ , 2016.
  * Kiefer & Wolfowitz [1952] Kiefer, J. and Wolfowitz, J.  Stochastic estimation of the maximum of a regression function.  _The Annals of Mathematical Statistics_ , 23(3):462–466, 1952.
  * Kifer [1988] Kifer, Y.  _Random Perturbations of Dynamical Systems_.  Birkhäuser, Boston, MA, 1988.
  * Kifer [1990] Kifer, Y.  A discrete-time version of the Wentzell-Freidlin theory.  _Annals of Probability_ , 18(4):1676–1692, October 1990.
  * Lan [2020] Lan, G.  _First-order and Stochastic Optimization Methods for Machine Learning_.  Springer, 2020.
  * Landau & Lifshitz [1976] Landau, L. D. and Lifshitz, E. M.  Statistical physics.  In _Course of Theoretical Physics_ , volume 5. Pergamon Press, Oxford, 1976.
  * Lee et al. [2019] Lee, J. D., Panageas, I., Piliouras, G., Simchowitz, M., Jordan, M. I., and Recht, B.  First-order methods almost always avoid strict saddle points.  _Mathematical Programming_ , 176(1):311–337, February 2019.
  * Li & Wang [2022a] Li, L. and Wang, Y.  A sharp uniform-in-time error estimate for stochastic gradient Langevin dynamics.  https://arxiv.org/abs/2207.09304, 2022a.
  * Li & Wang [2022b] Li, L. and Wang, Y.  On uniform-in-time diffusion approximation for stochastic gradient descent.  https://arxiv.org/abs/2207.04922, 2022b.
  * Li et al. [2017] Li, Q., Tai, C., and E, W.  Stochastic modified equations and adaptive stochastic gradient algorithms.  In _ICML ’17: Proceedings of the 34th International Conference on Machine Learning_ , 2017.
  * Li et al. [2019] Li, Q., Tai, C., and E, W.  Stochastic modified equations and dynamics of stochastic gradient algorithms I: Mathematical foundations.  _Journal of Machine Learning Research_ , 20(40):1–47, 2019.
  * Li et al. [2020] Li, Z., Lyu, K., and Arora, S.  Reconciling modern deep learning with traditional optimization analyses: The intrinsic learning rate.  In _NeurIPS ’20: Proceedings of the 34th International Conference on Neural Information Processing Systems_ , 2020.
  * Li et al. [2021] Li, Z., Wang, T., and Arora, S.  What happens after SGD reaches zero loss? –A mathematical framework.  In _ICLR ’21: Proceedings of the 2021 International Conference on Learning Representations_ , 2021.
  * Li et al. [2022] Li, Z., Wang, T., and Yu, D.  Fast mixing of stochastic gradient descent with normalization and weight decay.  In _NeurIPS ’22: Proceedings of the 36th International Conference on Neural Information Processing Systems_ , 2022.
  * Liu et al. [2020] Liu, S., Papailiopoulos, D., and Achlioptas, D.  Bad global minima exist and SGD can reach them.  In _NeurIPS ’20: Proceedings of the 34th International Conference on Neural Information Processing Systems_ , 2020.
  * Ljung [1977] Ljung, L.  Analysis of recursive stochastic algorithms.  _IEEE Trans. Autom. Control_ , 22(4):551–575, August 1977.
  * Lu et al. [2021] Lu, Y., Balasubramanian, K., Volgushev, S., and Erdogdu, M. A.  An analysis of constant step size SGD in the non-convex regime: Asymptotic normality and bias.  In _NeurIPS ’21: Proceedings of the 35th International Conference on Neural Information Processing Systems_ , 2021.
  * Lytras & Mertikopoulos [2024] Lytras, I. and Mertikopoulos, P.  Tamed Langevin sampling under weaker conditions.  https://arxiv.org/abs/2405.17693, 2024.
  * Majka et al. [2020] Majka, M. B., Mijatović, A., and Szpruch, Ł.  Nonasymptotic bounds for sampling algorithms without log-concavity.  _The Annals of Applied Probability_ , 30(4):1534–1581, 2020.
  * Mertikopoulos et al. [2020] Mertikopoulos, P., Hallak, N., Kavis, A., and Cevher, V.  On the almost sure convergence of stochastic gradient descent in non-convex problems.  In _NeurIPS ’20: Proceedings of the 34th International Conference on Neural Information Processing Systems_ , 2020.
  * Mertikopoulos et al. [2024] Mertikopoulos, P., Hsieh, Y.-P., and Cevher, V.  A unified stochastic approximation framework for learning in games.  _Mathematical Programming_ , 203:559–609, January 2024\.
  * Mignacco & Urbani [2022] Mignacco, F. and Urbani, P.  The effective noise of stochastic gradient descent.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2022(8):083405, 2022.
  * Mignacco et al. [2020] Mignacco, F., Krzakala, F., Urbani, P., and Zdeborová, L.  Dynamical mean-field theory for stochastic gradient descent in Gaussian mixture classification.  In _NeurIPS ’20: Proceedings of the 34th International Conference on Neural Information Processing Systems_ , 2020.
  * Mori et al. [2022] Mori, T., Ziyin, L., Liu, K., and Ueda, M.  Power-law escape rate of SGD.  In _ICML ’22: Proceedings of the 39th International Conference on Machine Learning_ , 2022.
  * Pavasovic et al. [2023] Pavasovic, K. L., Durmus, A., and Şimşekli, U.  Approximate heavy tails in offline (multi-pass) stochastic gradient descent.  In _NeurIPS ’23: Proceedings of the 37th International Conference on Neural Information Processing Systems_ , 2023.
  * Pemantle [1990] Pemantle, R.  Nonconvergence to unstable points in urn models and stochastic aproximations.  _Annals of Probability_ , 18(2):698–712, April 1990.
  * Raginsky et al. [2017] Raginsky, M., Rakhlin, A., and Telgarsky, M.  Non-convex learning via stochastic gradient Langevin dynamics: A nonasymptotic analysis.  In _COLT ’17: Proceedings of the 30th Annual Conference on Learning Theory_ , 2017.
  * Rigollet & Hütter [2023] Rigollet, P. and Hütter, J.-C.  High-dimensional statistics.  https://arxiv.org/abs/2310.19244, 2023.
  * Robbins & Monro [1951] Robbins, H. and Monro, S.  A stochastic approximation method.  _Annals of Mathematical Statistics_ , 22:400–407, 1951\.
  * Tutte [2001] Tutte, W. T.  _Graph theory_.  Cambridge University Press, 2001.
  * van den Dries & Miller [1996] van den Dries, L. and Miller, C.  Geometric categories and $o$-minimal structures.  _Duke Mathematical Journal_ , 84(2), August 1996.
  * Veiga et al. [2024] Veiga, R., Remizova, A., and Macris, N.  Stochastic gradient flow dynamics of test risk and its exact solution for weak features.  https://arxiv.org/abs/2402.07626, 2024.
  * Wainwright [2019] Wainwright, M. J.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge University Press, 2019.
  * Wang & Wang [2022] Wang, Y. and Wang, Z.  Three-stage evolution and fast equilibrium for SGD with non-degerate critical points.  In _ICML ’22: Proceedings of the 39th International Conference on Machine Learning_ , 2022.
  * Wojtowytsch [2023] Wojtowytsch, S.  Stochastic gradient descent with noise of machine learning type. Part I: Discrete time analysis.  _Journal of Nonlinear Science_ , 33(3):45, 2023\.
  * Xie et al. [2021] Xie, Z., Sato, I., and Sugiyama, M.  A diffusion theory for deep learning dynamics: Stochastic gradient descent exponentially favors flat minima.  In _ICLR ’21: Proceedings of the 2021 International Conference on Learning Representations_ , 2021.
  * Yang et al. [2021] Yang, J., Hu, W., and Li, C. J.  On the fast convergence of random perturbations of the gradient flow.  _Asymptotic Analysis_ , 122(3-4):371–393, 2021\.
  * Zhou et al. [2020] Zhou, Z., Mertikopoulos, P., Bambos, N., Boyd, S. P., and Glynn, P. W.  On the convergence of mirror descent beyond stochastic convex programming.  _SIAM Journal on Optimization_ , 30(1):687–716, 2020.
  * Ziyin et al. [2023] Ziyin, L., Li, H., and Ueda, M.  Law of balance and stationary distribution of stochastic gradient descent.  https://arxiv.org/abs/2308.06671, 2023.

  *[SGD]: stochastic gradient descent
  *[SDE]: stochastic differential equation
  *[SGLD]: stochastic gradient Langevin dynamics
  *[SME]: stochastic modified equation
  *[SFO]: stochastic first-order oracle
  *[i.i.d.]: independent and identically distributed
  *[LDP]: large deviation principle
  *[CGFs]: cumulant-generating function
  *[RHS]: right-hand side
  *[DMFT]: dynamic mean-field theory
  *[l.s.c.]: lower semi-continuous
  *[ODE]: ordinary differential equation
  *[u.s.c.]: upper semi-continuous
  *[r.v.]: random variable


# Quantifying the Carbon Emissions
of Machine Learning

Alexandre Lacoste
Element AI
<EMAIL_ADDRESS>
Alexandra Luccioni11footnotemark: 1
Mila, Université de Montréal
<EMAIL_ADDRESS>
Victor Schmidt11footnotemark: 1
Mila, Université de Montréal
<EMAIL_ADDRESS>
Thomas Dandres
Polytechnique Montréal, CIRAIG
<EMAIL_ADDRESS>
equal contribution

###### Abstract

From an environmental standpoint, there are a few crucial aspects of training
a neural network that have a major impact on the quantity of carbon that it
emits. These factors include: the location of the server used for training and
the energy grid that it uses, the length of the training procedure, and even
the make and model of hardware on which the training takes place. In order to
approximate these emissions, we present our _Machine Learning Emissions
Calculator_ , a tool for our community to better understand the environmental
impact of training ML models. We accompany this tool with an explanation of
the factors cited above, as well as concrete actions that individual
practitioners and organizations can take to mitigate their carbon emissions.

##  1 Introduction

While a decade ago, only a few ML pioneers were training neural networks on
GPUs (Graphical Processing Units), in recent years powerful GPUs have become
increasingly accessible and used by ML practitioners worldwide. Furthermore,
new models often need to beat existing challenges, which entails training on
more GPUs, with larger datasets, for a longer time. This expansion brings with
it ever-growing costs in terms of the energy needed to fuel it. This trend has
been the subject of recent studies aiming to evaluate the climate impact of
AI, which have predominantly put the focus on the environmental cost of
training large-scale models connected to grids powered by fossil fuels [1, 2].
While these models are not necessarily representative of common practice, we
believe that it is important to continue this conversation further and work
towards defining the tools and steps that we need to assess the carbon
emissions generated by the models we train, as well as to propose ways to
reduce those emissions.

In this work, we present our Machine Learning Emissions Calculator
(https://mlco2.github.io/impact/), a tool for our community to estimate the
amount of carbon emissions produced by training ML models. We accompany this
tool with a presentation of key concepts and an explanation of the factors
impacting emissions. Finally, we end our article with some recommendations of
best practices for the overall ML research community, as well as for
individual researchers.

##  2 Quantifying Carbon Emissions in Neural Network Training

In order to quantify carbon emissions, we use CO2-equivalents (CO2eq), which
is a standardized measure used to express the global-warming potential of
various greenhouse gases as a single number, i.e. as the amount of CO2 which
would have the equivalent global warming impact [3]. We will use this single
metric to compare the factors and choices that impact overall amount of
emissions produced by training an ML model in the sections below.

###  2.1 Type of Energy Used

Practically speaking, it is hard to estimate exactly the amount of CO2eq
emitted by a cloud server in a given location because the information
regarding the energy grid that it is connected to is rarely publicly
available. However, if we assume that all servers are connected to local grids
at their physical location, we are able to make an estimation of the amount of
CO2eq that they emit using public data sources [4, 5] . Therefore, in order to
create our emissions calculator, we gathered data regarding CO2eq emissions of
different grid locations and cross-referenced them with known GPU server
locations from the three major cloud providers: Google Cloud Platform,
Microsoft Azure and Amazon Web Services 111The data can be found at
https://github.com/mlco2/impact/tree/master/data. Our aim in doing this is to
illustrate the degree of variability that exists depending on the location of
a given server. For instance, in Figure 1, we show the distribution and
variation in carbon emissions depending on geographical region. It can be
noted that a large amount of variation can be found within a single region;
for instance, servers located in North America can emit anywhere between 20g
CO2eq/kWh in Quebec, Canada to 736.6g CO2eq/kWh in Iowa, USA [5].

Figure 1: Variation of the Average Carbon Intensity of Servers Worldwide, by
Region. (Vertical bars represent regions with a single available data point.)

###  2.2 Computing Infrastructure and Training Time

Another, more subtle, factor in carbon emitted by a neural network is the
computing infrastructure used and training time of the model. In terms of
performance the number of floating point operations per second (FLOPS) of GPUs
has been steadily increasing in recent years, from 100 Giga FLOPS per second
in 2004 to up to 15 Tera FLOPS per second in recent hardware [6]. However,
with neural network architectures becoming deeper and more complex, recent
state-of-the-art models are often trained on multiples GPUs for several weeks
(or months) to beat benchmark performance, requiring more and more energy.

Finally, when it comes to defining a training procedure for ML architectures,
there are several elements to consider: for starters, whether it is necessary
to train a model from scratch or whether fine-tuning is adequate for the task
at hand. Notably, recent research has shown that using pre-trained models with
task-specific fine-tuning performs as well as training from scratch, while
being more robust, for tasks in image recognition [7, 8] and NLP [9].
Furthermore, when it comes to hyperparameter search, it has been proven both
empirically and theoretically that random hyperparameter search is more
efficient than grid search for hyperparameter optimization [10], and there is
much research being done on ways to improve the efficiency of hyperparameter
optimization [11, 12], which makes it possible to continue choosing the right
hyperparameters for new models without incurring superfluous computing and
energy costs.

##  3 ML Emissions Calculator and Actionable Items

It is difficult to provide clear-cut guidelines for ML researchers to follow
in order to reduce the carbon emissions, or specific benchmarks for the
training time that a given model or task warrants. Nonetheless, we think that
there are certain best practices and actionable items that can be adopted by
our community to reduce environmental impact of the ML domain. We present some
of these, along with our ML emissions calculator, in the current section.

#### Quantify Your Emissions

Being informed regarding the factors that impact the quantity of carbon
emissions produced by ML research is the first step to making positive
changes. It is for this reason that we created our _ML Emissions Calculator_.
This tool, currently in its alpha version, takes as input the details
regarding the training of an ML model: the geographical zone of the server,
the type of GPU, and the training time, and gives as output the approximate
amount of CO2eq produced. We collected publicly available data for the 4 main
variables of this computation: (i) the energy consumption of hardware (see
"Choose More Efficient Hardware" below), (ii) the location of providers’
regions of compute – which we assumed to be connected to their local grid,
(iii) the region’s CO2eq emissions per kWh and (iv) potential offsets bought
by the provider.

We intend to adopt an open and transparent approach: the data we used is
publicly available, debatable and editable through Github issues and pull
requests. We are therefore open to updating data as more information becomes
available. Since this paper’s core goal is to raise awareness around the
carbon emissions of ML, we have also included two educational sections in the
website: one about learning the main notions and concepts related to this
domain (e.g. RECs, carbon neutrality, etc.), the other about actionable items
an individual or an organization can leverage to mitigate their carbon impact.

#### Choose Cloud Providers Wisely

In recent years, many cloud providers have defined ambitious sustainability
goals and are offsetting their emissions through Renewable Energy Certificates
(RECs) in an effort to become carbon neutral, a term used to indicate a net
zero carbon footprint of an organization. Each REC bought attests that 1 MWh
of renewable energy has been added to the energy grid and can be used to
offset an equivalent amount of non-renewable energy. For instance, Google
Cloud Platform is certified carbon neutral and funds solar and wind farms
directly on local grids through RECs [13]. Microsoft Azure is also certified
carbon neutral and 44% of its electricity consumption directly comes from
renewable energy, according to a 2016 estimate [14]. Finally, to the best of
our knowledge, while not yet 100% carbon neutral on an organizational level,
Amazon Web Services is also funding renewable energy projects and some of
their data centers are powered by renewable energy [15].

Another major energy consumption factor of server installations is the power
usage effectiveness (PUE) of the centers where the GPUs are hosted, which
represents the percentage of energy consumption that is used for cooling,
power conversion, and other auxiliary tasks, and can vary immensely. For
example, Google Cloud Services has an average PUE of 1.1, meaning that only
11% of their total energy usage is not used for the servers themselves, a
ratio that they have been steadily reducing using Reinforcement Learning [16,
17]. Finally, if you rely on a local private compute infrastructure, it is
also possible to engage with administrators about quantifying and offsetting
the emissions produced, as well as improving the efficiency of your grid –
this may help bring your organization toward carbon neutrality and have a
significant impact at scale.

#### Select Data Center Location

While many cloud providers are carbon neutral, some of their data centers may
still be carbon intensive due to the local grid that they are connected to,
whereas others will be low carbon and powered solely by renewable energy
sources. Hence, selecting the data center location where an algorithm will be
trained has a large impact on its direct carbon emissions. This choice can be
achieved by consciously selecting the server location before dispatching your
jobs. As we illustrated in previous sections, this single choice can make the
direct emissions of an algorithm vary by a factor of 40, from 20g CO2eq/kWh in
a location that uses renewable energy sources to 820g CO2eq/kWh in a location
that solely relies on fossil fuels [5]. For a model such as VGG [18] or BERT
[19], which are trained on multiple GPUs for several weeks, this can
correspond to avoiding emitting several hundreds of kilograms of CO2eq by
training on a server powered by hydroelectricity instead of fossil fuels.

#### Reduce Wasted Ressources

Grid search is still often used in practice, in spite of its low efficiency
both in terms of model performance and environmental impact. However, it has
been shown that random search (and others) not only is a straightforward
replacement but also has potential to significantly accelerate hyperparameter
search [20, 21, 22], consequently reducing carbon emissions. Also, while
failed experiments are a common part of ML research and are sometimes
unavoidable, their number can often be reduced with careful design such as
unit tests, integration tests, and extensive and early debugging.
Uninformative experiments are also frequent (sometimes unknowingly) – they can
be caused by unstable learning algorithms requiring averaging results over
many random seeds. Taking the time to carry out a literature review and to
understand the potential sources of noise before launching large-scale
hyperparameter searches increases the chance of obtaining reproducible and
statistically significant results. Hence, reducing the need to extend the
experiment cycles.

#### Choose More Efficient Hardware

The choice of computing hardware can also have a major impact on ML emissions.
To perform a comparison between different devices, their compute efficiency
can be estimated in FLOPS/W. This estimation is based on their theoretical
peak performance with respect to their Thermal Design Power (TDP)222Empirical
measurement of GFLOPS/W on various ML architecture would provide more accurate
numbers but we are only interested in approximate values to compare classes of
devices.. Using this approach, it can be found that CPUs can be 10 times less
efficient than GPUs while TPU 3 can be 4 to 8 times more efficient than GPUs
[23] (refer to Table Appendix B: Hardware Efficiency for details).
Interestingly, in contexts where low power consumption and efficiency are
important, e.g., for embedded applications, GPUs such as the Jetson AGX Xavier
can be 10 to 20 times more efficient than traditional GPUs.

##  4 Discussion

The factors that we discussed in the current work give ML practitioners a
certain amount of control over the environmental impact produced by the
training of their models. We are aware that these choices are not always
possible to make in practice – for instance, the choice of server location can
be limited due to privacy considerations in the case of applications in the
medical or financial domain, and large amounts of data may be needed to
produce most robust models. However, we find that our emissions calculator is
a good starting point to estimate the impact that small choices in model
training can have on direct carbon emissions resulting from ML research.

Despite our best efforts, our calculator remains simply an approximation of
the true emissions produced by ML training for several reasons: to start with,
there is the issue of global load balancing, i.e. if a majority of
practitioners choose to run their models in a low-carbon location, the servers
will get saturated and other servers will still need be used. In that
perspective, the global gain will not be a 40-fold reduction of emissions, but
much smaller. Furthermore, there is a lack of transparency with regards to the
true quantity of emissions produced by organizations, so while we use the
current best publicly-available sources, there is still a large margin of
error with regards to the exact quantity of energy consumed and carbon
produced – we remain open to additional data sources and numbers. Finally,
while in the current version of our tool, we focus on quantifying the
emissions of training ML models, there is still the issue of deploying them,
since the inference process is also energy-expensive, especially if done
continuously and on a large scale. This is something that should be taken into
account by ML practitioners in their products that are deployed in real-world
settings, for instance by using energy-efficient architectures [24] and
computing infrastructure.

There are also more far-reaching discussions to be had regarding the
environmental value of scientific knowledge in general and of ML research in
particular. On one hand, there is valuable research to be done in ML
especially with regards to tackling climate change [25, 26], whereas on the
other hand, the emissions of the field of ML are growing quickly [1]. We do
not propose the solution to this problem, but we believe that there are steps
to be taken, for instance by using efficiency as an evaluation criterion (as
proposed by [2]) or by taking concrete steps to reduce emissions (as proposed
by the current paper). We hope that our work, along with others, will open the
door for these conversations and debates to take place, to quantify the
environmental impact of our field, and for positive changes that can be made
to reduce it.

## References

  * [1] Emma Strubell, Ananya Ganesh, and Andrew McCallum.  Energy and policy considerations for deep learning in nlp.  arXiv preprint arXiv:1906.02243, 2019.
  * [2] Roy Schwartz, Jesse Dodge, Noah A Smith, and Oren Etzioni.  Green ai.  arXiv preprint arXiv:1907.10597, 2019.
  * [3] Simon Eggleston, Leandro Buendia, Kyoko Miwa, Todd Ngara, and Kiyoto Tanabe.  2006 IPCC guidelines for national greenhouse gas inventories, volume 5.  Institute for Global Environmental Strategies Hayama, Japan, 2006.
  * [4] WM To and Peter KC Lee.  Ghg emissions from electricity consumption: A case study of hong kong from 2002 to 2015 and trends to 2030.  Journal of cleaner production, 165:589–598, 2017.
  * [5] Matthew Brander, Aman Sood, Charlotte Wylie, Amy Haughton, and Jessica Lovell.  Electricity-specific emission factors for grid electricity.  Ecometrica, Emissionfactors. com, 2011.
  * [6] Ari Harju, Topi Siro, Filippo Federici Canova, Samuli Hakala, and Teemu Rantalaiho.  Computational physics on graphics processing units.  In Proceedings of the 11th international conference on Applied Parallel and Scientific Computing, pages 3–26. Springer-Verlag, 2012.
  * [7] Nima Tajbakhsh, Jae Y Shin, Suryakanth R Gurudu, R Todd Hurst, Christopher B Kendall, Michael B Gotway, and Jianming Liang.  Convolutional neural networks for medical image analysis: Full training or fine tuning?  IEEE transactions on medical imaging, 35(5):1299–1312, 2016.
  * [8] Keiji Yanai and Yoshiyuki Kawano.  Food image recognition using deep convolutional network with pre-training and fine-tuning.  In 2015 IEEE International Conference on Multimedia & Expo Workshops (ICMEW), pages 1–6. IEEE, 2015.
  * [9] Jeremy Howard and Sebastian Ruder.  Universal language model fine-tuning for text classification.  arXiv preprint arXiv:1801.06146, 2018.
  * [10] James Bergstra and Yoshua Bengio.  Random search for hyper-parameter optimization.  Journal of Machine Learning Research, 13(Feb):281–305, 2012.
  * [11] Stefan Falkner, Aaron Klein, and Frank Hutter.  Bohb: Robust and efficient hyperparameter optimization at scale.  arXiv preprint arXiv:1807.01774, 2018.
  * [12] Matthias Feurer and Frank Hutter.  Hyperparameter optimization.  In Automated Machine Learning, pages 3–33. Springer, 2019.
  * [13] Google.  Google environmental report 2018, 2018.
  * [14] Microsoft.  Beyond carbon neutral. white paper, 2018.
  * [15] Amazon Web Services.  Aws & sustainability, 2019.
  * [16] Jim Gao.  Machine learning applications for data center optimization, 2014.
  * [17] Google Data Centers efficiency: How we do it.  https://www.google.com/about/datacenters/efficiency/internal/, 2019\.  Accessed: 2019-08-23.
  * [18] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [19] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805, 2018.
  * [20] Lisha Li, Kevin Jamieson, Giulia DeSalvo, Afshin Rostamizadeh, and Ameet Talwalkar.  Hyperband: A novel bandit-based approach to hyperparameter optimization.  arXiv preprint arXiv:1603.06560, 2016.
  * [21] Liam Li, Kevin Jamieson, Afshin Rostamizadeh, Ekaterina Gonina, Moritz Hardt, Benjamin Recht, and Ameet Talwalkar.  Massively parallel hyperparameter tuning.  arXiv preprint arXiv:1810.05934, 2018.
  * [22] Stefan Falkner, Aaron Klein, and Frank Hutter.  Bohb: Robust and efficient hyperparameter optimization at scale.  arXiv preprint arXiv:1807.01774, 2018.
  * [23] Paul Teich.  Tearing apart google’s tpu 3.0 ai coprocessor.  https://www.nextplatform.com/2018/05/10/tearing-apart-googles-tpu-3-0-ai-coprocessor/, 2018\.
  * [24] Ermao Cai, Da-Cheng Juan, Dimitrios Stamoulis, and Diana Marculescu.  Neuralpower: Predict and deploy energy-efficient convolutional neural networks.  arXiv preprint arXiv:1710.05420, 2017.
  * [25] David Rolnick, Priya L Donti, Lynn H Kaack, Kelly Kochanski, Alexandre Lacoste, Kris Sankaran, Andrew Slavin Ross, Nikola Milojevic-Dupont, Natasha Jaques, Anna Waldman-Brown, et al.  Tackling climate change with machine learning.  arXiv preprint arXiv:1906.05433, 2019.
  * [26] Victor Schmidt, Alexandra Luccioni, S. Karthik Mukkavilli, Narmada Balasooriya, Kris Sankaran, Jennifer Chayes, and Yoshua Bengio.  Visualizing the consequences of climate change using cycle-consistent adversarial networks.  CoRR, abs/1905.03709, 2019.

## Appendix A: Energy Grid Data Used for the ML Emissions Calculator

For clarity purposes, the data presented in this appendix contains fewer
columns than what can be found in our public database:
https://github.com/mlco2/impact/tree/master/data. For instance we did not
include the sources in these tables. For future reference, the data’s commit
hash at the time of publication is e692e28.

### Google Cloud Platform

Region | Country | City | Estimated gCO2e/kWh)
---|---|---|---
asia-east1 | Taiwan | Changhua County | 557
asia-east2 | China | Hong Kong | 702
asia-northeast1 | Japan | Tokyo | 516
asia-northeast2 | Japan | Osaka | 516
asia-south1 | India | Mumbai | 920
asia-southeast1 | Singapore | Jurong West | 419
australia-southeast1 | Australia | Sydney | 802
europe-north1 | Finland | Hamina | 211
europe-west1 | Belgium | St. Ghislain | 267
europe-west2 | United Kingdom | London | 623
europe-west3 | Germany | Frankfurt | 615
europe-west4 | Netherlands | Eemshaven | 569
europe-west6 | Switzerland | Zürich | 16
northamerica-northeast1 | Canada | Montréal | 20
southamerica-east1 | Brazil | São Paulo | 205
us-central1 | USA | Council Bluffs | 566.3
us-east1 | USA | Moncks Corner | 367.8
us-east4 | USA | Ashburn | 367.8
us-west1 | USA | The Dalles | 297.6
us-west2 | USA | Los Angeles | 240.6

### Amazon Web Services

Region | Country | City | gCO2e/kWh
---|---|---|---
us-east-2 | USA | Columbus | 568.2
us-east-1 | USA | Ashburn | 367.8
us-west-1 | USA | San Francisco | 240.6
us-west-2 | USA | Portland | 297.6
ap-east-1 | China | Hong Kong | 702
ap-south-1 | India | Mumbai | 920
ap-northeast-3 | Japan | Osaka | 516
ap-northeast-2 | South Korea | Seoul | 517
ap-southeast-1 | Singapore | Singapore | 419
ap-southeast-2 | Australia | Sydney | 802
ap-northeast-1 | Japan | Tokyo | 516
ca-central-1 | Canada | Montreal | 20
cn-north-1 | China | Beijing | 680
cn-northwest-1 | China | Zhongwei | 680
eu-central-1 | Germany | Frankfurt am Main | 615
eu-west-1 | Ireland | Dublin | 617
eu-west-2 | United Kingdom | London | 623
eu-west-3 | France | Paris | 105
eu-north-1 | Sweden | Stockholm | 47
sa-east-1 | Brazil | Sao Paulo | 205
us-gov-east-1 | USA | Dublin | 568.2
us-gov-west-1 | USA | Seattle | 297.6

### Microsoft Azure

Region | Country | City | gCO2e/kWh
---|---|---|---
eastasia | Hong Kong | Wan Chai | 702
southeastasia | Singapore | Singapore | 419
centralus | USA | Des Moines | 736.6
eastus | USA | Blue Ridge | 367.8
eastus2 | USA | Boydton | 367.8
westus | USA | San Francisco | 240.6
northcentralus | USA | Chicago | 568.2
southcentralus | USA | San Antonio | 460.4
northeurope | Ireland | Dublin | 617
westeurope | Netherlands | Amsterdam | 569
japanwest | Japan | Osaka-shi | 516
japaneast | Japan | Tokyo | 516
brazilsouth | Brazil | Sao Paulo | 205
australiaeast | Australia | Sydney | 802
australiasoutheast | Australia | Melbourne | 805
southindia | India | Pallavaram | 920
centralindia | India | Lohogaon | 920
westindia | India | Mumbai | 920
canadacentral | Canada | Toronto | 69.3
canadaeast | Canada | Quebec | 20
uksouth | United Kingdom | Midhurst | 623
ukwest | United Kingdom | Wallasey | 623
westcentralus | USA | Mountain View | 297.6
westus2 | USA | Quincy | 297.6
koreacentral | South Korea | Seoul | 517
koreasouth | South Korea | Busan | 517
francecentral | France | Huriel | 105
francesouth | France | Realmont | 105
australiacentral | Australia | Forrest | 900
australiacentral2 | Australia | Forrest | 900
southafricanorth | South Africa | Pretoria | 1009
southafricawest | South Africa | Stellenbosch | 1009

## Appendix B: Hardware Efficiency

Name | Watt (TDP) | TFLOPS32 | TFLOPS16 | GFLOPS32/W | GFLOPS16/W
---|---|---|---|---|---
RTX 2080 Ti | 250 | 13.45 | 26.90 | 53.80 | 107.60
RTX 2080 | 215 | 10.00 | 20.00 | 46.51 | 93.02
GTX 1080 Ti | 250 | 11.34 | 0.17 | 45.36 | 0.68
GTX 1080 | 180 | 8.00 | 0.13 | 44.44 | 0.72
AMD RX480 | 150 | 5.80 | 5.80 | 38.67 | 38.67
Titan V | 250 | 14.90 | 29.80 | 59.60 | 119.20
Tesla V100 | 300 | 15.00 | 30.00 | 50.00 | 100.00
TPU2 | 250 | 22.00 | 45.00 | 88.00 | 180.00
TPU3 | 200 | 45.00 | 90.00 | 225.00 | 450.00
Intel Xeon E5-2699 | 145 | 0.70 | 0.70 | 4.83 | 4.83
AGX Xavier | 30 | 16.00 | 32.00 | 533.33 | 1066.67