# A Better Angle on Hadron Transverse Momentum Distributions at the EIC

Anjie Gao<EMAIL_ADDRESS>Johannes K. L. Michel<EMAIL_ADDRESS>Iain W.
Stewart<EMAIL_ADDRESS>Zhiquan Sun<EMAIL_ADDRESS>Center for Theoretical
Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

(September 22, 2022)

###### Abstract

We propose an observable $q_{*}$ sensitive to transverse momentum dependence
(TMD) in $eN\to ehX$, with $q_{*}/E_{N}$ defined purely by lab-frame angles.
In 3D measurements of confinement and hadronization this resolves the
crippling issue of accurately reconstructing small transverse momentum
$P_{hT}$. We prove factorization for $\mathrm{d}\sigma_{h}/\mathrm{d}q_{*}$
for $q_{*}\ll Q$ with standard TMD functions, enabling $q_{*}$ to substitute
for $P_{hT}$. A double-angle reconstruction method is given which is exact to
all orders in QCD for $q_{*}\ll Q$. $q_{*}$ enables an order-of-magnitude
improvement in the expected experimental resolution at the EIC.

††preprint:  MIT-CTP 5462

##  I Introduction

A deeper understanding of the emergent properties of the nucleon, such as
confinement and hadronization, has been a frontier of nuclear and particle
physics research since the inception of Quantum Chromodynamics (QCD) five
decades ago. An important one-dimensional view of the nucleon is provided by
the deep-inelastic scattering (DIS) process $e^{-}(\ell)+N(P)\to
e^{-}(\ell^{\prime})+X$, where the scattering is mediated by an off-shell
photon of momentum $q=\ell-\ell^{\prime}$ (with $Q^{2}\equiv-q^{2}>0$).
Confinement is probed by measurements of $x=Q^{2}/(2P\cdot q)$, the momentum
fraction carried by the colliding parton inside the nucleon $N$. A more
intricate view is obtained by identifying a hadron $h$ in semi-inclusive DIS
(SIDIS),

$\displaystyle e^{-}(\ell)+N(P)\to e^{-}(\ell^{\prime})+h(P_{h})+X\,.$ (1)

Here measurements of the longitudinal momentum fraction $z=(P\cdot
P_{h})/(P\cdot q)$ that the hadron retains when forming from the struck quark
give insight into the complex dynamics of hadronization. Measuring the
hadron’s _transverse_ momentum $\vec{P}_{hT}$ relative to $\vec{q}$ gives
access to a three-dimensional view of the confinement and hadronization
processes for $N$ and $h$, together with spin correlations that probe these
processes.

The region most sensitive to these dynamics occurs for small transverse
momentum, $P_{hT}\sim\Lambda_{\mathrm{QCD}}\ll Q$, where
$\Lambda_{\mathrm{QCD}}$ is the QCD confinement scale. Here the cross section
obeys a rigorous factorization theorem Collins (2011), with the confinement
and hadronization dynamics encoded in universal transverse momentum-dependent
(TMD) parton distribution functions (PDFs) and fragmentation functions (FFs).
SIDIS cross sections have been extensively studied experimentally at
HERMESAirapetian  _et al._ (2013, 2019, 2020), COMPASSAlekseev  _et al._
(2009); Aghasyan  _et al._ (2018); Parsamyan (2018), RHICAschenauer  _et al._
(2015); Adamczyk  _et al._ (2016), and JLabAvakian  _et al._ (2004); Jawalkar
_et al._ (2018); Moran  _et al._ (2022) and together with the Drell-Yan
process have enabled extractions of TMD PDFs and FFs (TMDs) by various groups,
e.g. Scimemi and Vladimirov (2020); Bacchetta  _et al._ (2020); Bury  _et al._
(2021). A key scientific goal of the upcoming Electron-Ion-Collider (EIC)Abdul
Khalek  _et al._ (2022) is to study SIDIS with enormous beam luminosities to
determine TMD PDFs and FFs with unprecedented precision. Progress has also
been made towards calculations of TMD PDFs from lattice QCD Ebert  _et al._
(2019); Shanahan  _et al._ (2020, 2021); Schlemmer  _et al._ (2021); Li  _et
al._ (2022); Chu  _et al._ (2022).

A key challenge in experimental studies of TMDs is that measurements of
$\vec{P}_{hT}$ require reconstructing the photon momentum (or, equivalently,
the so-called Breit frame) to great accuracy to avoid loss of precision on
$P_{hT}=\lvert\vec{P}_{hT}\rvert\ll Q$. More precisely, a misreconstruction of
$\vec{\ell}^{\prime}$ by $\mathcal{O}(\Delta)$ leads to a misreconstruction of
$\vec{q}$ and therefore $\vec{P}_{hT}$ by $\mathcal{O}(\Delta)$, which is a
large uncertainty for small $P_{hT}\ll Q$. For example, for a nominal
measurement at $P_{hT}/z=1\,{\rm GeV}$ with $Q=20\,{\rm GeV}$, a typical
detector resolution of $\Delta=0.5\,{\rm GeV}$ leads to a 50% uncertainty.
This puts in peril the EIC physics program to unveil the dynamics of
hadronization and confinement in the kinematic region with the largest
sensitivity.

In this paper we construct a novel SIDIS observable, $q_{*}$, designed to be
maximally resilient against resolution effects while delivering the same
sensitivity to TMD dynamics as $\vec{P}_{hT}$. The key insight is that while
the magnitude of the electron and hadron three momentum is subject to limited
detector resolution and requires careful calibration, modern tracking
detectors deliver near-perfect resolution on the _angles_ of charged particle
tracks. We will therefore construct $q_{*}$ to satisfy the following three
criteria: (i) it is purely defined in terms of lab-frame angles and the beam
energies; (ii) at small values $q_{*}\ll Q$, the differential cross section
$\mathrm{d}\sigma/\mathrm{d}q_{*}$, including all possible spin correlations,
still satisfies a rigorous factorization theorem in terms of the standard TMD
PDFs and FFs; (iii) it does not dilute the statistical power of the available
event sample. Our construction is inspired by, but features key differences
to, the Drell-Yan $\phi^{*}_{\eta}$ observable in hadron-hadron collisions
Banfi  _et al._ (2011a), which has enabled extraordinarily successful
measurements of the TMD region at the permil level at the Tevatron and at the
LHC Abazov  _et al._ (2011); Aaij  _et al._ (2013); Aad  _et al._ (2016); Aaij
_et al._ (2015); Aad  _et al._ (2020). A novel feature of $q_{*}$ compared to
$\phi^{*}_{\eta}$ is that it can be negative, allowing one to measure
asymmetries $\mathrm{d}\sigma(q_{*}>0)-\mathrm{d}\sigma(q_{*}<0)$ to obtain
direct access to spin correlations.

Below we define $q_{*}$ in detail, prove the factorization theorem for $q_{*}$
with standard TMDs, and evaluate the expected detector resolution, statistical
power, and resilience against systematic biases of $q^{*}$ versus
$\vec{P}_{hT}$.

##  II Constructing $q_{*}$

Figure 1: (a) Definition of $\phi_{\mathrm{acop}}^{\mathrm{rest}}$ needed to
construct $q^{*}$. Momenta are not drawn to scale. We define the conventional
Trento frame Bacchetta  _et al._ (2004) for SIDIS, as well as the target rest
frame and the EIC frame. $\phi_{h}$ is the azimuthal separation between
$\vec{P}_{h}$ and $\vec{\ell}$ in the Trento frame. The acoplanarity angle in
the target rest frame is
$\phi_{\mathrm{acop}}^{\mathrm{rest}}\equiv\pi-\Delta\phi^{\mathrm{rest}}$,
where $\Delta\phi^{\mathrm{rest}}$ is the azimuthal separation between
$\vec{P}_{h}$ and $\vec{\ell}^{\,\prime}$ in the target rest frame. (b) In-
plane geometry for leading power kinematics $P_{h,\perp}/Q\ll 1$, where we can
approximate that $\vec{P}_{h}$ is along the same direction as $\vec{q}$. This
geometry gives rise to the double angle formula in Eq. (II), which allows us
to determine DIS kinematic variables $Q$, $x$, $y$ with angular measurements.

Consider the target rest frame shown in Fig. 1a where the nucleus $N$ is at
rest and the $z$-axis is along the incoming lepton beam. The lepton momenta
$\vec{\ell}$ and $\vec{\ell}^{\prime}$ define the lepton plane as the $x$-$z$
plane. We wish to take advantage of the high-precision reconstruction of polar
angles (equivalent to rapidities) and azimuthal angles in the EIC lab frame.
Here we give results in terms of EIC frame rapidities in the light target mass
limit $M\ll Q$. Results with full $M$ dependence are given in Supplement IV.1.
The acoplanarity angle in the target rest frame,
$\phi_{\mathrm{acop}}^{\mathrm{rest}}$, is defined by
$\tan\phi_{\mathrm{acop}}^{\mathrm{rest}}=-P_{h,y}/P_{h,x}$, where $P_{h,x}$
and $P_{h,y}$ are components of $P_{h}$. From the illustration of the 3D
kinematics in Fig. 1 a, it is obvious that
$\tan\phi_{\mathrm{acop}}^{\mathrm{rest}}$ $\propto\sin\phi_{h}\,P_{hT}$,
where $\phi_{h}$ is the azimuthal angle of $\vec{P}_{h}$ in the Trento frame.
We may thus use $\phi_{\mathrm{acop}}$ as a precision probe of the hadron
transverse momentum $P_{hT}$.

To work out the full relation between $\phi_{\mathrm{acop}}^{\mathrm{rest}}$
and $P_{hT}$, consider now the leading-power (LP) kinematics illustrated in
Fig. 1 b, where $\lambda\sim P_{hT}/(zQ)\ll 1$. We find

$\displaystyle\tan\phi_{\mathrm{acop}}^{\mathrm{rest}}$
$\displaystyle=\frac{\sin\phi_{h}\,P_{hT}}{zQ\sqrt{1-y}}+\mathcal{O}(\lambda^{2})\,.$
(2)

We now wish to express $Q$ and $y=(P\cdot q)/(P\cdot\ell)$ in terms of final-
state angles in the EIC frame, which is defined by a $180^{\circ}$ rotation
about our rest frame $y$ axis and then a boost along the $z$-axis, so
$\phi_{\mathrm{acop}}^{\mathrm{EIC}}=-\phi_{\mathrm{acop}}^{\mathrm{rest}}$.
From Fig. 1 b, momentum conservation gives
$\ell_{x^{\prime}}=\ell^{\prime}_{x^{\prime}}$, $q_{x}=-\ell^{\prime}_{x}$,
and $P_{hT}\ll Q$ implies $\theta_{h}+\theta_{e}+\alpha=\pi/2$. We find
$y=1-\sin\theta_{h}/\cos\alpha$ and
$Q^{2}=(\ell^{0}_{\mathrm{rest}})^{2}[\frac{\sin^{2}\theta_{e}}{\cos^{2}\alpha}-(1-\frac{\sin\theta_{h}}{\cos\alpha})^{2}]$.
Boosting to the lab frame and taking the $M\ll Q$ limit, we have:

$\displaystyle Q^{2}\\!$
$\displaystyle=\\!(2P^{0}_{\mathrm{EIC}})^{2}\frac{e^{\eta_{e}+\eta_{h}}}{1+e^{\Delta\eta}}\\!+\\!\mathcal{O}(\lambda)\,,\,\,\
y\\!=\\!\frac{1}{1+e^{\Delta\eta}}\\!+\\!\mathcal{O}(\lambda^{2})\,,$
$\displaystyle x$
$\displaystyle=\bigl{[}(2P^{0}_{\mathrm{EIC}})^{2}/s\bigr{]}\,e^{\eta_{e}+\eta_{h}}+{\cal
O}(\lambda)\,,$ (3)

where $\eta_{i}$ are the EIC frame pseudorapidities of the outgoing lepton
$i=e$ and hadron $i=h$, $\Delta\eta\equiv\eta_{h}-\eta_{e}$, and
$s=(P+\ell)^{2}$. This construction agrees with the double angle formula in
Ref. Bentvelsen  _et al._ (1992). However, Ref. Bentvelsen  _et al._ (1992)
use the struck quark angle in a tree-level picture, while our Eq. (II) uses
the hadron angle and holds to all orders in $\alpha_{s}$, and up to power
corrections in $P_{hT}/(zQ)$ which control the distance to the Born limit. The
${\cal O}(\lambda)$ corrections to Eq. (II) are given in Supplement IV.2, and
can be used to assess the size of these power corrections.

To exploit the proportionality in Eq. (2) to probe $P_{hT}$ more directly, we
define an optimized observable:

$\displaystyle q_{*}$ $\displaystyle\equiv
2P^{0}_{\mathrm{EIC}}\frac{e^{\eta_{h}}}{1+e^{\Delta\eta}}\tan\phi_{\mathrm{acop}}^{\mathrm{EIC}}\,.$
(4)

Expanding in $P_{hT}\ll zQ$ it has a simple leading power (LP) limit

$\displaystyle
q_{*}\stackrel{{\scriptstyle\text{LP}}}{{=}}-\sin\phi_{h}\frac{P_{hT}}{z}\,.$
(5)

Thus for TMD analyses, $Q^{2}$, $x$, $y$, and $q^{*}$ can all be measured from
the beam energy $P^{0}_{\mathrm{EIC}}$ and angular variables. We may also
define a dimensionless variable,

$\displaystyle\phi^{*}_{\mathrm{SIDIS}}$
$\displaystyle=\sqrt{\frac{e^{\Delta\eta}}{1+e^{\Delta\eta}}}\tan\phi_{\mathrm{acop}}^{\mathrm{EIC}}=\,\frac{q_{*}}{Q}+{\cal
O}(\lambda^{2})\,.$ (6)

This is analogous to the setup for the $\phi^{*}_{\eta}$ observable in Drell-
Yan Banfi  _et al._ (2011a). We expect the purely angular observables $q_{*}$
and $\phi^{*}_{\mathrm{SIDIS}}$ to be measured to much higher relative
precision compared to the transverse momentum $P_{hT}$.

##  III Factorization

Consider the standard factorization theorem for polarized SIDIS Ji  _et al._
(2005, 2004); Bacchetta  _et al._ (2007); Collins (2011),

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}^{2}\vec{P}_{hT}}=\sigma_{0}\Bigl{\\{}W_{UU,T}+\lambda_{e}S_{L}\sqrt{1-\epsilon^{2}}\,W_{LL}$
$\displaystyle\>\;+\epsilon\cos(2\phi_{h})W_{UU}^{\cos(2\phi_{h})}+S_{L}\epsilon\sin(2\phi_{h})W_{{UL}}^{\sin(2\phi_{h})}$
$\displaystyle\>\;+S_{T}\sin(\phi_{h}\\!-\\!\phi_{S})W_{UT,T}^{\sin(\phi_{h}\\!-\phi_{S})}+\epsilon
S_{T}\Bigl{[}\sin(\phi_{h}\\!+\\!\phi_{S})$ $\displaystyle\>\;\>\;\times
W_{UT}^{\sin(\phi_{h}\\!+\phi_{S})}\\!\\!+\sin(3\phi_{h}\\!-\\!\phi_{S})W_{UT}^{\sin(3\phi_{h}\\!-\phi_{S})}\Bigr{]}$
$\displaystyle\>\;+\lambda_{e}S_{T}\sqrt{1-\epsilon^{2}}\cos(\phi_{h}\\!-\\!\phi_{S})W_{LT}^{\cos(\phi_{h}\\!-\phi_{S})}\Bigr{\\}}\,,$
(7)

where $\sigma_{0}\equiv\alpha^{2}\pi y\kappa_{\gamma}/[zQ^{2}(1-\epsilon)]$,
$\alpha$ is the fine-structure constant, $\lambda_{e}$ is the lepton beam
helicity, $S^{\mu}=(0,S_{T}\cos\phi_{S},S_{T}\sin\phi_{S},-S_{L})$ is the
nucleon spin vector decomposed in the Trento frame Bacchetta  _et al._ (2007),
$\epsilon=(1-y)/(1-y+y^{2}/2)$, and $\kappa_{\gamma}=1$ for $M\ll Q$. We have
only kept the structure functions $W$ that have nonvanishing contributions at
leading power in $\lambda$. They can be written in terms of the hard function
$\mathcal{H}(Q^{2})$ and various $b_{T}$-space TMD PDFs $\tilde{g}(x,b_{T})$
and TMD FFs $\tilde{D}(z,b_{T})$ Boer  _et al._ (2011):

$\displaystyle W_{PP^{\prime}}^{\rm ang}$
$\displaystyle\propto\mathcal{F}\bigl{[}\mathcal{H}\,\tilde{g}^{(n)}\\!\tilde{D}^{(m)}\bigr{]}$
(8) $\displaystyle\equiv
2z\int_{0}^{\infty}\\!\frac{\mathrm{d}b_{T}b_{T}}{2\pi}\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{g}^{(n)}\tilde{D}^{(m)}\bigr{]}J_{n+m}(b_{T}P_{hT}/z)$

with
$\mathcal{I}[\mathcal{H}\,\tilde{g}^{(n)}\tilde{D}^{(m)}]\equiv(Mb_{T})^{n}(-M_{h}b_{T})^{m}\sum_{f}\\!\mathcal{H}_{f}\tilde{g}_{f}^{(n)}\tilde{D}_{f}^{(m)}$,
where $f$ sums over quarks and antiquarks. For example,
$W_{UU}^{\cos(2\phi_{h})}=-\mathcal{F}\bigl{[}\mathcal{H}\,\tilde{h}_{1}^{\perp(1)}\tilde{H}_{1}^{\perp(1)}\bigr{]}$,
where $\tilde{h}_{1}^{\perp(1)}$ and $\tilde{H}_{1}^{\perp(1)}$ are the
standard Boer-Mulders Boer and Mulders (1998) and Collins Collins (1993)
functions. We refer the reader to Ref. Ebert  _et al._ (2022) for details on
our notation and conventions.

To compute the spectrum differential in $x,y,z$ and $q_{*}$, we insert the
leading-power measurement $\delta(q_{*}+\sin\phi_{h}P_{hT}/z)$ and
analytically perform the integral over
$\mathrm{d}^{2}\vec{P}_{hT}=\mathrm{d}P_{hT}\,P_{hT}\mathrm{d}\phi_{h}$. As an
explicit example, we work out the contribution from
$W_{UU}^{\cos(2\phi_{h})}$. Using Eq. (8):

$\displaystyle\int_{0}^{\infty}\\!\mathrm{d}P_{hT}\,P_{hT}\int_{0}^{2\pi}\\!\mathrm{d}\phi_{h}\,\delta\bigl{(}q_{*}+\sin\phi_{h}P_{hT}/z\bigr{)}\cos(2\phi_{h})$
$\displaystyle\times(-2z)\int\\!\frac{\mathrm{d}b_{T}b_{T}}{2\pi}\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{h}_{1}^{\perp(1)}\\!\tilde{H}_{1}^{\perp(1)}\bigr{]}\,J_{2}(b_{T}P_{hT}/z)$
$\displaystyle=$
$\displaystyle-\frac{2z^{3}}{\pi}\int\\!\mathrm{d}b_{T}\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{h}_{1}^{\perp(1)}\\!\tilde{H}_{1}^{\perp(1)}\bigr{]}$
$\displaystyle\times\int_{0}^{2\pi}\\!\\!\\!\frac{\mathrm{d}\phi_{H}}{\sin^{2}\\!\phi_{h}}\,\Theta\Bigl{(}-\frac{q_{*}}{\sin\phi_{h}}\Bigr{)}\cos(2\phi_{h})\,\frac{b_{T}\lvert
q_{*}\rvert}{2}\,J_{2}\Bigl{(}\frac{b_{T}q_{*}}{\sin\phi_{h}}\Bigr{)}$
$\displaystyle=$
$\displaystyle-\frac{2z^{3}}{\pi}\int\mathrm{d}b_{T}\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{h}_{1}^{\perp(1)}\\!\tilde{H}_{1}^{\perp(1)}\bigr{]}\,\cos(q_{*}b_{T})\,.$
(9)

The $\phi_{h}$ integral, which is specific to the structure function, can only
depend on $q_{*}b_{T}$ by dimensional analysis, and in this case yields a
simple $\cos(q_{*}b_{T})$.

Carrying out calculations for all contributions, the leading-power SIDIS cross
section differential in $q_{*}$ is:

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*}}=\frac{2z^{3}}{\pi}\sigma_{0}\int_{0}^{\infty}\\!\\!\\!\mathrm{d}b_{T}\,\Bigl{\\{}\cos(q_{*}b_{T})\Bigl{(}\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{f}_{1}\,\tilde{D}_{1}\bigr{]}$
$\displaystyle\quad-\epsilon\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{h}_{1}^{\perp(1)}\\!\tilde{H}_{1}^{\perp(1)}\bigr{]}+\lambda_{e}S_{L}\sqrt{1\\!-\\!\epsilon^{2}}\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{g}_{1L}\,\tilde{D}_{1}\bigr{]}\Bigr{)}$
$\displaystyle\>\;+\cos\phi_{S}\sin(q_{*}b_{T})\,S_{T}\Bigl{(}\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{f}_{1T}^{\perp(1)}\tilde{D}_{1}\bigr{]}+\epsilon\,\mathcal{I}[\mathcal{H}\,\tilde{h}_{1}\,\tilde{H}_{1}^{\perp(1)}]$
$\displaystyle\quad+\frac{\epsilon}{4}\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{h}_{1T}^{\perp(2)}\tilde{H}_{1}^{\perp(1)}\bigr{]}\Bigr{)}$
$\displaystyle\>\;-\sin\phi_{S}\sin(q_{*}b_{T})\,\lambda_{e}S_{T}\sqrt{1-\epsilon^{2}}\,\mathcal{I}\bigl{[}\mathcal{H}\,\tilde{g}_{1T}^{\perp(1)}\tilde{D}_{1}\bigr{]}\Bigr{\\}}\,.$
(10)

We stress that the TMD PDFs and FFs
$\tilde{f}_{1},\tilde{D}_{1},\tilde{h}_{1}^{\perp(1)},\dots$ are the _same_
functions as in the standard factorization for the $P_{hT}$ spectrum and TMD
spin correlations. This is analogous to the role of the unpolarized Banfi  _et
al._ (2009, 2011b) and Boer-Mulders Ebert  _et al._ (2021) TMD PDFs in the
Drell-Yan $\phi^{*}_{\eta}$. The factorization theorem can equivalently be
written in terms of momentum-space TMDs, see Supplement IV.4. Crucially,
definite subsets of these TMDs contribute to the even and odd parts of the
spectrum under $q_{*}\to-q_{*}$. The odd parts are accessible through the
asymmetry $\mathrm{d}\sigma(q_{*}>0)-\mathrm{d}\sigma(q_{*}<0)$. Contributions
can be further disentangled experimentally through their unique dependence on
$\lambda_{e}$, $S^{\mu}$ and $\epsilon$, i.e., by taking asymmetries with
opposite beam polarizations and by measuring cross sections as a function of
$y$.111We recommend reconstructing $S^{\mu}$ using a rotation by $\theta_{h}$
to maintain a purely angular measurement, which is justified at LP as in Eq.
(II), see Supplement IV.3. Note that the transversity and pretzelosity PDFs
have a degenerate contribution $\epsilon S_{T}(h_{1}+h_{1T}^{\perp}/4)$ to
$q_{*}$ in Eq. (III), while the worm-gear $L$ function $h_{1L}^{\perp}$ drops
out. We leave the creation of other purely angular observables that retain
sensitivity to these TMDs to future work. E.g., the double asymmetry for
$q_{*}\to-q_{*}$ and $\lambda_{e}\to-\lambda_{e}$ as a function of $x$ and
$\lvert q_{*}\rvert$ gives direct access to the worm-gear $T$ function
$\tilde{g}_{1T}^{\perp}(x,b_{T})$.

##  IV Experimental sensitivity

Figure 2: Expected event-level detector resolution $\sigma_{\mathcal{O}}$ for
different SIDIS TMD observables $\mathcal{O}$. We show relative resolutions as
a function of the magnitude of $\mathcal{O}=q_{*}$ (solid red), $P_{hT}/z$
(dashed blue), and $\sin\phi_{h}P_{hT}/z$ (dotted green).

To show the improvement that $q_{*}$ makes for TMD analyses, we use Pythia
8.306 Bierlich  _et al._ (2022) to simulate $e^{-}p\to e^{-}X$ at $18\times
275\,\mathrm{GeV}^{2}$, disabling QED corrections. We select on $h=\pi^{+}$
and apply the following cuts (“SIDIS cuts”):

$\displaystyle x$ $\displaystyle>0.001\,,\quad$ $\displaystyle 0.01$
$\displaystyle<y<0.95\,,\quad z>0.05\,,$ $\displaystyle Q^{2}$
$\displaystyle>16\,\mathrm{GeV}^{2}\,,\;$ $\displaystyle W^{2}$
$\displaystyle=(P+q)^{2}>100\,\mathrm{GeV}^{2}\,.$ (11)

Scaled to an integrated EIC luminosity of $10\,\mathrm{fb}^{-1}$, this results
in a sample with $N_{\pi^{+}}=4.18\times 10^{8}$. We first assess the expected
event-level detector resolution of $q_{*}$ compared to $P_{hT}$. To model the
detector resolution, we apply Gaussian smearing to the final-state electron
and hadron momenta, assuming a tracking detector that matches the performance
given in Ref. Abdul Khalek  _et al._ (2022): a resolution of
$\sigma_{p}/p=0.05\%\,p/\\!\,\mathrm{GeV}\oplus 0.5\%$ on the momentum
$p=\lvert\vec{p\mkern-2.0mu}\mkern 2.0mu\rvert$ of charged particles in the
central barrel region $\lvert\eta\rvert<1$, $0.05\%\,p/\\!\,\mathrm{GeV}\oplus
1\%$ in the inner endcap $1<\lvert\eta\rvert<2.5$, and
$0.1\%\,p/\\!\,\mathrm{GeV}\oplus 2\%$ in the outer endcap
$2.5<\lvert\eta\rvert<3.5$. A particle is forward (backward) if it has
$1<|\eta|<3.5$ and $\eta>0$ ($\eta<0$). We assume a fixed angular resolution
of $\sigma_{\theta}=\sigma_{\phi}=0.001$. We ignore the electromagnetic
calorimeter as its $e^{-}$ energy resolution is expected to be a factor of two
worse than the tracking resolution for the relevant energies and central or
inner-endcap detector Abdul Khalek  _et al._ (2022). Our results for the
detector resolution on $q_{*}$ compared to $P_{hT}$ are shown in Fig. 2 for
the case of a backward $e^{-}$ and a central $h$, which accounts for the
largest share ($41\%$) of the event sample. We see that $q_{*}$ improves over
the resolution of $P_{hT}/z$ by an order of magnitude across the full strongly
confined TMD region $q_{*},P_{hT}/z\lesssim 2\,\mathrm{GeV}$. Similar results
are obtained for other cases, see Supplement IV.5. It is also interesting to
compare $q_{*}$ to a direct measurement of $\sin\phi_{h}P_{hT}/z$, which it
reduces to at leading power, and which has the same factorized cross section
in Eq. (III). It also has improved resolution over $P_{hT}/z$ at small values
thanks to picking up on the same acoplanarity of the event, which is stable
against the electron momentum resolution, but it does not exhibit the same
improved resolution as $q_{*}$ since it is not a pure angular variable.

Figure 3: Statistical sensitivity to TMD nonperturbative model coefficients at
the $10\,\mathrm{fb}^{-1}$ EIC when measuring $P_{hT}$ (blue) or $q_{*}$ (red)
compared to the prior (MAPTMD22 fit Bacchetta  _et al._ (2022), gray). Despite
its superior resolution, $q_{*}$ enjoys comparable statistical sensitivity to
$P_{hT}$.

To verify the statistical sensitivity of $q_{*}$ to TMD physics we perform a
Bayesian reweighting analysis of the unpolarized cross section
$\mathrm{d}\sigma/(\mathrm{d}x\,\mathrm{d}z\,\mathrm{d}Q^{2}\,\mathrm{d}\mathcal{O})$
for $\mathcal{O}=P_{hT}/z$ and $\lvert q_{*}\rvert$. We assume that the
$\mathcal{O}$ spectrum is measured in twenty equally spaced bins between
$0\leq\mathcal{O}\leq 4\,\mathrm{GeV}$ inside 1000 three-dimensional bins in
$x,z,Q^{2}$, with equal total statistics $N_{\pi^{+}}/1000$ in each, and for
definiteness consider a bin centered on $x=0.1$, $z=0.15$,
$Q=20\,\mathrm{GeV}$ in the following. To account for the fact that the
factorized dependence on $x$ ($z$, $Q^{2}$) is determined from all bins at
equal $x$ ($z$, $Q^{2}$) simultaneously, we multiply the available statistics
by another factor $100$, arriving at an effective sample size of
$N_{\mathrm{eff}}=N_{\pi^{+}}/10=4.18\times 10^{7}$. At fixed $x,z,Q^{2}$, one
common model for the nonperturbative TMD physics is Bacchetta  _et al._ (2022)

$\displaystyle\tilde{f}_{1}^{\,\mathrm{NP}}$
$\displaystyle=e^{-\omega_{1}b_{T}^{2}}\,,$
$\displaystyle\tilde{D}_{1}^{\mathrm{NP}}$
$\displaystyle=\alpha\,e^{-\omega_{2}b_{T}^{2}}+(1-\alpha)(1-\omega_{3}b_{T}^{2})\,e^{-\omega_{3}b_{T}^{2}}\,,$
(12)

where $\omega_{i}$ encode the Gaussian or weighted-Gaussian shape of the TMDs
at large distances. We are interested in how much better the three free
parameters $\omega_{i}$, can be determined at the EIC using either $q_{*}$ or
$P_{hT}$. (We hold the relative size $\alpha$ of the two contributions fixed
for simplicity.) We take a Gaussian prior probability density
$\pi(\omega_{i})$ based on the central values and standard deviations from
Bacchetta  _et al._ (2022). We combine Eq. (IV) with leading-logarithmic TMD
evolution and tree-level matching in SCETlib Ebert  _et al._ (2018), and
insert it into Eqs. (III) and (III) to generate EIC pseudodata $d_{n}$ for bin
$n$ in $\mathcal{O}$, using the central $\omega_{i}$. In the same way, we
generate theory replicas $t_{n}(\omega_{i})$ distributed according to the
prior. By normalizing $\sum d_{n}=\sum t_{n}=1$ to the sum over bins at fixed
$x$ and $z$, collinear PDFs and FFs drop out at this order. (Details on the
theory calculation are given in Supplement IV.6.) We then sample the posterior
parameter probability distribution

$\displaystyle\pi(\omega_{i}\,|\,d_{n})\propto\exp\biggl{[}-\sum_{n}\Bigl{(}\frac{d_{n}-t_{n}(\omega_{i})}{\sigma_{n}}\Bigr{)}^{2}\biggr{]}\,\pi(\omega_{i})$
(13)

using a standard $\chi^{2}$ likelihood function, where
$\sigma_{n}=\sqrt{d_{n}/N_{\mathrm{eff}}}$ is the Poisson error on pseudodata
bin $n$. Our results for the mean and variance of the posterior distribution
are shown in Fig. 3 compared to those of the prior. Comparing $P_{hT}$ and
$q_{*}$, we find that the superior experimental resolution of $q_{*}$ only
requires giving up a minor amount of statistical sensitivity to the
$\omega_{i}$. In particular, there is more than a factor 10 improvement in
uncertainty on the dominant fragmentation parameter $\omega_{3}$ in either
case. The choice of binning at small $q_{*}$ should be optimized in the future
to exploit its excellent resolution, but we stress that we have not done so
here.

Figure 4: Strength of bias on TMD nonperturbative model coefficients from
potential sources of systematic uncertainty. By construction $q_{*}$ (red) is
robust against the large calibration uncertainty $\delta_{p_{e}}$ that impacts
$P_{hT}$ (blue). Both $q_{*}$ and $P_{hT}$ exhibit similar susceptibility to
non-uniform detector response, modeled by $\Delta\epsilon_{X}$.

The same setup can be used to assess the robustness of $q_{*}$ against
systematic uncertainties. We use Pythia Bierlich  _et al._ (2022) to generate
biased pseudodata $d_{n}^{\mathrm{bias}}$ subject to either (i) a momentum
miscalibration, $p\to(1+\delta_{p})\,p$, or (ii) a shape effect from a non-
uniform detector response (encoded e.g. in an efficiency) that changes at a
slow rate $\Delta\epsilon_{X}$ as a function of
$X=\\{p_{e},p_{h},\eta_{e},\eta_{h}\\}$ across the $x,z,Q^{2}$ bin at hand.
(The absolute value of the efficiency cancels out in the normalized
$d_{n}^{\mathrm{bias}}$.) Repeating the reweighting analysis, we evaluate the
partial derivatives of the posterior’s mean $\omega_{i}$ with respect to the
bias parameters, which we dub the “strength” of the bias, as shown in Fig. 4.
As anticipated, we find that an analysis using $P_{hT}$ is severely
susceptible to the electron momentum calibration $\delta_{p_{e}}$, while the
calibration uncertainty using $q_{*}$ vanishes exactly due to its purely
angular nature. (Both $P_{hT}/z$ and $q_{*}$ are independent of
$\delta_{p_{h}}$ at observable level for $M_{h}\ll Q$.) Figure 4 also shows
that $P_{hT}$ and $q_{*}$ have roughly comparable susceptibility to non-
uniform detector response, despite the exponential factors of $\eta_{e,h}$
appearing in $q_{*}$, demonstrating the robustness of $q_{*}$ against these
sources of bias.

By replacing measurements of
$\mathrm{d}\sigma/\mathrm{d}P_{hT}\mathrm{d}\phi_{h}$ by
$\mathrm{d}\sigma/\mathrm{d}q_{*}$ with angular reconstruction of $Q,x,y$, the
prospects for precisely mapping the 3D structure of hadronization and
confinement with TMDs are bright. We anticipate a follow-up campaign to aid
this endeavor by discovering other useful angular observables and by studying
theoretical ingredients, like the convergence of the known higher-order QCD
corrections for these cross sections. Many of the improvements provided by
$q_{*}$ will also carry over to the case where the identified hadron is
replaced by a jet, in particular if the jet $\phi_{\rm acop}$ is determined by
the sum of the signed acoplanarity angles of the particles inside the jet to
preserve the angular nature of the measurement.

###### Acknowledgements.

##### Acknowledgments

This work was supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, from DE-SC0011090. I.S. was also supported in part
by the Simons Foundation through the Investigator grant 327942.

## References

  * Collins (2011) J. Collins,  _Foundations of perturbative QCD_, Cambridge monographs on particle physics, nuclear physics, and cosmology (Cambridge Univ. Press, New York, NY, 2011).
  * Airapetian  _et al._ (2013) A. Airapetian _et al._ (HERMES), Phys. Rev. D87, 074029 (2013), arXiv:1212.5407 [hep-ex] .
  * Airapetian  _et al._ (2019) A. Airapetian _et al._ (HERMES), Phys. Lett. B 797, 134886 (2019), arXiv:1903.08544 [hep-ex] .
  * Airapetian  _et al._ (2020) A. Airapetian _et al._ (HERMES), JHEP 12, 010 (2020), arXiv:2007.07755 [hep-ex] .
  * Alekseev  _et al._ (2009) M. Alekseev _et al._ (COMPASS), Phys. Lett. B673, 127 (2009), arXiv:0802.2160 [hep-ex] .
  * Aghasyan  _et al._ (2018) M. Aghasyan _et al._ (COMPASS), Phys. Rev. D97, 032006 (2018), arXiv:1709.07374 [hep-ex] .
  * Parsamyan (2018) B. Parsamyan, PoS DIS2017, 259 (2018), arXiv:1801.01488 [hep-ex] .
  * Aschenauer  _et al._ (2015) E.-C. Aschenauer, A. Bazilevsky, M. Diehl, J. Drachenberg, K. O. Eyser, _et al._ , (2015), arXiv:1501.01220 [nucl-ex] .
  * Adamczyk  _et al._ (2016) L. Adamczyk _et al._ (STAR), Phys. Rev. Lett. 116, 132301 (2016), arXiv:1511.06003 [nucl-ex] .
  * Avakian  _et al._ (2004) H. Avakian _et al._ (CLAS), Phys. Rev. D 69, 112004 (2004), arXiv:hep-ex/0301005 .
  * Jawalkar  _et al._ (2018) S. Jawalkar _et al._ (CLAS), Phys. Lett. B 782, 662 (2018), arXiv:1709.10054 [nucl-ex] .
  * Moran  _et al._ (2022) S. Moran _et al._ (CLAS), Phys. Rev. C 105, 015201 (2022), arXiv:2109.09951 [nucl-ex] .
  * Scimemi and Vladimirov (2020) I. Scimemi and A. Vladimirov, JHEP 06, 137 (2020), arXiv:1912.06532 [hep-ph] .
  * Bacchetta  _et al._ (2020) A. Bacchetta, V. Bertone, C. Bissolotti, G. Bozzi, F. Delcarro, F. Piacenza, and M. Radici, JHEP 07, 117 (2020), arXiv:1912.07550 [hep-ph] .
  * Bury  _et al._ (2021) M. Bury, A. Prokudin, and A. Vladimirov, JHEP 05, 151 (2021), arXiv:2103.03270 [hep-ph] .
  * Abdul Khalek  _et al._ (2022) R. Abdul Khalek _et al._ , Nucl. Phys. A 1026, 122447 (2022), arXiv:2103.05419 [physics.ins-det] .
  * Ebert  _et al._ (2019) M. A. Ebert, I. W. Stewart, and Y. Zhao, Phys. Rev. D99, 034505 (2019), arXiv:1811.00026 [hep-ph] .
  * Shanahan  _et al._ (2020) P. Shanahan, M. Wagman, and Y. Zhao, Phys. Rev. D 102, 014511 (2020), arXiv:2003.06063 [hep-lat] .
  * Shanahan  _et al._ (2021) P. Shanahan, M. Wagman, and Y. Zhao, Phys. Rev. D 104, 114502 (2021), arXiv:2107.11930 [hep-lat] .
  * Schlemmer  _et al._ (2021) M. Schlemmer, A. Vladimirov, C. Zimmermann, M. Engelhardt, and A. Schäfer, JHEP 08, 004 (2021), arXiv:2103.16991 [hep-lat] .
  * Li  _et al._ (2022) Y. Li _et al._ , Phys. Rev. Lett. 128, 062002 (2022), arXiv:2106.13027 [hep-lat] .
  * Chu  _et al._ (2022) M.-H. Chu _et al._ (LPC), (2022), arXiv:2204.00200 [hep-lat] .
  * Banfi  _et al._ (2011a) A. Banfi, S. Redford, M. Vesterinen, P. Waller, and T. R. Wyatt, Eur. Phys. J. C 71, 1600 (2011a), arXiv:1009.1580 [hep-ex] .
  * Abazov  _et al._ (2011) V. M. Abazov _et al._ (D0), Phys. Rev. Lett. 106, 122001 (2011), arXiv:1010.0262 [hep-ex] .
  * Aaij  _et al._ (2013) R. Aaij _et al._ (LHCb), JHEP 02, 106 (2013), arXiv:1212.4620 [hep-ex] .
  * Aad  _et al._ (2016) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 76, 291 (2016), arXiv:1512.02192 [hep-ex] .
  * Aaij  _et al._ (2015) R. Aaij _et al._ (LHCb), JHEP 08, 039 (2015), arXiv:1505.07024 [hep-ex] .
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 80, 616 (2020), arXiv:1912.02844 [hep-ex] .
  * Bacchetta  _et al._ (2004) A. Bacchetta, U. D’Alesio, M. Diehl, and C. A. Miller, Phys. Rev. D 70, 117504 (2004), arXiv:hep-ph/0410050 .
  * Bentvelsen  _et al._ (1992) S. Bentvelsen, J. Engelen, and P. Kooijman, in  _Workshop on Physics at HERA_ (1992).
  * Ji  _et al._ (2005) X.-d. Ji, J.-p. Ma, and F. Yuan, Phys. Rev. D 71, 034005 (2005), arXiv:hep-ph/0404183 .
  * Ji  _et al._ (2004) X.-d. Ji, J.-P. Ma, and F. Yuan, Phys. Lett. B 597, 299 (2004), arXiv:hep-ph/0405085 .
  * Bacchetta  _et al._ (2007) A. Bacchetta, M. Diehl, K. Goeke, A. Metz, P. J. Mulders, and M. Schlegel, JHEP 02, 093 (2007), arXiv:hep-ph/0611265 .
  * Boer  _et al._ (2011) D. Boer, L. Gamberg, B. Musch, and A. Prokudin, JHEP 10, 021 (2011), arXiv:1107.5294 [hep-ph] .
  * Boer and Mulders (1998) D. Boer and P. J. Mulders, Phys. Rev. D 57, 5780 (1998), arXiv:hep-ph/9711485 .
  * Collins (1993) J. C. Collins, Nucl. Phys. B 396, 161 (1993), arXiv:hep-ph/9208213 .
  * Ebert  _et al._ (2022) M. A. Ebert, A. Gao, and I. W. Stewart, JHEP 06, 007 (2022), arXiv:2112.07680 [hep-ph] .
  * Banfi  _et al._ (2009) A. Banfi, M. Dasgupta, and R. M. Duran Delgado, JHEP 12, 022 (2009), arXiv:0909.5327 [hep-ph] .
  * Banfi  _et al._ (2011b) A. Banfi, M. Dasgupta, and S. Marzani, Phys. Lett. B701, 75 (2011b), arXiv:1102.3594 [hep-ph] .
  * Ebert  _et al._ (2021) M. A. Ebert, J. K. L. Michel, I. W. Stewart, and F. J. Tackmann, JHEP 04, 102 (2021), arXiv:2006.11382 [hep-ph] .
  * Bierlich  _et al._ (2022) C. Bierlich _et al._ , (2022), arXiv:2203.11601 [hep-ph] .
  * Bacchetta  _et al._ (2022) A. Bacchetta, V. Bertone, C. Bissolotti, G. Bozzi, M. Cerutti, F. Piacenza, M. Radici, and A. Signori, (2022), arXiv:2206.07598 [hep-ph] .
  * Ebert  _et al._ (2018) M. A. Ebert, J. K. L. Michel, F. J. Tackmann, _et al._ , DESY-17-099 (2018), webpage: http://scetlib.desy.de.

## Supplemental material

###  IV.1 Constructing $q_{*}$ with finite target mass

In the main text, we present our results in the light target mass limit $M\ll
Q$. Here, we give the corresponding results when fully retaining the target
mass $M$, which can be important when $Q$ is not sufficiently large or when
$M$ is large (such as for a nucleus). This amounts to including in our
construction the dependence on

$\displaystyle\gamma=\frac{2xM}{Q}\,.$ (S1)

We also have the following variable generalizations:

$\displaystyle\epsilon=\frac{1-y-\frac{1}{4}y^{2}\gamma^{2}}{1-y+\frac{1}{4}y^{2}\gamma^{2}+\frac{1}{2}y^{2}}\,,\qquad\kappa_{\gamma}=\Bigl{[}1-\Bigl{(}\frac{m_{hT}\,\gamma}{zQ}\Bigr{)}^{2}\Bigr{]}^{-1/2}\,,\qquad
m_{hT}^{2}=M_{h}^{2}+P_{hT}^{2}\,.$ (S2)

Note that we still take $P_{hT},M_{h}\ll Q$ (appropriate for example when $N$
is a proton or ion and $h$ is a pion), and hence can still take
$\kappa_{\gamma}=1$. Only $\epsilon$ receives mass corrections through
$\gamma$. We will show that the above corrections do not change any of the
conclusions regarding the utility of $q_{*}$.

We continue to work with the leading power (LP) kinematics $\lambda\sim
P_{hT}/(zQ)\ll 1$. For the relationship between the acoplanarity angle
$\phi_{\mathrm{acop}}^{\mathrm{rest}}$ and the hadron transverse momentum
$P_{hT}$, we now have

$\displaystyle\tan\phi_{\mathrm{acop}}^{\mathrm{rest}}$
$\displaystyle=\frac{\sin\phi_{h}\,P_{hT}}{zQ}\sqrt{\frac{1+4M^{2}x^{2}/Q^{2}}{1-M^{2}x^{2}y^{2}/Q^{2}-y}}+\mathcal{O}(\lambda^{2})$
$\displaystyle=\frac{\sin\phi_{h}\,P_{hT}}{zQ}\sqrt{\frac{1+\gamma^{2}}{1-\gamma^{2}y^{2}/4-y}}+\mathcal{O}(\lambda^{2})\,.$
(S3)

We now wish to construct an optimized observable $q^{M}_{*}$ for a massive
target with $M\sim Q$, such that $q^{M}_{*}\stackrel{{\scriptstyle M\ll
Q}}{{=}}q_{*}$ while retaining the desired leading power relation
$q^{M}_{*}\stackrel{{\scriptstyle\text{LP}}}{{=}}-\sin\phi_{h}\frac{P_{hT}}{z}$.

First of all, we emphasize that the following relations presented in the main
text are independent of $M$:

$\displaystyle y=1-\frac{\sin\theta_{h}}{\cos\alpha}\,,\qquad
Q^{2}=(\ell^{0}_{\mathrm{rest}})^{2}\Bigl{[}\frac{\sin^{2}\theta_{e}}{\cos^{2}\alpha}-\bigl{(}1-\frac{\sin\theta_{h}}{\cos\alpha}\bigr{)}^{2}\Bigr{]}\,.$
(S4)

Furthermore, we notice that $\gamma^{2}$ can also be written in terms of
angles in a manner that is independent of $M$:

$\displaystyle\gamma^{2}=\frac{4M^{2}x^{2}}{Q^{2}}=\frac{M^{2}Q^{2}}{(P\cdot
q)^{2}}=\Bigl{(}\frac{\sin\theta_{e}}{\cos\alpha-\sin\theta_{h}}\Bigr{)}^{2}-1\,.$
(S5)

This allows us to define $q^{M}_{*}$ in terms of _target rest frame_
quantities completely free of $M$ dependence:

$\displaystyle
q^{M}_{*}\equiv-\ell^{0}_{\mathrm{rest}}\frac{[\cos(\theta_{e}+\theta_{h})+\cos\theta_{h}]\tan(\frac{\theta_{e}}{2})\sin\theta_{h}}{\sin(\theta_{e}+\theta_{h})}\,\tan\phi_{\mathrm{acop}}^{\mathrm{rest}}\stackrel{{\scriptstyle\text{LP}}}{{=}}-\sin\phi_{h}\frac{P_{hT}}{z}\,.$
(S6)

This relation is especially useful for fixed target experiments, where all the
quantities can be readily measured. For collider experiments like the EIC, the
above equation may be expressed in terms of lab frame quantities simply by
substituting
$\tan\frac{\theta_{i}}{2}=\exp\bigl{[}\cosh^{-1}(P^{0}_{\mathrm{EIC}}/M)+\eta_{i}\bigr{]}$
for $i=e,\,h$, where $\eta_{i}$ are the lab frame pseudorapidities,
$\ell^{0}_{\mathrm{rest}}=\exp\bigl{[}\cosh^{-1}(P^{0}_{\mathrm{EIC}}/M)\bigr{]}\ell^{0}_{\mathrm{EIC}}$,
and
$\phi_{\mathrm{acop}}^{\mathrm{rest}}=-\phi_{\mathrm{acop}}^{{\mathrm{EIC}}}$.
We remind the reader that the last relation is due to the different convention
for the orientation of the $z$ axis at the EIC, see Fig. 1.

We emphasize that since $q^{M}_{*}$ has the desired LP limit
$q^{M}_{*}\stackrel{{\scriptstyle\text{LP}}}{{=}}-\sin\phi_{h}\frac{P_{hT}}{z}$,
the factorization formula for
$\mathrm{d}\sigma/(\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*}^{M})$
is the same as written in Eq. (III), does not receive mass corrections, and is
valid for both fixed-target and collider experiments.

###  IV.2 Power corrections to double angle formulas for $Q^{2}$, $y$, and $x$

In the main text, we give the expression of $Q^{2}$ and $y$ in Eq. (II) in
terms of lab frame angles to leading order in $\lambda\sim P_{hT}/(zQ)$. Here,
we derive the leading power correction to these kinematic relations. These
results can be used to get an idea of the size of power corrections to an
analysis, including both i) the size of corrections to the double angle
construction for $x$, $y$, and $Q^{2}$, and ii) power corrections to the
factorization formula for
$\mathrm{d}\sigma/(\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*})$.
Note that in this section we still work in the massless target limit $M\ll Q$.

We have the following relations:

$\displaystyle s$ $\displaystyle=(P+\ell)^{2}=2P\cdot\ell=\frac{Q^{2}}{xy}\,,$
(S7a) $\displaystyle\frac{(2P^{0}_{\mathrm{EIC}})^{2}}{s}\,\exp(-2\eta_{e})$
$\displaystyle=\frac{P\cdot\ell^{\prime}}{\ell\cdot\ell^{\prime}}=\frac{P\cdot
q-P\cdot\ell}{\ell\cdot q}=\frac{1}{x}\frac{1-y}{y}\,,$ (S7b)
$\displaystyle\frac{(2P^{0}_{\mathrm{EIC}})^{2}}{s}\exp(-2\eta_{h})$
$\displaystyle=\frac{P\cdot P_{h}}{\ell\cdot P_{h}}=\frac{P\cdot
q}{x(\ell\cdot P)+\ell\cdot q+\frac{1}{z}\,\ell\cdot P_{hT}}$
$\displaystyle=\frac{y}{x(1-y)}\Bigl{(}1+\frac{2\cos\phi_{h}}{\sqrt{1-y}}\,\frac{P_{hT}}{zQ}\Bigr{)}+\mathcal{O}(\lambda^{2})\,,$
(S7c)

where $\eta_{e}$ and $\eta_{h}$ are lab frame pseudorapidities. Note that we
have indicated that Eq. (S7c) is the only equation here that receives
corrections in $\lambda\sim P_{hT}/(zQ)$ when expanding in this ratio. Solving
the above relations for $x$, $y$, and $Q^{2}$, we get

$\displaystyle x$
$\displaystyle=\frac{(2P^{0}_{\mathrm{EIC}})^{2}}{s}e^{\eta_{e}+\eta_{h}}\,\biggl{(}1+\cos\phi_{h}\sqrt{1+e^{-\Delta\eta}}\,\frac{P_{hT}}{zQ}\biggr{)}+\mathcal{O}(\lambda^{2}),$
(S8a) $\displaystyle y$
$\displaystyle=\frac{1}{1+e^{\Delta\eta}}+\mathcal{O}(\lambda^{2}),$ (S8b)
$\displaystyle Q^{2}$
$\displaystyle=(2P^{0}_{\mathrm{EIC}})^{2}\frac{e^{\eta_{e}+\eta_{h}}}{1+e^{\Delta\eta}}\,\biggl{(}1+\cos\phi_{h}\sqrt{1+e^{-\Delta\eta}}\,\frac{P_{hT}}{zQ}\biggr{)}+\mathcal{O}(\lambda^{2})\,,$
(S8c)

where $\Delta\eta=\eta_{h}-\eta_{e}$ is boost invariant along the
$z$-direction. Notice that in Eq. (S8b), $y$ does not receive linear
corrections in $\lambda$. An application of Eq. (S8) is to test the size of
the power corrections in the expressions for $x$ and $Q^{2}$ in a given data
set, and thus apply cuts to restrict the data to TMD regions where the leading
term is dominant.

We can also invert the formulae in Eq. (S8) to define a set of variables
$x_{*},\,y_{*}$ and $Q_{*}$ that use lab frame measurements. The variables
$x_{*},\,y_{*}$ and $Q_{*}$ agree with the kinematic invariants $x,\,y$ and
$Q$ up to the determined ${\cal O}(\lambda)$ power corrections:

$\displaystyle x_{*}$
$\displaystyle\equiv\frac{\bigl{(}2P_{\mathrm{EIC}}^{0}\bigr{)}^{2}}{s}e^{\eta_{e}+\eta_{h}}\,=x\biggl{(}1-\cos\phi_{h}\sqrt{\frac{1}{1-y}}\,\frac{P_{hT}}{zQ}\biggr{)}+\mathcal{O}(\lambda^{2}),$
(S9a) $\displaystyle y_{*}$
$\displaystyle\equiv\frac{1}{1+e^{\Delta\eta}}=y+\mathcal{O}(\lambda^{2}),$
(S9b) $\displaystyle Q_{*}^{2}$
$\displaystyle\equiv\bigl{(}2P_{\mathrm{EIC}}^{0}\bigr{)}^{2}\frac{e^{\eta_{e}+\eta_{h}}}{1+e^{\Delta\eta}}\,=Q^{2}\biggl{(}1-\cos\phi_{h}\sqrt{\frac{1}{1-y}}\,\frac{P_{hT}}{zQ}\biggr{)}+\mathcal{O}(\lambda^{2})\,.$
(S9c)

The all-order definition of $q_{*}$ in the main text, Eq. (4), can be written
as $q_{*}\equiv Q_{*}\sqrt{1-y_{*}}\tan\phi_{\mathrm{acop}}^{\mathrm{EIC}}$.
This allows us to easily determine the leading power correction to the formula
for $q_{*}$ obtained by expanding in $P_{hT}\ll zQ$:

$\displaystyle
q_{*}=-\sin\phi_{h}\frac{P_{hT}}{z}\biggl{(}1-\frac{\cos\phi_{h}}{2}\sqrt{\frac{1}{1-y}}\,\frac{P_{hT}}{zQ}\biggr{)}+\mathcal{O}(\lambda^{3})\,.$
(S10)

This kinematic correction to the relationship between variables is the only
power correction that would give non-trivial dependence on $y$ and
$\cos\phi_{h}$ to the factorization formula in Eq. (III). Hence it can be
unambiguously included in the factorization analysis by using this more
complicated relationship in the $\delta(q_{*}+\ldots)$ when switching
variables and integrating over $P_{hT}$ and $\phi_{h}$ (cf. the example given
in Eq. (III) without these corrections). However, this still does not capture
the dynamic hadronic power corrections, which arise from using the $P_{hT}\ll
zQ$ expansion when deriving the original factorization theorem for
$\mathrm{d}\sigma/(\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}^{2}\vec{P}_{hT})$.

###  IV.3 Leading-power formulas for target spin vector from angular
measurements

As mentioned in the main text, the $S_{L}$, $S_{T}$ and $\phi_{S}$ that appear
in the factorization formula are defined in the Trento frame by writing the
nucleon spin vector as
$S^{\mu}=(0,S_{T}\cos\phi_{S},S_{T}\sin\phi_{S},-S_{L})_{\mathrm{Trento}}$.
Here we give leading-power expressions for these variables in terms of the
target rest frame components of
$S^{\mu}=(0,S_{x},S_{y},S_{z})_{\mathrm{rest}}$, and the EIC lab frame hadron
pseudorapidity $\eta_{h}$. Note we do not assume that the nucleon is in a pure
spin state.

We start with expressing $S_{L}$, $S_{T}$ and $\phi_{S}$ using the polar angle
$\theta_{q}$ of $\vec{q}$,

$\displaystyle S_{L}=-S_{z}\cos\theta_{q}+S_{x}\sin\theta_{q}\,,\qquad
S_{T}=\sqrt{(S_{x}\cos\theta_{q}+S_{z}\sin\theta_{q})^{2}+S_{y}^{2}}\,.$ (S11)

$\phi_{S}$ is given by

$\displaystyle\sin\phi_{S}=\frac{S_{y}}{S_{T}}\,,\qquad\cos\phi_{S}=\frac{S_{x}\cos\theta_{q}+S_{z}\sin\theta_{q}}{S_{T}}\,.$
(S12)

At leading power in $\lambda\sim P_{hT}/(zQ)\ll 1$, we may replace
$\theta_{q}$ by $\theta_{h}$, which can be written in terms of $\eta_{h}$. We
have

$\displaystyle S_{L}$
$\displaystyle=-S_{z}\,\frac{1-A^{2}\,e^{2\eta_{h}}}{1+A^{2}\,e^{2\eta_{h}}}+S_{x}\,\frac{2A\,e^{\eta_{h}}}{1+A^{2}\,e^{2\eta_{h}}}+\mathcal{O}\bigl{(}\lambda\bigr{)}\,,$
(S13) $\displaystyle S_{T}$
$\displaystyle=\sqrt{\left[\frac{(1-A^{2}\,e^{2\eta_{h}})S_{x}+2A\,e^{\eta_{h}}S_{z}}{1+A^{2}\,e^{2\eta_{h}}}\right]^{2}+S_{y}^{2}}\,+\mathcal{O}\bigl{(}\lambda\bigr{)}\,,$
$\displaystyle\sin\phi_{S}$
$\displaystyle=\frac{S_{y}}{S_{T}}\,,\qquad\cos\phi_{S}=\frac{(1-A^{2}\,e^{2\eta_{h}})S_{x}+2A\,e^{\eta_{h}}S_{z}}{(1+A^{2}\,e^{2\eta_{h}})S_{T}}+\mathcal{O}\bigl{(}\lambda\bigr{)}\,.$
(S14)

Here we define $A$ to be

$\displaystyle
A=\exp\bigl{[}\cosh^{-1}(P^{0}_{\mathrm{EIC}}/M)]\stackrel{{\scriptstyle M\ll
P^{0}_{\mathrm{EIC}}}}{{\approx}}\frac{2P^{0}_{\mathrm{EIC}}}{M}\,.$ (S15)

###  IV.4 Momentum space factorization formula for the $q_{*}$ spectrum

In the main, text our factorization theorem for
$\mathrm{d}\sigma/(\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*})$
was written in terms of $b_{T}$ space TMDs. For completeness, here we give the
factorization theorem written in terms of momentum space TMDs, which are
Fourier conjugate to those in $b_{T}$ space.

We start with the standard leading power momentum space TMD factorization
formulae Bacchetta  _et al._ (2007) for the structure functions appearing in
Eq. (III):

$\displaystyle W_{UU,T}$
$\displaystyle=\tilde{\mathcal{F}}\left[\mathcal{H}^{(0)}\,f_{1}D_{1}\right]\,,$
$\displaystyle W_{UU}^{\cos 2\phi_{h}}$
$\displaystyle=\tilde{\mathcal{F}}\left[\frac{-2(\hat{h}\cdot{\vec{k}}_{T})(\hat{h}\cdot{\vec{p}}_{T})+{\vec{k}}_{T}\cdot{\vec{p}}_{T}}{MM_{h}}\mathcal{H}^{(0)}\,h_{1}^{\perp}H_{1}^{\perp}\right]\,,$
$\displaystyle W_{UL}^{\sin 2\phi_{h}}$
$\displaystyle=\tilde{\mathcal{F}}\left[\frac{-2(\hat{h}\cdot{\vec{k}}_{T})(\hat{h}\cdot{\vec{p}}_{T})+{\vec{k}}_{T}\cdot{\vec{p}}_{T}}{MM_{h}}\mathcal{H}^{(0)}\,h_{1L}^{\perp}H_{1}^{\perp}\right]\,,$
$\displaystyle W_{LL}$
$\displaystyle=\tilde{\mathcal{F}}\left[\mathcal{H}^{(0)}\,g_{1L}D_{1}\right]\,,$
$\displaystyle W_{UT,T}^{\sin\left(\phi_{h}-\phi_{S}\right)}$
$\displaystyle=\tilde{\mathcal{F}}\left[-\frac{\hat{h}\cdot{\vec{k}}_{T}}{M}\mathcal{H}^{(0)}\,f_{1T}^{\perp}D_{1}\right]\,,$
$\displaystyle W_{UT}^{\sin\left(\phi_{h}+\phi_{S}\right)}$
$\displaystyle=\tilde{\mathcal{F}}\left[-\frac{\hat{h}\cdot{\vec{p}}_{T}}{M_{h}}\mathcal{H}^{(0)}\,h_{1}H_{1}^{\perp}\right]\,,$
$\displaystyle W_{UT}^{\sin\left(3\phi_{h}-\phi_{S}\right)}$
$\displaystyle=\tilde{\mathcal{F}}\left[\frac{2(\hat{h}\cdot{\vec{k}}_{T})\,(\vec{k}_{T}\cdot\vec{p}_{T})+k_{T}^{2}\,(\hat{h}\cdot{\vec{p}}_{T})-4(\hat{h}\cdot{\vec{k}}_{T})^{2}\,(\hat{h}\cdot{\vec{p}}_{T})}{2M^{2}M_{h}}\,\mathcal{H}^{(0)}h_{1T}^{\perp}H_{1}^{\perp}\right]\,,$
$\displaystyle W_{LT}^{\cos\left(\phi_{h}-\phi_{S}\right)}$
$\displaystyle=\tilde{\mathcal{F}}\left[\frac{\hat{h}\cdot{\vec{k}}_{T}}{M}\mathcal{H}^{(0)}\,g_{1T}D_{1}\right]\,,$
(S16)

where $\tilde{\mathcal{F}}$ is defined as

$\displaystyle\tilde{\mathcal{F}}[\omega\,\mathcal{H}\,g\,D]$
$\displaystyle\equiv
2z\sum_{f}\mathcal{H}_{f}(Q^{2})\int\\!\mathrm{d}^{2}{\vec{k}}_{T}\,\mathrm{d}^{2}{\vec{p}}_{T}\,\delta^{(2)}\bigl{(}{\vec{q}}_{T}+{\vec{k}}_{T}-{\vec{p}}_{T}\bigr{)}\times\omega(\vec{k}_{T},\vec{p}_{T})\,g_{f}(x,k_{T})\,D_{f}(z,p_{T})\,.$
(S17)

Here $\vec{q}_{T}=-\vec{P}_{hT}/z$, $\hat{h}=\vec{P}_{hT}/P_{hT}$, and
$\omega({\vec{k}}_{T},{\vec{p}}_{T})$ denotes the weight prefactors in Eq.
(IV.4) that depend on ${\vec{k}}_{T}$ and ${\vec{p}}_{T}$.

The leading-power SIDIS cross section differential in $q_{*}$ is

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*}}=\int_{0}^{\infty}\\!\mathrm{d}P_{hT}\,P_{hT}\int_{0}^{2\pi}\\!\mathrm{d}\phi_{h}\,\delta\bigl{(}q_{*}+\sin\phi_{h}P_{hT}/z\bigr{)}\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}^{2}\vec{P}_{hT}}\,.$
(S18)

Plugging Eq. (III) into Eq. (S18), since the structure functions themselves
have no dependence on $\phi_{h}$, terms which have angular prefactors that are
odd under $\phi_{h}\to\pi-\phi_{h}~{}$ vanish under the integral of
$\mathrm{d}\phi_{h}$. Thus we are left with

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*}}$
$\displaystyle=\sigma_{0}\int_{0}^{\infty}\\!\mathrm{d}P_{hT}\,P_{hT}\int_{0}^{2\pi}\\!\mathrm{d}\phi_{h}\,\delta\bigl{(}q_{*}+\sin\phi_{h}P_{hT}/z\bigr{)}\Bigl{\\{}W_{UU,T}+\lambda_{e}S_{L}\sqrt{1-\epsilon^{2}}\,W_{LL}$
$\displaystyle\>\;+\epsilon\cos(2\phi_{h})W_{UU}^{\cos(2\phi_{h})}+S_{T}\sin\phi_{h}\cos\phi_{S}\left(W_{UT,T}^{\sin(\phi_{h}\\!-\phi_{S})}+\epsilon
W_{UT}^{\sin(\phi_{h}\\!+\phi_{S})}\right)$
$\displaystyle\>\;+S_{T}\,\cos\phi_{S}\sin(3\phi_{h})W_{UT}^{\sin(3\phi_{h}\\!-\phi_{S})}+\lambda_{e}S_{T}\sqrt{1-\epsilon^{2}}\sin\phi_{h}\sin\phi_{S}W_{LT}^{\cos(\phi_{h}\\!-\phi_{S})}\Bigr{\\}}\,.$
(S19)

Especially notice that $W_{UL}^{\sin 2\phi_{h}}\sim
h_{1L}^{\perp}H_{1}^{\perp}$ does not appear in Eq. (IV.4) since its prefactor
$\sin(2\phi_{h})$ is odd under $\phi_{h}\to\pi-\phi_{h}$.

For each term in Eq. (S18), writing the $\phi_{h}$ dependent coefficient as
$\kappa(\phi_{h})$, we have

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*}}[\kappa(\phi_{h})\,W(\omega\mathcal{H}gD)]$
$\displaystyle=2z\sum_{f}\mathcal{H}_{f}(Q^{2})\int\\!\mathrm{d}^{2}\vec{P}_{hT}\,\mathrm{d}^{2}{\vec{k}}_{T}\,\mathrm{d}^{2}{\vec{p}}_{T}\,\delta^{(2)}\bigl{(}{\vec{q}}_{T}+{\vec{k}}_{T}-{\vec{p}}_{T}\bigr{)}\,\delta\bigl{(}q_{*}+\sin\phi_{h}P_{hT}/z\bigr{)}$
$\displaystyle\hskip
72.26999pt\times\kappa(\phi_{h})\,\omega(\vec{k}_{T},\vec{p}_{T})\,g_{f}(x,k_{T})\,D_{f}(z,p_{T})$
$\displaystyle=2z^{3}\sum_{f}\mathcal{H}_{f}(Q^{2})\int\\!\mathrm{d}^{2}{\vec{k}}_{T}\,\mathrm{d}^{2}{\vec{p}}_{T}\,\delta\bigl{(}q_{*}+\hat{y}\cdot({\vec{k}}_{T}-{\vec{p}}_{T})\bigr{)}$
$\displaystyle\hskip
72.26999pt\times\kappa({\vec{k}}_{T},{\vec{p}}_{T})\,\omega(\vec{k}_{T},\vec{p}_{T})\,g_{f}(x,k_{T})\,D_{f}(z,p_{T})\,,$
(S20)

where $\hat{y}$ is the unit vector in the Trento frame. In the second line,
all appearances of $\hat{h}$ in $\omega({\vec{k}}_{T},{\vec{p}}_{T})$ are
replaced by

$\displaystyle\hat{h}=\frac{{\vec{k}}_{T}-{\vec{p}}_{T}}{\lvert{\vec{k}}_{T}-{\vec{p}}_{T}\rvert}\,,$
(S21)

and $\kappa(\phi_{h})$ is replaced by a function of ${\vec{k}}_{T}$ and
${\vec{p}}_{T}$ according to the rule

$\displaystyle\cos\phi_{h}\to\hat{h}\cdot\hat{x}^{\prime}\,,\qquad\sin\phi_{h}\to\hat{h}\cdot\hat{y}\,,$
(S22)

where $\hat{x}^{\prime}$ and $\hat{y}$ are again Trento frame unit vectors.

Then we can write Eq. (IV.4) as

$\displaystyle\frac{1}{\sigma_{0}}\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z\,\mathrm{d}q_{*}}=\tilde{\mathcal{F}}_{*}[f_{1}D_{1}]+\lambda_{e}S_{L}\sqrt{1-\epsilon^{2}}\,\tilde{\mathcal{F}}_{*}[g_{1L}D_{1}]$
$\displaystyle\qquad+\epsilon\,\tilde{\mathcal{F}}_{*}\\!\left[\left[(\hat{h}\cdot\hat{x}^{\prime})^{2}-(\hat{h}\cdot\hat{y})^{2}\right]\frac{-2(\hat{h}\cdot{\vec{k}}_{T})(\hat{h}\cdot{\vec{p}}_{T})+{\vec{k}}_{T}\cdot{\vec{p}}_{T}}{MM_{h}}\,h_{1}^{\perp}H_{1}^{\perp}\right]$
$\displaystyle\qquad+S_{T}\cos\phi_{S}\Biggl{(}\tilde{\mathcal{F}}_{*}\left[(\hat{h}\cdot\hat{y})\frac{-\hat{h}\cdot{\vec{k}}_{T}}{M}f_{1T}^{\perp}D_{1}\right]+\epsilon\,\tilde{\mathcal{F}}_{*}\\!\left[(\hat{h}\cdot\hat{y})\frac{-\hat{h}\cdot{\vec{p}}_{T}}{M_{h}}h_{1}H_{1}^{\perp}\right]\Biggr{)}$
$\displaystyle\qquad+S_{T}\cos\phi_{S}\,\tilde{\mathcal{F}}_{*}\\!\left[\left[3(\hat{h}\cdot\hat{y})-4(\hat{h}\cdot\hat{y})^{3}\right]\frac{2(\hat{h}\cdot{\vec{k}}_{T})\,(\vec{k}_{T}\cdot\vec{p}_{T})+k_{T}^{2}\,(\hat{h}\cdot{\vec{p}}_{T})-4(\hat{h}\cdot{\vec{k}}_{T})^{2}\,(\hat{h}\cdot{\vec{p}}_{T})}{2M^{2}M_{h}}\,h_{1T}^{\perp}H_{1}^{\perp}\right]$
$\displaystyle\qquad+\lambda_{e}S_{T}\sqrt{1-\epsilon^{2}}\sin\phi_{S}\,\tilde{\mathcal{F}}_{*}\\!\left[(\hat{h}\cdot\hat{y})\frac{\hat{h}\cdot{\vec{k}}_{T}}{M}g_{1T}D_{1}\right]\,.$
(S23)

Note that we have used $\sin(3\phi_{h})=3\sin\phi_{h}-4\sin^{3}\phi_{h}$. Here
we define $\tilde{\mathcal{F}}_{*}$ as

$\displaystyle\tilde{\mathcal{F}}_{*}[\omega_{*}\,\mathcal{H}\,g\,D]\equiv
2z^{3}\sum_{f}\mathcal{H}_{f}(Q^{2})\int\\!\mathrm{d}^{2}{\vec{k}}_{T}\,\mathrm{d}^{2}{\vec{p}}_{T}\,\delta\bigl{(}q_{*}+\hat{y}\cdot{\vec{k}}_{T}-\hat{y}\cdot{\vec{p}}_{T}\bigr{)}\,\omega_{*}(\vec{k}_{T},\vec{p}_{T})\,g_{f}(x,k_{T})\,D_{f}(z,p_{T})\,,$
(S24)

where the weight function $\omega_{*}({\vec{k}}_{T},{\vec{p}}_{T})$ includes
the full prefactor structures, and can be understood as product of
$\kappa({\vec{k}}_{T},{\vec{p}}_{T})$ and
$\omega({\vec{k}}_{T},{\vec{p}}_{T})$ in Eq. (IV.4).

###  IV.5 Resolution curves for all detector regions

Fig. 2 in the main manuscript was restricted to the case of a backward
electron and central pion, which features the largest share of the total pion
sample for our selection cuts. In Fig. S1, we provide additional results for
the expected detector resolution of different SIDIS TMD observables in all
other relevant detector regions. Note the change in vertical scale compared to
Fig. 2. We find that a clear improvement in resolution from using $q_{*}$
persists across all detector regions. The case of electrons in the forward
detector region (i.e., from backscattered electrons at very large $Q^{2}\to
s$) has a negligible contribution to the total rate. We note that the
improvement of $q_{*}$ in resolution compared to $P_{hT}$ deteriorates
slightly in the cases where the hadron is backward for $P_{hT}<5\,{\rm GeV}$,
and actually features worse resolution for $P_{hT}>5\,{\rm GeV}$. This is
expected as the fixed angular resolution
$\sigma_{\theta_{h}^{\mathrm{EIC}}}=0.001$ translates to a wider range in
pseudorapidity as $\theta_{h}^{\mathrm{EIC}}\to\pi$.

Figure S1: Expected event-level detector resolution $\sigma_{\mathcal{O}}$ for
different SIDIS TMD observables $\mathcal{O}$. We show relative resolutions as
a function of the magnitude of $\mathcal{O}=q_{*}$ (solid red), $P_{hT}/z$
(dashed blue), and $P_{hT}/z\sin\phi$ (dotted green), for backward electrons
(left column) and central electrons (right column) as well as forward pions in
the top row, central pions in the center row, and backward pions in the bottom
row. The center left panel corresponds to Fig. 2 in the main manuscript.

###  IV.6 Details on Bayesian reweighting

####  IV.6.1 Theory templates

In this section we describe in detail how the pseudodata $d_{i}$ and theory
replicas $t_{i}(\omega_{i})$ in the main text are constructed. They are
defined as the normalized bin-integrated spectrum for the observable
$\mathcal{O}=P_{hT}/z$ or $\mathcal{O}=q_{*}$ being tested,

$\displaystyle\\{d_{n},t_{n}(w_{i})\\}=\frac{1}{\mathrm{d}\sigma/(\mathrm{d}x\,\mathrm{d}z\,\mathrm{d}Q^{2})}\int_{\mathcal{O}_{a,n}}^{\mathcal{O}_{b,n}}\\!\mathrm{d}\mathcal{O}\,\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}z\,\mathrm{d}Q^{2}\,\mathrm{d}\mathcal{O}}\,,$
(S25)

where the right-hand side is evaluated with the appropriate values of
$\omega_{i}$ (either central or chosen according to the Monte-Carlo replica)
inserted into Eq. (IV). Restricting Eqs. (III) and (III) to the unpolarized
contribution from $W_{UU,T}$, the differential leading-power spectrum
explicitly reads

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}z\,\mathrm{d}Q^{2}\,\mathrm{d}\mathcal{O}}=\frac{2yz^{3}}{Q^{2}}\,\sigma_{0}\int_{0}^{\infty}\\!\mathrm{d}b_{T}\,K_{\mathcal{O}}(\mathcal{O}b_{T})\sum_{f}\mathcal{H}_{f}(Q^{2},\mu)\,\tilde{f}_{1f}(x,b_{T},\mu,\zeta)\,\tilde{D}_{1f}(z,b_{T},\mu,\zeta)\,,$
(S26)

where $\mu$ is the $\overline{\text{MS}}$ scale that the TMD PDF and FF are
being evolved to, $\zeta=Q^{2}$ is the Collins-Soper scale, and

$\displaystyle K_{P_{hT}/z}(x)=xJ_{0}(x)\,,\qquad
K_{q_{*}}=\frac{2\cos(x)}{\pi}\,,$ (S27)

is the integral kernel that depends on the respective observable. At tree-
level, the hard function simply reduces to the electric charge,
$\mathcal{H}_{f}=e_{f}^{2}$. Using tree-level matching onto collinear PDFs and
FFs, the evolved TMDs are given by

$\displaystyle\tilde{f}_{1f}(x,b_{T},\mu,\zeta)$
$\displaystyle=f_{1f}(x,\mu)\,U(b_{T},\mu_{0},\zeta_{0},\mu,\zeta)\,\tilde{f}_{1}^{\,\mathrm{NP}}(x,b_{T})\,,$
$\displaystyle\tilde{D}_{1f}(z,b_{T},\mu,\zeta)$
$\displaystyle=D_{1f}(z,\mu)\,U(b_{T},\mu_{0},\zeta_{0},\mu,\zeta)\,\tilde{\tilde{D}}_{1}^{\mathrm{NP}}(z,b_{T})\,,$
(S28)

where $\tilde{f}_{1}^{\,\mathrm{NP}}$ and
$\tilde{\tilde{D}}_{1}^{\mathrm{NP}}$ are the nonperturbative model function
of interest and satisfy
$\tilde{f}_{1}^{\,\mathrm{NP}}(x,b_{T}),\tilde{D}_{1}^{\mathrm{NP}}(z,b_{T})=1+\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})$.
Here we have used that at our leading-logarithmic working order, we are free
to evaluate the collinear PDF $f_{1f}(x,\mu)$ and collinear FF $D_{1f}(z,\mu)$
at the high scale $\mu\sim Q$. The evolution factor $U$ accounts for the
virtuality and Collins-Soper evolution of the TMD PDF and at leading-
logarithmic order is given by

$\displaystyle
U(b_{T},\mu_{0},\zeta_{0},\mu,\zeta)=\exp\Bigl{[}-\int_{\mu_{0}}^{\mu}\\!\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\frac{\alpha_{s}(\mu^{\prime})}{2\pi}\,\Gamma_{0}^{q}\ln\frac{\sqrt{\zeta}}{\mu^{\prime}}\Bigr{]}\,\exp\Bigl{[}\frac{1}{2}\tilde{\gamma}_{\zeta}^{\mathrm{NP}}(b_{T})\ln\frac{\zeta}{\zeta_{0}}\Bigr{]}\,,$
(S29)

where $\Gamma_{0}^{q}=4C_{F}$ is the one-loop quark cusp anomalous dimension
and
$\tilde{\gamma}_{\zeta}^{\mathrm{NP}}(b_{T})=\mathcal{O}(\Lambda_{\mathrm{QCD}}^{2}b_{T}^{2})$
is the nonperturbative contribution to the Collins-Soper kernel. The initial
scales $\mu_{0}$ and $\zeta_{0}$ define the boundary condition at which the
nonperturbative TMDs are defined. We use the prescription Bacchetta  _et al._
(2020, 2022)

$\displaystyle\mu_{0}(b_{T})=\frac{b_{0}}{b_{\mathrm{max}}}\Bigl{(}\frac{1-e^{-b_{T}^{4}/b_{\mathrm{max}}^{4}}}{1-e^{-b_{T}^{4}/b_{\mathrm{min}}^{4}}}\Bigr{)}^{-1/4}\,,\quad
b_{0}=2e^{\gamma_{E}}\approx 1.123\,,\quad
b_{\mathrm{max}}=b_{0}\,\mathrm{GeV}^{-1}\,,\quad
b_{\mathrm{min}}=\frac{b_{0}}{\mu}\,,\quad\zeta_{0}=1\,\mathrm{GeV}^{2}\,.$
(S30)

This ensures that the coupling $\alpha_{s}(\mu^{\prime})$ under the integral
in Eq. (S29) remains perturbative due to $\mu^{\prime}\geq
b_{0}/b_{\mathrm{max}}=1\,\mathrm{GeV}$. In addition, $\mu_{0}(b_{T}\to 0)\to
b_{0}/b_{\mathrm{min}}=\mu$ together with $\gamma_{\zeta}^{\mathrm{NP}}\to 0$
ensures that $U\to 1$ for $b_{T}\to 0$. We use N3LL $\alpha_{s}$ evolution
with $n_{f}=5$ active flavors and evolve to $\mu=\sqrt{\zeta}=Q$ throughout
our predictions.

Several further simplifications apply to Eq. (S25). We may analytically
perform the bin integral over
$\mathcal{O}_{a}\leq\mathcal{O}\leq\mathcal{O}_{b}$,

$\displaystyle\tilde{K}_{\mathcal{O}}(\mathcal{O}_{b}b_{T})-\tilde{K}_{\mathcal{O}}(\mathcal{O}_{a}b_{T})\equiv
b_{T}\int_{\mathcal{O}_{a}}^{\mathcal{O}_{b}}\\!\mathrm{d}\mathcal{O}\,K_{\mathcal{O}}(\mathcal{O}b_{T})\,,\qquad\tilde{K}_{P_{hT}/z}(x)=xJ_{1}(x)\,,\qquad\tilde{K}_{q_{*}}(x)=\frac{2\sin(x)}{\pi}\,.$
(S31)

An important limit is the total integral over $\mathcal{O}$, which involves
the integral kernel

$\displaystyle\lim_{\mathcal{O}_{b}\to\infty}\frac{1}{b_{T}}\tilde{K}_{P_{hT}/z}(\mathcal{O}_{b}b_{T})=\lim_{\mathcal{O}_{b}\to\infty}\frac{1}{b_{T}}\tilde{K}_{q_{*}}(\mathcal{O}_{b}b_{T})=\delta(b_{T})\,.$
(S32)

This allows us to analytically evaluate the total cross section that appears
in the denominator in Eq. (S25),

$\displaystyle\frac{\mathrm{d}\sigma}{\mathrm{d}x\,\mathrm{d}z\,\mathrm{d}Q^{2}}=\frac{2yz^{3}}{Q^{2}}\,\sigma_{0}\sum_{f}e_{f}^{2}\,f_{1f}(x,\mu)\,D_{1f}(z,\mu)\,,$
(S33)

which recovers the tree-level SIDIS total cross section. Taking the ratio in
Eq. (S25), many factors drop out, including in particular the sum over
collinear PDFs and FFs, leaving behind only the flavor-independent model
functions and TMD evolution,

$\displaystyle\\{d_{n},t_{n}(\omega_{i})\\}=\int_{0}^{\infty}\\!\frac{\mathrm{d}b_{T}}{b_{T}}\,\Bigl{[}\tilde{K}_{\mathcal{O}}(b_{T}\mathcal{O}_{b,n})-\tilde{K}_{\mathcal{O}}(b_{T}\mathcal{O}_{a,n})\Bigr{]}\,\bigl{[}U(b_{T},\mu_{0},\zeta_{0},\mu,\zeta)\bigr{]}^{2}\,\tilde{f}_{1}^{\,\mathrm{NP}}(x,b_{T})\,\tilde{D}_{1}^{\mathrm{NP}}(z,b_{T})\,,$
(S34)

for $\mathcal{O}=P_{hT}/z$ and $\mathcal{O}=q_{*}$. It is easy to see from Eq.
(S34) and Eq. (S32) that the pseudodata and theory templates indeed satisfy
$\sum_{n}d_{n}=\sum_{n}t_{n}(\omega_{i})=1$, where we note that we include an
overflow bin $\mathcal{O}\geq 4\,\mathrm{GeV}$ in all sums over $n$, including
in particular the likelihood function used in the main text.

The simple form of the normalized $P_{hT}/z$ and $q_{*}$ spectrum in Eq. (S34)
in terms of the underlying nonperturbative TMD distributions ensures that our
conclusions about the TMD sensitivity are fully general, even though we
consider fixed $x$ and $z$. This is because $x$ and $z$ only enter through
$\tilde{f}_{1}^{\,\mathrm{NP}}(x,b_{T})$ and
$\tilde{D}_{1}^{\mathrm{NP}}(z,b_{T})$, which for fixed $x$ and $z$ map onto
specific values of $\omega_{\\{1,2,3\\}}$ and $\alpha$ in Eq. (IV) as
discussed below, and hence only affect the initial values of our study. We
stress again that the analytic relation we exploited between the cross
sections in TMD and collinear factorization is due to our leading-logarithmic
working order and the specific perturbative scheme choices we made, but the
simplifications above are fully justified for an exploratory study like this
one. In particular, Eq. (S34) still retains all the key characteristic
features of a realistic TMD study such as the correct amount of Sudakov
suppression at long distances.

####  IV.6.2 Nonperturbative model parameters and Bayesian priors

We work with the nonperturbative model used in the MAPTMD22 global fit of TMD
PDFs and FFs from Drell-Yan and SIDIS data Bacchetta  _et al._ (2022). The
model for the CS kernel reads

$\displaystyle\tilde{\gamma}_{\zeta}^{\mathrm{NP}}(b_{T})=g_{K}(b_{T}^{2})=-g_{2}^{2}\frac{b_{T}^{2}}{2}$
(S35)

Since our reweighting study is performed at fixed $Q$, we do not expect to be
sensitive to CS evolution between different values of $Q$ and hold the
parameter $g_{2}=(0.248\pm 0.008)\,\mathrm{GeV}$ fixed at its central value.
(A full four-dimensional measurement differential in $x$, $z$, $Q$, and
$q_{*}$ would of course exhibit the usual sensitivity to the CS kernel.)

For fixed $x$ and $z$, the nonperturbative models of Bacchetta  _et al._
(2022) for the TMD PDF and FF relate to Eq. (IV) as

$\displaystyle\omega_{1}=\frac{g_{1}(x)}{4}\,,\qquad\omega_{2}=\frac{g_{3}(z)}{4z^{2}}\,,\qquad\omega_{3}=\frac{g_{3B}(z)}{4z^{2}}\,,\qquad\alpha=\frac{g_{3}(z)}{g_{3}(z)+\lambda_{F}\,g_{3B}^{2}(z)/z^{2}}\,,$
(S36)

where $g_{\\{1,3,3B\\}}$ are functions of $x$ and $z$ defined in terms of the
underlying model parameters in Bacchetta  _et al._ (2022). In the comparison,
we have used that the fit result of Bacchetta  _et al._ (2022) down to $x\geq
0.1$ is compatible with a single Gaussian for the TMD PDF to a good
approximation, and have set $\lambda=\lambda_{2}=0$. For $x=0.1$ and $z=0.15$
as considered in the main text, the first three parameters evaluate to

$\displaystyle\omega_{1}=(0.0791\pm
0.0063)\,\mathrm{GeV}^{2}\,,\qquad\omega_{2}=(0.0167\pm
0.0059)\,\mathrm{GeV}^{2}\,,\qquad\omega_{3}=(0.8153\pm
0.0637)\,\mathrm{GeV}^{2}\,,$ (S37)

In our reweighting study, we take the variance of the Gaussian prior
probability distributions for $\\{\omega_{1},\omega_{2},\omega_{3}\\}$ to be
numerically equal to the $1\sigma$ confidence intervals above. (We ignore
nondiagonal entries in the covariance matrix of Bacchetta  _et al._ (2022),
which were found to be small there.) We hold $\alpha=0.079\pm 0.029$ fixed at
its central value for simplicity. We note that we have also performed the
reweighting study using a simplified version of the nonperturbative model used
in Scimemi and Vladimirov (2020), which features an approximately exponential
term at large distances, and have arrived at similar conclusions.

####  IV.6.3 Parametrizing systematic bias

$X$ | $p_{e}/\\!\,\mathrm{GeV}$ | $\eta_{e}$ | $p_{h}/\\!\,\mathrm{GeV}$ | $\eta_{h}$
---|---|---|---|---
$\langle X\rangle$ | 19.86 | -0.472 | 4.23 | 0.892
$\Delta X$ | 0.32 | 0.011 | 2.18 | 0.275
Table S1: Mean and variance of kinematic distributions in
$X=\\{p_{e},p_{h},\eta_{e},\eta_{h}\\}$ for a bin with $0.085\leq x\leq
0.115,0.13\leq z\leq 0.17,400\,\mathrm{GeV}^{2}\leq Q^{2}\leq
401\,\mathrm{GeV}^{2}$ as relevant for the biased event weights in Eq. (S39).

Here we describe how the bias strengths shown in Fig. 4 in the main text are
evaluated. To model the effect of a momentum miscalibration $\delta_{p}$, we
use the Pythia sample of pions described in the main text, restrict to a bin
$0.085\leq x\leq 0.115,0.13\leq z\leq 0.17,400\,\mathrm{GeV}^{2}\leq Q^{2}\leq
401\,\mathrm{GeV}^{2}$ centered on our choices for $x,z,Q^{2}$ in the main
text, and evaluate the $\mathcal{O}=P_{hT}/z$ or $\mathcal{O}=q_{*}$ spectrum
after replacing

$\displaystyle p_{e}\to(1+\delta_{p_{e}})\,p_{e}$ (S38)

in the event record. Note that we only use the biased event record to
calculate $\mathcal{O}$, but continue to cut on the true values of $x$, $z$,
$Q^{2}$ since we expect the impact of $\delta_{p_{e}}$ on the reference Born
kinematics to be subleading compared to the direct effect on the reconstructed
$\mathcal{O}$. Similarly, to model the effect of a (generic) non-uniform
detector response that changes as a function of
$X=\\{p_{e},p_{h},\eta_{e},\eta_{h}\\}$, we apply an additional weight
$\epsilon(X)$ to each event calculated as

$\displaystyle\epsilon(X)=1+\Delta\epsilon_{X}\frac{X-\langle X\rangle}{\Delta
X}\,,\qquad\Delta X=\sqrt{\langle X^{2}\rangle-\langle X\rangle^{2}}\,,$ (S39)

where the $\langle\dots\rangle$ refers to an average over all events in the
current $x,z,Q^{2}$ bin. Normalizing the effect of $\Delta\epsilon_{X}$ to the
variance of the distribution ensures that the impacts of $\Delta\epsilon_{X}$
for different $X$ are comparable even when individual $X$ have more or less
narrow underlying distributions or different mass dimension. The explicit
values of $\langle X\rangle$, $\Delta X$ for the bin we consider are collected
in Table S1 for reference.

To obtain the final biased theory templates, we take one of the bias
parameters
$B=\\{\Delta\epsilon_{p_{e}},\Delta\epsilon_{\eta_{e}},\Delta\epsilon_{p_{h}},\Delta\epsilon_{\eta_{h}},\delta_{p_{e}}\\}$
to be nonzero at a time and evaluate the biased normalized pseudodata as

$\displaystyle
d_{n}^{\mathrm{bias}}=d_{n}\frac{h_{n}^{\mathrm{bias}}}{h_{n}}\biggl{(}\sum_{m}d_{m}\frac{h_{m}^{\mathrm{bias}}}{h_{m}}\biggr{)}^{-1}\,,$
(S40)

where $h_{n}$ ($h_{n}^{\mathrm{bias}}$) is the Pythia result for the
normalized $\mathcal{O}$ spectrum in the same bins with vanishing (nonzero)
bias. We then insert the biased pseudodata from Eq. (S40) into the reweighting
analysis, using unbiased theory templates, and evaluate the biased mean values
$\omega_{i}^{\mathrm{bias}}$ of the nonperturbative physics parameters of
interest. Repeating this for several values of the bias parameters, we can
evaluate the required partial derivatives (“bias strengths”)
$\partial\omega_{i}^{\mathrm{bias}}/\partial B$ with respect to each of the
bias parameters $B$ by finite differences, leading to the results in Fig. 4.