aainstitutetext: Departamento de Física Teórica and IPARCOS, Universidad
Complutense de Madrid (UCM), E-28040 Madrid, Spainbbinstitutetext: Theoretical
Division, MS B283, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
ccinstitutetext: Institute for Theoretical Physics Amsterdam and Delta
Institute for Theoretical Physics University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlandsddinstitutetext: Nikhef, Theory Group,
Science Park 105, 1098 XG, Amsterdam, The Netherlands

# Probing Transverse-Momentum Distributions With Groomed Jets

Daniel Gutierrez-Reyes b Yiannis Makris b Varun Vaidya a Ignazio Scimemi c,d
and Lorenzo Zoppi<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present the transverse momentum spectrum of groomed jets in di-jet events
for $e^{+}e^{-}$ collisions and semi-inclusive deep inelastic scattering
(SIDIS). The jets are groomed using a soft-drop grooming algorithm which helps
in mitigating effects of non-global logarithms and underlying event. At the
same time, by reducing the final state hadronization effects, it provides a
clean access to the non-perturbative part of the evolution of transverse
momentum dependent (TMD) distributions. In SIDIS experiments we look at the
transverse momentum of the groomed jet measured w.r.t. the incoming hadron in
the Breit frame. Because the final state hadronization effects are
significantly reduced, the SIDIS case allows to probe the TMD parton
distribution functions. We discuss the sources of non-perturbative effects in
the low transverse momentum region including novel (but small) effects that
arise due to grooming. We derive a factorization theorem within SCET and resum
any large logarithm in the measured transverse momentum up to NNLL accuracy
using the $\zeta$-prescription as implemented in the artemide package and
provide a comparison with simulations.

††preprint: Nikhef: 2019-034, LA-UR-19-31102

##  1 Introduction

The understanding of the structure of nucleons is one of the most important
and interesting research subject in modern nuclear physics. The ultimate goal
would be to have a complete description of quarks/gluons position and momenta
inside a hadron, which is not easy because of the entanglement of
initial/final states in all hadronic processes. In order to properly define
all hadron constituent contributions, the cross sections should be factorized
in some region of the phase space into properly defined hadronic matrix
elements. Here we will consider the transverse momentum dependent
distributions (TMD), which appear in the factorization of several processes
like Drell-Yan, semi-inclusive deep-inelastic-scattering (SIDIS) and
$e^{+}e^{-}$ hadron production Collins:2011zzd ; GarciaEchevarria:2011rb ;
Echevarria:2014rua ; Gaunt:2014ska ; Vladimirov:2017ksc . Drell-Yan processes
directly test the TMD parton distribution functions (TMDPDF), while in SIDIS
cross sections the TMDPDFs are coupled to a TMD fragmentation function in the
final state. Finally in $e^{+}e^{-}$ hadron production only the TMD
fragmentation is present. Because of the factorization theorem, the TMDs have
several universal features like rapidity and renormalization scale evolution,
which should be also tested including their (universal) non-perturbative part.
Recently some of us have considered the possibility to define a jet-TMD,
replacing a final state hadron with a jet Neill:2016vbi ; Gutierrez-
Reyes:2018qez ; Gutierrez-Reyes:2019vbx  in SIDIS and $e^{+}e^{-}$ processes.
The check of this possibility has revealed that standard jet definitions are
compatible with a factorization theorem only in the case of small enough
radii, which is a not obvious experimental condition in the planned electron-
hadron collider like EIC or LHeC. Instead large jet-radii need a specific
definition of jet, which allows soft radiation to be independent of radius. In
Gutierrez-Reyes:2018qez ; Gutierrez-Reyes:2019vbx  this was achieved using the
winner-take-all (WTA) axis Bertolini:2013iqa , and the perturbative
calculations were done with a precision similar to the case of fragmenting
hadrons.

In this work we consider the possibility of groomed jets in SIDIS or
$e^{+}e^{-}\to\text{2\;jets}$. Developments in jet substructure have shown
that applying a grooming algorithm to a jet, using for example the so called
“soft-drop” procedure, robustly removes the contamination from both underlying
event and non-global correlations. Since this process essentially removes wide
angle soft radiation, retaining only a collinear core, it also dramatically
reduces hadronization effects (see fig. 1), thus allowing an easier access to
the TMD non-perturbative physics which we want to probe. Groomed jets with an
identified light/heavy hadron in the jet were also proposed as probes of TMD
evolution and distribution in Makris:2017arq ; Makris:2018npl . The residual
non-perturbative effects contain pieces that depend on the soft-drop grooming
procedure and require careful analysis as was pointed out in Hoang:2019ceu .
In addition, with the use of soft-drop we can derive factorization theorems
for large jet radius ($R\sim 1$), which we consider to be the relevant case
for low energy experiments, such as EIC. In order to focus on collimated jet
configurations, we also impose an upper cutoff in the groomed jet invariant
mass.111Note that the small transverse momentum constraint does not
necessarily ensure collimated configurations since topologies with two or more
widely separated sub-jets are also permitted. This constraint allows us to
derive a factorization theorem involving the same universal soft function that
appears in traditional hadronic TMD, and which is independent of the jet
radius for $R\sim 1$. This is a key feature for groomed jets and it is
necessary for the universality of TMDs and for this reason, in this paper, we
only consider $R\sim 1$. The cutoff is imposed using groomed jet-thrust,
$e\equiv(m/Q)^{2}$, where $m$ is the groomed invariant mass and $Q$ is the
center of mass energy. This allows us to introduce a single cutoff parameter,
$e_{\text{cut}}$, independent of the jet energy or transverse momentum.

Figure 1: Hadronization effects in a typical $e^{+}e^{-}\to\text{2\;jets}$
from Pythia 8 Sjostrand:2006za ; Sjostrand:2007gs . at center of mass values,
Left: $Q=100$ GeV, Right: $Q=50$ GeV.

The paper is organized as follows. In sec. 2, we give a review of soft-drop
and discuss the factorization of the $e^{+}e^{-}\to\text{2\;jets}$, transverse
momentum decorrelation within SCET and give detailed comparisons of our NNLL
accurate prediction with simulations for this observable. In sec. 3, we
consider the factorization for the corresponding observable in DIS. We
carefully enumerate all the non-perturbative corrections and discuss their
universality in sec. 4. We conclude in sec. 5. The details of the one loop
calculations, resummation, and evolution are provided in the appendix.

##  2 Di-jet events in electron-positron colliders

In this section we discuss the measurement of momentum de-correlation in
electron-positron colliders. We identify events with two final state jets and
we consider the transverse momentum of one jet w.r.t. the other. The
measurement that we are considering in this work is a generalization of the
di-hadron momentum de-correlation,

${\boldsymbol{q}}_{T}=\frac{{\boldsymbol{p}}_{Th_{1}}}{z_{1}}+\frac{{\boldsymbol{p}}_{Th_{2}}}{z_{2}}$
(1)

where one or both of the identified hadrons is replaced by a jet, defined
through an infrared-safe jet algorithm. Here ${\boldsymbol{p}}_{Th_{i}}$ and
$z_{i}$ are the transverse momentum and energy fraction of the hadron $i$
respectively. The factorization theorem is usually written for this normalized
vector sum of the transverse momenta rather than just the sum of the
transverse momenta. It can be verified by momentum conservation and simple
geometry that the quantity ${\boldsymbol{p}}_{Th_{1}}/z_{1}$ represents the
transverse momentum of the radiation recoiling against the hadron w.r.t the
axis defined by the hadron itself. This makes it convenient to write a
factorization theorem which matches onto the standard hadron fragmentation
function as explained in Makris:2017arq .

We consider three possible scenarios as illustrated in fig. 2 and we refer to
them as di-hadron, hadron-jet, and di-jet momentum de-correlation. To simplify
the discussion we focus on the case of di-jets (fig. 2c) and we briefly
comment how our results are generalized for the case of hadron-jet de-
correlation. For the case of groomed jets the observable
${\boldsymbol{q}}_{T}$ is defined with the groomed quantities, i.e.,
$p_{\text{jet}}^{\mu}$ is the _groomed_ jet four-momentum and
$z_{\text{jet}}=2p^{0}_{\text{jet}}/Q$. The transverse component
${\boldsymbol{p}}_{T\;\text{jet}}$ is measured with respect to an axis close
to the full or groomed jet axes. The exact choice of the axis only differs by
power corrections. For concreteness in the results that follow we make the
choice of the axis to lie along one of the groomed jets.

Figure 2: Three possible transverse momentum de-correlation measurements in
$e^{+}e^{-}$ annihilation: (a) Identify two hadrons $h_{1}$ and $h_{2}$ with
momenta $p_{h_{1}},p_{h_{2}}$ and energy fractions $z_{h_{1}},z_{h_{2}}$
respectively, (b) Identify a jet and a hadron with momenta $p_{\text{jet}}$,
$p_{h}$ with energy fractions $z_{\text{jet}},z_{h}$, (c) Identify two jets
with momenta $p_{\text{jet}_{1}},p_{\text{jet}_{2}}$ and energy fractions
$z_{\text{jet}_{1}},z_{\text{jet}_{2}}$.

Since we want to probe the non-perturbative physics, we wish to work in the
small transverse momentum regime ($q_{T}\ll\sqrt{s}$ where
$q_{T}\equiv|\boldsymbol{q}_{T}|$). There are various ways one can define the
jet axis and the choice of definition will impact the form of factorization.
It was discussed in ref. Gutierrez-Reyes:2018qez  that the standard jet axis
choice suffers from factorization breakdown for large jet radius. This
breakdown is due to energetic emissions at relatively wide angles. Such
configurations will contribute to the small transverse momentum region when
the energetic subjets are clustered in a single large radius jet. To avoid
this problem in refs. Gutierrez-Reyes:2018qez ; Gutierrez-Reyes:2019vbx  the
winner-take-all (WTA) axis was used instead. This way ensures that wide angle
energetic emissions induce large transverse momentum ($q_{T}\sim\sqrt{s}$)
pushing the $q_{T}$ measurement away from the observable region.

In this paper we propose, alternatively, the use of groomed jet-substructure
to isolate the collimated configurations and choose the jet axis to be the
groomed jet axis which is insensitive to jet boundary effects. Particularly we
consider the normalized jet mass as the relevant jet-substructure observable,

$e\equiv\Big{(}\frac{m_{J}}{Q}\Big{)}^{2}\ .$ (2)

We shall see that imposing this constraint still allows us to capture a
majority of events and hence does not significantly impact the cross-section.

###  2.1 Soft-drop: a brief review

The grooming procedure that we use is the soft-drop algorithm. We give here a
brief review of the soft-drop groomer and eventually discuss the various
hierarchies, the relevant modes and the factorization of the cross section in
the next sections.

Soft-drop grooming Larkoski:2014wba  removes contaminating soft radiation from
the jet by constructing an angular ordered tree of the jet, and removing the
branches at the widest angles which fail an energy requirement. The angular
ordering of the jet is constructed through the Cambridge/Aachen (C/A)
clustering algorithm Ellis:1993tq ; Catani:1993hr ; Dokshitzer:1997in ;
Wobisch:1998wt ; Wobisch:2000dk . As soon as a branch is found that passes the
test, it is declared the groomed jet, and all the constituents of the branch
are the groomed constituents. At the end of the grooming procedure only the
narrow energetic core remains from the original jet. Since at large angles all
collinear energetic radiation is to be found at the center of the jet, no cone
is actually imposed to enclose this core. One simply finds the branch whose
daughters are sufficiently energetic. Formally the daughters could have any
opening angle, though their most likely configuration is collinear.

The strict definition of the algorithm is as follows. Given an ungroomed jet
(which itself is identified first using a suitable algorithm such as the
anti-$k_{T}$, Cacciari:2008gp ), first we build the clustering history by
starting with a list of particles in the jet. At each stage we merge the two
particles within the list that are closest in angle222This merging is usually
taken to be summing the momenta of the particles, though one could use winner-
take-all schemes Salam:WTAUnpublished ; Bertolini:2013iqa ; Larkoski:2014uqa
.. This gives a pseudo-particle, and we remove the two daughters from the
current list of particles, replacing them with the merged pseudo-particle.
This is repeated until all particles are merged into a single parent. Then we
open the tree back up working backwards so that at each stage of the
declustering, we have two branches available, label them $i$ and $j$. We
require:

$\displaystyle\frac{\text{min}\\{E_{i},E_{j}\\}}{E_{i}+E_{j}}>z_{\text{cut}}\Big{(}\frac{\theta_{ij}}{R}\Big{)}^{\beta},$
(3)

where $z_{\text{cut}}$ is the modified mass drop parameter, $\beta$ is the
parameter which controls the angularities, $\theta_{ij}$ is the angle between
$i^{th}$ and $j^{th}$ particle, $R$ is the jet radius and $E_{i}$ is the
energy of the branch $i$. If the two branches fail this requirement, the
softer branch is removed from the jet, and we decluster the harder branch,
once again testing eq. (3) within the hard branch. The pruning continues until
we have a branch that when declustered passes the condition eq. (3). All
particles contained within this branch whose daughters are sufficiently
energetic constitute the groomed jet. Intuitively we have identified the first
genuine collinear splitting.

For a hadron-hadron collision, one uses the transverse momentum $(p_{T})$ with
respect to the beam for the condition of eq. (3),

$\displaystyle\frac{\text{min}\\{p_{Ti},p_{Tj}\\}}{p_{Ti}+p_{Tj}}>z_{\text{cut}}\Big{(}\frac{\theta_{ij}}{R}\Big{)}^{\beta}.$
(4)

We formally adopt the power counting $z_{\text{cut}}\ll 1$, though typically
one chooses $z_{\text{cut}}\sim 0.1$. See Marzani:2017mva  for a study on the
magnitude of the power corrections with respect to $z_{\text{cut}}$ for jet
mass distributions. To be specific, in this paper we consider only the case
$\beta=0$. Note that for $\beta=0$ the energy of the groomed jet constituents
is a collinear unsafe observable Larkoski:2014wba ; Larkoski:2015lea ,
however, the additional constraint of the measured transverse momentum
${\boldsymbol{q}}_{T}$ provides a physical collinear cutoff in a similar way a
jet shape measurement does. For detailed discussion on this we refer to
appendix D.

###  2.2 Hierarchies, modes, and factorization

In order to compute the transverse momentum de-correlation
${\boldsymbol{q}}_{T}$, defined in eq. (1), for two groomed jets in di-jet
events in $e^{+}e^{-}$ annihilation (fig. 2 (c)) we are going to impose a
normalized jet mass measurement as defined in eq. (2) on both jets. The other
parameters that enter our cross section are the soft-drop parameters
$z_{\text{cut}}\sim 0.1$, $\beta=0$. Ultimately we are going to integrate over
the jet mass measurement up to an appropriate (but still small) cut-off value
$e_{\text{cut}}$.

We have a rich spectrum of possible hierarchies of momenta, which are all
consistent with maintaining $q_{T}/Q,e_{\text{cut}}$, $z_{\text{cut}}\ll 1$.
We have that $q_{T}/Q,e_{\text{cut}}$, $z_{\text{cut}}$ are now expansion
parameters in the effective field theory (EFT), and they should be taken into
account in the factorization of the process. We first list and briefly discuss
these hierarchies and the corresponding factorization theorems within an EFT.
The general modes that we will consider will fall into three classes. Modes
that explicitly pass soft drop (usually the highly energetic collinear modes),
modes that explicitly fail soft-drop (the global soft function modes) and
finally those which can live on the border and need to be tested, as to
whether they pass or fail. Only the modes that pass soft-drop will contribute
to $e$, while ${\boldsymbol{q}}_{T}$ receives contributions from all radiation
that fails soft-drop.

The first regime in which we are interested is $Q\gg Qz_{\text{cut}}\gg
q_{T}\gtrsim Q\sqrt{e}\gg Q\sqrt{ez_{\text{cut}}}$. Here we have low values of
$q_{T}$ which are of the order of $Q\sqrt{e}$. We identify the following modes
to be relevant to the cross section:

soft: $\displaystyle\;\;\;p_{s}^{\mu}\sim q_{T}(1,1,1);$ collinear:
$\displaystyle\;\;\;p_{c}^{\mu}\sim
Q(\lambda_{c}^{2},1,\lambda_{c}),\;\lambda_{c}=\sqrt{e},$ (5)

and the factorization of the cross section in this region is schematically

$\frac{d\sigma}{de_{1}de_{2}d\boldsymbol{q}_{T}}=H^{ij}_{2}(Q;\mu)\times
S(\boldsymbol{q}_{T})\otimes\mathcal{J}^{\perp}_{i}(e_{1},Q,z_{\text{cut}},\boldsymbol{q}_{T})\otimes\mathcal{J}^{\perp}_{j}(e_{2},Q,z_{\text{cut}},\boldsymbol{q}_{T}).$
(6)

Apart from the hard factor $H$ all the other terms in this equation are
affected by rapidity divergences. The global soft function $S$ that appears in
the factorization theorem in eq. (6) (and later in the SIDIS case eq. (32)) is
the universal function that is also present in the factorization theorem of
Drell-Yan, di-hadron production in electron-positron annihilation, and semi-
inclusive DIS with TMDs. The operator definition of the soft function (see
refs. Collins:2011zzd ; GarciaEchevarria:2011rb ; Chiu:2012ir ) is given by

$S(\boldsymbol{q}_{T})=\frac{1}{N_{R}}\text{tr}\;{\langle}[S_{n}^{{\dagger}}S_{\bar{n}}](0)\delta^{(2)}(\boldsymbol{q}_{T}-\boldsymbol{\mathcal{P}}_{\perp})[S_{\bar{n}}^{{\dagger}}S_{n}](0){\rangle}\;,$
(7)

where $N_{R}=N_{c}$ for $S_{n/\bar{n}}$ in the fundamental and $N_{c}^{2}-1$
for the adjoint representation of $SU(N_{c})$. This function has been
calculated at NNLO in Echevarria:2015byo . This function is responsible for
the TMD evolution which is actually known up to third order Li:2016ctv ;
Vladimirov:2016dll . The power corrections to the evolution have been studied
in Scimemi:2016ffw . Because of the universality of this soft function the
non-perturbative corrections that it generates in the TMD-evolution factor are
process independent Collins:2011zzd ; GarciaEchevarria:2011rb ;
Scimemi:2016ffw .

The soft factor provides finally a rapidity renormalization factor for the
jets which is totally analogous to the TMD case, see ref. Echevarria:2016scs ,
so that in this sense we can re-write eq. (6) as

$\frac{d\sigma}{de_{1}de_{2}d\boldsymbol{q}_{T}}=H^{ij}_{2}(Q;\mu)\times\mathcal{J}^{\perp}_{i}(e_{1},Q,z_{\text{cut}},\boldsymbol{q}_{T};\mu,\zeta_{A})\otimes\mathcal{J}^{\perp}_{j}(e_{2},Q,z_{\text{cut}},\boldsymbol{q}_{T};\mu,\zeta_{B})\;,$
(8)

with $\zeta_{A}\zeta_{B}=Q^{4}z_{\mathrm{cut}}^{4}$ , which recalls clearly
the all-order factorization for the di-hadron fragmentation case using TMD.
The hadronization corrections to eq. (6-8) are discussed in more detail in
sec. 4.

The jet-TMD of eq. (8) can be re-factorized depending on the relative
magnitudes of the effective scales which define it so that one can identify
the more complete set of modes

soft: $\displaystyle\;\;\;p_{s}^{\mu}\sim q_{T}(1,1,1);$ collinear:
$\displaystyle\;\;\;p_{c}^{\mu}\sim
Q(\lambda_{c}^{2},1,\lambda_{c}),\;\lambda_{c}=\sqrt{e};$ soft-collinear:
$\displaystyle\;\;\;p_{sc}^{\mu}\sim
Qz_{\text{cut}}(\lambda_{sc}^{2},1,\lambda_{sc}),\;\lambda_{sc}=q_{T}/(Qz_{\text{cut}});$
collinear-soft: $\displaystyle\;\;\;p_{cs}^{\mu}\sim
Qz_{\text{cut}}(\lambda_{cs}^{2},1,\lambda_{cs}),\;\lambda_{cs}=\sqrt{e/z_{\text{cut}}}$
(9)

and we illustrate this in fig. 3. We start considering the limit $q_{T}\gtrsim
Q\sqrt{e}\gg Q\sqrt{ez_{\text{cut}}}$, which corresponds to region II in fig.
3, when the unintegrated and unsubtracted jet function,
$\mathcal{J}^{\perp}_{i}$, in eq. (6) can be re-factorized into three terms,

$\mathcal{J}^{\perp}_{i}(e,Q,z_{\text{cut}},\boldsymbol{q}_{T})=S_{sc,i}^{\perp}(Qz_{\text{cut}},\boldsymbol{q}_{T})\times\int
de^{\prime}\;S_{cs,i}(e-e^{\prime},Qz_{\text{cut}})J_{i}(e^{\prime},Q)$ (10)

where all the rapidity divergent part and transverse momentum dependence is
contained in the calculable $S_{sc,i}^{\perp}$. The subtracted and
unsubtracted jet-TMD are related by

$\displaystyle\mathcal{J}^{\perp}_{i}(e,Q,z_{\text{cut}},\boldsymbol{b},\mu,\zeta)=\sqrt{S(\boldsymbol{b})}\mathcal{J}^{\perp}_{i}(e,Q,z_{\text{cut}},\boldsymbol{b})$
(11)

where we have expressed all the subtraction in $\boldsymbol{b}$-space.
333Throughout the paper we will interchange between
$\boldsymbol{q}_{T},\boldsymbol{b}$ spaces for the transverse spectrum and
between $e,s$ spaces for the jet mass. We use the same symbol for any function
in either space. The variable we are working in should be clear from the
argument of the function.

Figure 3: Three possible hierarchies for $q_{T}$. Shaded region is one that
fails Soft-Drop. (I) Largest $q_{T}\sim Qz_{\text{cut}}$. The cross section is
factorized into 3 function s, cs and c. (II) The soft function s splits into
two s and sc.(III) The sc function merges with the cs function.

For smaller values of $q_{T}$: $Q\gg Qz_{\text{cut}}\gtrsim Q\sqrt{e}\gg
q_{T}\sim Q\sqrt{ez_{\text{cut}}}$, the collinear-soft and soft-collinear
merge into the same mode which we still refer to as collinear-soft. The soft
and collinear modes remain unchanged in their scaling compared to region II.
The form of factorization theorem in eq. (6) does not change but now the
corresponding jet TMDs are re-factorized as (see region III in fig. 3),

$\mathcal{J}^{\perp}_{i}(e,Q,z_{\text{cut}},\boldsymbol{q}_{T})=\int
de^{\prime}\;S_{cs,i}^{\perp}(e-e^{\prime},Qz_{\text{cut}},\boldsymbol{q}_{T})J_{i}(e^{\prime},Q).$
(12)

Several of the parameters in the differential cross-secion in eq. (8) are in
practice integrated in experiments, so that it is convenient to explicitly
write the cumulant (or partially integrated) distribution

$\frac{d\sigma}{d\boldsymbol{q}_{T}}(e_{\text{cut}})=\int_{0}^{e_{\text{cut}}}de_{1}de_{2}\frac{d\sigma}{de_{1}de_{2}d\boldsymbol{q}_{T}}\;.$
(13)

For this cross section we work with the integrated jet function which depends
on $e_{\text{cut}}$ rather than $e$,

$\mathcal{J}^{\perp}_{j}(e_{\text{cut}},Q,z_{\text{cut}},{\boldsymbol{q}}_{T};\mu,\zeta)=\int_{0}^{e_{\text{cut}}}de\;\mathcal{J}^{\perp}_{j}(e,Q,z_{\text{cut}},{\boldsymbol{q}}_{T};\mu,\zeta)\;.$
(14)

and the factorization theorem for electron-positron annihilation is

$\frac{d\sigma}{d\boldsymbol{q}_{T}}(e_{\text{cut}})=H^{ij}_{2}(Q;\mu)\int\frac{d\boldsymbol{b}}{4\pi}e^{i\boldsymbol{b}\cdot\boldsymbol{q}_{T}}\mathcal{J}^{\perp}_{i}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{b};\mu,\zeta)\mathcal{J}^{\perp}_{j}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{b};\mu,\zeta)\;.$
(15)

The resummation of logarithms inside the jet-TMD implied by eq. (10-12) is
taken into account defining the cumulant jet function as

$\displaystyle\mathcal{J}^{\perp}_{i}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{b};\mu,\zeta)$
$\displaystyle=\sqrt{S(\boldsymbol{b})}\mathcal{J}^{\perp}_{i}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{b})\
,$ (16)
$\displaystyle\mathcal{J}^{\perp}_{i}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{b})$
$\displaystyle=S_{sc,i}^{\perp}(Qz_{\text{cut}},\boldsymbol{b})\mathcal{J}_{i}(e_{\text{cut}},Q,z_{\text{cut}};\mu),$
(17) $\displaystyle\mathcal{J}_{i}(e_{\text{cut}},Q,z_{\text{cut}};\mu)$
$\displaystyle=\int_{0}^{e_{\text{cut}}}de\int
de^{\prime}\;S_{cs,i}(e-e^{\prime},Q,z_{\text{cut}};\mu)J_{i}(e^{\prime},Q;\mu)$
(18)

and we recall that the rapidity divergences are present only in $S$ and
$S_{sc,i}^{\perp}$, canceling in their product in eq. (16). With the exception
of the soft-collinear function, $S_{sc}^{\perp}$, all other ingredients of the
factorization are already known at least up to NLO accuracy. In app. B we
report the defining matrix elements of each function, we summarize the NLO
results and we perform the NLO calculation of $S_{sc}^{\perp}$. We have
performed the calculation using rapidity regulator. The connection between
rapidity regulator and $\zeta$-parameter is outlined in app. C.2.

Finally we observe that using monte-carlo simulations (particularly Pythia 8
Sjostrand:2006za ; Sjostrand:2007gs ) most of the events fall in the kinematic
regime

$Q\gg Qz_{\text{cut}}\gg q_{T}\sim Q\sqrt{e_{\text{cut}}}\;.$ (19)

An important consequence of the jet function refactorization in eq. (10) is
that the transverse momentum dependent elements decouple from the jet mass
elements. This suggests that, as long as we remain within the hierarchy of eq.
(19), then the exact mass cutoff on the invariant mass will only influence the
overall normalization and not the shape of the TMD distribution. We test this
observation against the monte-carlo simulations by comparing the normalized
TMD distributions for various values of $e_{\text{cut}}$. We show the results
in fig. 4 (left). The jet algorithm is implemented through FastJet-3
Cacciari:2011ma .

Figure 4: Left: The normalized cross sections for different values of the jet
mass cutoff parameter $e_{\text{cut}}$. We also include the corresponding
ratios with respect to the case $e_{\text{cut}}=0.01$. Right: The _relatively_
normalized cross section for fixed $e_{\text{cut}}=0.01$ and for different
value of the jet radius $R$. The corresponding ratios are with respect to
$R=1$.

In addition we note that as long as we measure $q_{T}\ll Qz_{\text{cut}}$ and
for $R\sim 1$ the shape and normalization of the cross section is independent
of the choice of $R$. We also demonstrate this with the help of simulations.
We simulate events at $Q=50$ GeV and we analyze them for different values of
$R\gtrsim 1$. We show the resulting distributions in fig. 4 (right). Note that
for that plot we preserve the relative normalizations of the curves.

###  2.3 Renormalization group evolution

The two main quantities involved in the factorization procedure carried out in
previous section are the subtracted jet-TMD for which we have

$\displaystyle\mu\frac{d}{d\mu}\mathcal{J}^{\perp}(e,Q,z_{\rm
cut},\boldsymbol{b},\mu,\zeta)=\gamma_{F}^{q}(\mu,\zeta)\mathcal{J}^{\perp}(e,Q,z_{\rm
cut},\boldsymbol{b},\mu,\zeta),$ (20)
$\displaystyle\zeta\frac{d}{d\zeta}\mathcal{J}^{\perp}(e,Q,z_{\rm
cut},\boldsymbol{b},\mu,\zeta)=-\mathcal{D}^{q}(\mu,\boldsymbol{b})\mathcal{J}^{\perp}(e,Q,z_{\rm
cut},\boldsymbol{b},\mu,\zeta),$ (21)

where on the r.h.s. we have considered just quark initiated jets and we have
Fourier transformed with respect to ${\boldsymbol{q}}_{T}$ the jet functions
appearing in eq. (8). Of course this result recalls literally the standard TMD
case.

However, because of the re-factorization of $\mathcal{J}^{\perp}$ (see eq.
(10-12)) this resummation is not complete and large logarithms can still spoil
the convergence of the perturbative series. Defining $s$ as the variable
conjugate to $e$ in Laplace space (see app. A) and

$\displaystyle G\in\left\\{S_{sc}^{sub}(Qz_{\rm
cut},\boldsymbol{b}),\,S_{cs}(s,Qz_{cut}),\,J(s,Q)\right\\};\quad
S_{sc}^{sub}(Qz_{\rm
cut},\boldsymbol{b})=\sqrt{S(\boldsymbol{b})}S_{sc}(Qz_{\rm
cut},\boldsymbol{b})\;,$ (22)

we have

$\mu\frac{d}{d\mu}G=\gamma^{G}(\mu,\alpha_{s})G=\left(\Gamma^{G}[\alpha_{S}]\textbf{l}_{m^{2}_{G}}+\Delta\gamma^{G}[\alpha_{S}]\right)G,$
(23)

which are formally similar to the TMD case and the values of $m_{G}$ are
reported in the appendix in tab. 1. The only function in $G$ which has a
rapidity evolution equation is $S_{sc}^{sub}$ and it scales like
$\mathcal{J}^{\perp}$ in eq. (21). The cusp part of eq. (23) is proportional
to the standard cusp anomalous dimension

$\Gamma^{G}_{\mu}[\alpha_{s}]=\frac{\Gamma^{G}_{0}}{\Gamma^{\text{cusp}}_{0}}\Gamma^{\text{cusp}}=\frac{\Gamma^{G}_{0}}{\Gamma^{\text{cusp}}_{0}}\sum_{n=0}^{\infty}\left(\frac{\alpha_{s}}{4\pi}\right)^{1+n}\Gamma^{\text{cusp}}_{n},$
(24)

For the non-cusp part we have also a perturbative expansion

$\Delta\gamma^{G}[\alpha_{S}]=\sum_{n=0}^{\infty}\left(\frac{\alpha_{s}}{4\pi}\right)^{1+n}\gamma^{G}_{n}.$
(25)

The anomalous dimensions that enter in the calculations for each case are
given in app. C. The evolution in rapidity and factorization scales of all
quantities can be implemented using the the $\zeta$-prescription whose general
framework can be found in ref. Scimemi:2018xaf . We provide some details for
the present case in the appendix.

The resummation of potentially large logarithms inside the jet-TMD is done
performing the evolution in Laplace space and then integrating such that we
get the cumulant before we take the inverse transform. In this way we resum
logarithms which are associated to $e_{\text{cut}}$. All this works as
follows. Starting from eq. (18), then taking the Laplace and consecutively the
inverse transform with respect to $e$ we find

$\mathcal{J}_{i}(e_{\text{cut}},Q,z_{\text{cut}};\mu)=\frac{1}{2\pi
i}\int_{\gamma-i\infty}^{\gamma+i\infty}ds\frac{\exp(se_{\text{cut}})-1}{s}S_{cs,i}(s,Q,z_{\text{cut}};\mu)J_{i}(s,Q;\mu)\
.$ (26)

Then solving the RGE equations for the collinear-soft and jet function as
described in app. C.1, and performing the last remaining integral over the
Laplace conjugate variable $s$ we get

$\mathcal{J}_{i}(e_{\text{cut}},Q,z_{\text{cut}};\mu)=\exp\Big{(}K_{cs}(\mu,\mu_{cs})+K_{J}(\mu,\mu_{J})\Big{)}S_{cs,i}(L_{cs}\to\partial_{\omega_{cs}};\mu_{cs})J_{i}(L_{J}\to\partial_{\omega_{J}};\mu_{J})\\\
\Big{(}\frac{\mu_{cs}}{Q\sqrt{z_{\text{cut}}e_{\text{cut}}}}\Big{)}^{2\omega_{cs}(\mu,\mu_{cs})}\Big{(}\frac{\mu_{J}}{Q\sqrt{e_{\text{cut}}}}\Big{)}^{2\omega_{J}(\mu,\mu_{J})}\frac{\exp(\gamma_{E}(\omega_{cs}(\mu,\mu_{cs})+\omega_{J}(\mu,\mu_{J})))}{\Gamma(1-\omega_{cs}(\mu,\mu_{cs})-\omega_{J}(\mu,\mu_{J}))}\
.$ (27)

This is our final result for the resummed cumulant jet function. The order of
logarithmic accuracy is then determined by the order of which the kernels
$K_{F}$, $\omega_{F}$, and the fixed order collinear-soft and jet functions
are evaluated. At this stage of the calculation the canonical scales,
$\mu_{cs}$ and $\mu_{J}$, are not yet fixed. This allows us to choose the
scales such that potentially large logarithms are minimized in momentum space.
From the above is clear that the canonical choice of scales such as the fixed
order logarithms are minimized are,

$\displaystyle\mu_{cs}$
$\displaystyle=Q\sqrt{z_{\text{cut}}e_{\text{cut}}}\;,$ $\displaystyle\mu_{J}$
$\displaystyle=Q\sqrt{e_{\text{cut}}}\;.$ (28)

In numerical applications one needs to perform variations around these scales
in order to obtain an estimate of the theoretical uncertainty.

###  2.4 Numerical results for $e^{+}e^{-}$

In this section, we provide the results of our calculation for
$e^{+}e^{-}\to\text{2\;jets}$ computed up to NNLL accuracy. The implementation
necessarily needs a choice for the rapidity scales and we have done it using
the $\zeta$-prescription as described in ref. Scimemi:2018xaf  and adapting
the code artemide to the present case. This consisted of performing the
evolution of the transverse momentum dependent components within the artemide
framework, while for all other scales not involved in the rapidity evolution,
i.e., the hard and jet functions, see app. C.1.

There are some important modifications to the $\zeta$-prescription framework
for our case which affect the numerics. One of this is that now
$\zeta_{A}\zeta_{B}\sim Q^{4}z_{\text{cut}}^{4}$ compared to the di-hadron
decorrelation case where $\zeta_{A}\zeta_{B}\sim Q^{4}$. This means that the
effective hard scale to which the distributions are sensitive is lower.
Because the TMD factorization is valid when $q_{T}$ is much lower than the
hard scale of the process, one needs that the product $Qz_{\mathrm{cut}}$ be
sufficiently high. In our plots we have considered the case $q_{T}\lesssim
Qz_{\mathrm{cut}}$. Then the evolution of the jet-TMD given in eq. (23) is
also slightly different from the standard hadron TMD, although the changes are
implemented easily in the artemide code. A one-loop check of all anomalous
dimensions is provided in app. B.

Figure 5: Comparison of the NLL result against the partonic shower of Pythia 8
for $R=1$ and $e_{\text{cut}}=0.01$ for two different center of mass energies,
Left: 50 GeV, Right: 100 GeV.

In fig. 5 we compare our analytic result for NLL cross section (normalized)
against Pythia simulations for $Q=50$ and 100 GeV. For the purposes of
comparison we turn hadronization off in the simulation and we compare against
our purely perturbative result. The perturbative calculation depends on the
parameter $B_{\text{NP}}$ which in practice implements a cutoff in the inverse
Laplace transform such that the soft scale, that behaves as $1/b$, does not
hit the Landau pole. As long as we choose this parameter such that convergence
of the integral is reached before the cutoff, then the perturbative result is
not much sensitive to the value of $B_{\text{NP}}$. Although, as we now
discuss, the theoretical uncertainty of the cross section for these energies
at NLL is quite large, we find very good agreement with the simulations for
the canonical choice of scales (i.e., central line in fig. 5).

In fig. 6 we give the NNLL results including a theoretical uncertainty band.
We compare against the NLL cross section and although the error bands seem to
be larger than what is typically expected we can clearly see that the result
convergences and the theory error decreases by approximately factor of two. To
estimate the theoretical uncertainty we first vary all the factorization
scales of a factor 2 (0.5) around their canonical value, then we separately
take the envelope of the variations involved in rapidity evolution,
$\mu,\mu_{sc}$, and of the ones involved only in the virtuality evolution of
the jet function, $\mu_{cs},\mu_{J}$. The final error bands we show are the
quadrature of the two contributions. The reason for this prescription is that
rapidity and virtuality evolutions are in principle uncorrelated. The
uncertainty is somewhat larger than what one might expect for a NNLL
calculation, and is practically dominated by the variations in the jet
function. This is attributed to the small values of the collinear-soft scale,
$\mu_{cs}\sim Q\sqrt{e_{\text{cut}}z_{\text{cut}}}$, which approaches the non-
perturbative regime even for values of $Q\sim m_{Z}$. One might attempt to
reduce the uncertainty by increasing either $e_{\text{cut}}$ or
$z_{\text{cut}}$, but caution is needed not to invalidate the corresponding
hierarchy. We will see later that when only the mass of one jet is measured
(e.g., in DIS or hadron-jet decorrelation) then the error band decreases
significantly.

Figure 6: Transverse momentum de-correlation for $e^{+}e^{-}\rightarrow$
dijets with center of mass energy at the Z mass.

##  3 Jets in DIS

The advent of new colliders like EIC and LHeC makes the measurement of jets
interesting also in semi-inclusive deep inelastic scattering (SIDIS)
experiments. Actually we want to explore the possibility of using jets to
study the TMDPDF.

For the present case we demand that the hard scattering of the lepton on the
proton produces a single jet. In the Breit frame we measure the transverse
component, $\boldsymbol{q}_{T}$, of the transferred momentum,
$q^{\mu}={k^{\prime}}^{\mu}-{k}^{\mu}$ with respect to the single groomed jet.
As before, we impose a jet mass cut-off $e_{\text{cut}}$ and the grooming
parameter $z_{\text{cut}}$. In this framework the initial state proton is
moving along the $-z$ direction and the final state jet is moving in the
opposite $+z$-direction, so that we can assign the directions $n$ and
$\bar{n}$ to the beam and jet definition. The contribution to this transverse
momentum measurement comes from the initial state radiation which forms part
of the TMDPDF and the radiation that fails soft-drop in the final state jet.
We demand that there is a single energetic jet with $E_{J}\sim
Q/2=\sqrt{-q^{2}}/2$ with accompanying soft radiation.

It is instructive to setup some of the notation that we are using for
describing the kinematics in the Breit frame. The virtual photon is assumed to
be completely space-like and it has only the $z$ component of the momentum.
Defining our light-cone directions $n^{\mu}=(1,0,0,+1)$ and
$\bar{n}^{\mu}=(1,0,0,-1)$, the resulting photon momentum $q^{\mu}$ is

$q^{\mu}=\frac{Q}{2}(n^{\mu}-\bar{n}^{\mu})\ ,$ (29)

where $Q^{2}=-q^{2}$ is a positive quantity. We assume that at the partonic
level, a single quark carrying $x$ fraction of the proton longitudinal
momentum undergoes a hard interaction with the virtual photon. In this frame,
the proton is moving along the $-z$ direction and its momentum can be written
as:

$P^{\mu}=\frac{Q}{2x}\bar{n}^{\mu}\ .$ (30)

At tree level and by momentum conservation the final state parton will carry
momentum

$xP^{\mu}+q^{\mu}=\frac{Q}{2}(n^{\mu}-\bar{n}^{\mu})+\frac{Q}{2}\bar{n}^{\mu}=\frac{Q}{2}n^{\mu}\;,$
(31)

which is exactly opposite in direction to the incoming beam. Of course this
will be modified beyond tree-level when initial and final state radiation is
included.

###  3.1 Schematics for factorization

Since we are working with two back-to-back directions, our usual definition of
the soft function holds: in other words the change from future pointing to
past pointing Wilson lines does not affect its value Collins:2011zzd ;
GarciaEchevarria:2011rb ; Echevarria:2014rua ; Gaunt:2014ska ;
Vladimirov:2017ksc .

Since we still impose the same jet mass measurement on the final state jet, we
have all the modes that we had in the $e^{+}e^{-}$ case. The main difference
is that now the initial hadronic state is a TMDPDF. The form of the factorized
cross section follows again the hierarchy $Q\gg Qz_{\text{cut}}\gg
q_{T},\;R\sim 1$ and

$\frac{d\sigma}{dxdQ^{2}d\boldsymbol{q}_{T}}=\mathcal{N}(x,Q)H_{2}(Q,\mu)\times
S(\boldsymbol{q}_{T})\otimes B_{i\leftarrow
h}(x,Q,\boldsymbol{q}_{T})\otimes\mathcal{J}_{j}^{\perp}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{q}_{T})\
,$ (32)

where $x=-q^{2}/(2P\cdot k)$, $k$ is the momentum of the incoming electron,
and $\mathcal{N}(x,Q)$ is the over-all normalization which we give later in
this section. The un-subtracted TMDPDF is $B_{i\leftarrow h}$. In our rapidity
regularization scheme the (subtracted) TMDPDF is defined as

$\displaystyle F_{i\leftarrow
h}(x,\boldsymbol{b};\mu,\zeta)=\sqrt{S(\boldsymbol{b})}B_{i\leftarrow
h}(x,Q,\boldsymbol{b}).$ (33)

At perturbative values of $q_{T}$, the $F_{i\leftarrow h}$ can be matched onto
the collinear PDF. The matching coefficients at NNLO are evaluated in
Gehrmann:2014yya ; Echevarria:2016scs  and in the appendix we review some one-
loop results. Once the subtracted quantities are included we can write

$\frac{d\sigma}{dxdQ^{2}d\boldsymbol{q}_{T}}=\mathcal{N}(x,Q)H_{2}(Q,\mu)\int\frac{d\boldsymbol{b}}{4\pi^{2}}e^{i\boldsymbol{b}\cdot\boldsymbol{q}_{T}}F_{i\leftarrow
h}(x,Q,\boldsymbol{b},\mu,\zeta_{A})\mathcal{J}_{j}^{\perp}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{b};\mu,\zeta_{B})\
.$ (34)

The evolution under renormalization group equations for the TMDPDF is widely
known (see e.g. Scimemi:2018xaf ; DAlesio:2014mrz ; Echevarria:2012pw ;
Chiu:2011qc ) and we recall a few characteristics here. One has

$\displaystyle\mu\frac{d}{d\mu}F_{f\leftarrow
f^{\prime}}(x,\boldsymbol{b},\mu,\zeta)=\gamma_{F}^{f}(\mu,\zeta)F_{f\leftarrow
f^{\prime}}(x,\boldsymbol{b},\mu,\zeta),$
$\displaystyle\zeta\frac{d}{d\zeta}F_{f\leftarrow
f^{\prime}}(x,\boldsymbol{b},\mu,\zeta)=-\mathcal{D}^{f}(\mu,\boldsymbol{b})F_{f\leftarrow
f^{\prime}}(x,\boldsymbol{b},\mu,\zeta),$ (35)

where $\mathcal{D}_{f}$ and $\gamma_{F}^{f}$ are the rapidity and UV anomalous
dimensions, respectively. The integrability requirement of this couple of
equation results in

$\displaystyle\mu\frac{d}{d\mu}\left(-\mathcal{D}^{f}(\mu,\boldsymbol{b})\right)=\zeta\frac{d}{d\zeta}\gamma^{f}_{F}(\mu,\zeta)=-\Gamma^{\rm
cusp}_{f}$ (36)

where $\Gamma^{\rm cusp}_{f}$ is the cusp anomalous dimension. The UV
anomalous dimension is written in these terms as

$\displaystyle\gamma_{F}^{f}=\Gamma^{\rm
cusp}_{f}\textbf{l}_{\zeta}-\gamma_{V}^{f},$ (37)

$\gamma_{V}^{f}$ being the non-cusp part of the anomalous dimension and
$\textbf{l}_{\zeta}=\ln\left(\mu^{2}/\zeta\right)$. The $\gamma_{V}$ and
$\mathcal{D}$ anomalous dimensions are known up to $\mathcal{O}(a_{s}^{3})$
Moch:2004pa ; Moch:2005tm ; Baikov:2009bg ; Vladimirov:2016dll ; Li:2016ctv .
A numerical calculation for the four-loop cusp anomalous dimension was
recently given in Vogt:2018miu . All the evolution equations are the same for
the case of TMD fragmentation functions, and we do not discuss them any more
here.

###  3.2 Derivation of the factorized cross section using jets

In this section we provide some details for the factorization of the SIDIS
cross section in eq. (32, 34). The scattering amplitude for the process $ep\to
ef$ where $f$ is the final state is given by:

$iM(ep\to
ef)=(-ie^{2})\bar{u}(k^{\prime})\gamma_{\mu}u(k)\frac{1}{q^{2}}{\langle}f{|}J^{\mu}(0){|}p(P){\rangle}\;,$
(38)

and thus the corresponding cross section is given by

$d\sigma(ep\to
ef)=\frac{e^{4}}{4(s-m^{2})}\int\frac{d^{3}k^{\prime}}{2(2\pi)^{3}E_{k^{\prime}}}\;\text{tr}\Big{[}\not{k}\gamma_{\mu}\not{k^{\prime}}\gamma_{\nu}\Big{]}\\\
\sum_{f}\int
d\Pi_{f}\;{\langle}p(P){|}J^{{\dagger}\mu}(0){|}f{\rangle}{\langle}f{|}J^{\nu}(0){|}p(P){\rangle}(2\pi)^{4}\delta^{(4)}(q+P-p_{f})\;,$
(39)

where $q=k^{\prime}-k$. We can use the standard parametrization of the final
electron phase-space to write:

$\int\frac{d^{3}k^{\prime}}{2(2\pi)^{3}E_{k^{\prime}}}=dxdy\frac{ys}{(4\pi)^{2}}\;,$
(40)

where $y=(2P\cdot q)/(2P\cdot k)$ and $s$ is the hadronic Mandelstam variable.
We then get,

$\frac{d\sigma}{dxdy}(ep\to ef)=L_{\mu\nu}(k,k^{\prime})\sum_{f}\int
d^{4}re^{iq\cdot r}\int
d\Pi_{f}\;{\langle}p(P){|}J^{{\dagger}\mu}(0){|}f{\rangle}{\langle}f{|}J^{\nu}(x){|}p(P){\rangle}\;,$
(41)

where $r^{\mu}$ is Fourier conjugate of the momenta $q^{\mu}$ and $L^{\mu\nu}$
is the leptonic tensor,

$L_{\mu\nu}(k,k^{\prime})\equiv\frac{\alpha^{2}ys}{4(s-m^{2})}\text{tr}\Big{[}\not{k}\gamma_{\mu}\not{k^{\prime}}\gamma_{\nu}\Big{]}\;.$
(42)

The next step is to project the hadronic final state $|f\rangle$ onto the one
that corresponds to the measurement that we are proposing, i.e.,

$\int d\Pi_{f}{|}f{\rangle}{\langle}f{|}\;\;\to\int d\boldsymbol{q}_{T}zdz\int
d\Pi_{f[\text{g-jet}(z\boldsymbol{q}_{T},z)]}{|}f{\rangle}{\langle}f{|}\;.$
(43)

We can now match the full theory hadronic current $J_{\mu}(x)$ onto the SCET+
Bauer:2011uc  current working in the Breit frame,

$J^{\mu}(x)=C^{\mu\nu}(Q)\Big{[}\bar{\chi}_{n,Q}S_{n}^{{\dagger}}W_{t}^{{\dagger}}U_{n}\gamma_{\nu}S_{\bar{n}}\chi_{\bar{n},Q}\Big{]}+\mathcal{O}(\lambda)\;,$
(44)

where $\lambda$ is the power counting parameter of our EFT which will turn out
to be $q_{\perp}/Q\sim e_{\text{cut}}/Q$. Note that in the same step, through
BPS field redefinition, we decoupled the collinear soft modes from the
collinear modes and hence the presence of the $U_{n}$ Wilson lines. In the
matching we also have the soft Wilson lines $S_{n}$. From the kinematic
constraints of the measurement and since all the modes that are present in the
projected final state are decoupled from each other at the level of the
Lagrangian, (we assume that contributions from Glauber gluon exchanges cancel)
it is possible to factorize the final state as follows,

${|}f{\rangle}\;\;\to\;\;{|}X_{\bar{n}}{\rangle}{|}X_{n}{\rangle}{|}X_{s}{\rangle}{|}X_{sc}{\rangle}\;,$
(45)

where we have included in $X_{n}$ all possible modes that contribute to the
invariant mass measurement. Refactorization of the $n-$collinear sector
follows from the same steps as in the case of electron-positron annihilation
presented in ref. Frye:2016aiz . We are now ready to factorize the cross
section into individual SCET matrix elements. In the final result one needs to
be careful regarding all the index contractions and the tensor structures.
This was carefully considered in ref. Kang:2013nha . In addition we are
considering the case where the frame we are working is rotated such that the
transverse momentum of the groomed jet is zero. After all rearrangements we
get,

$\frac{d\sigma}{dxdydzd\boldsymbol{q}_{T}}(ep\to ef)=\sigma_{0}(x,Q)\times
H_{2}(Q)\int d^{4}re^{iq\cdot r}\\\
\frac{1}{N_{c}}\sum_{X_{s}}{\langle}0{|}S_{n}S_{\bar{n}}^{{\dagger}}(r_{\perp}){|}X_{s}{\rangle}{\langle}X_{s}{|}S_{\bar{n}}S_{n}^{{\dagger}}(0){|}0{\rangle}\\\
\times\sum_{X_{\bar{n}}}{\langle}p(P){|}\bar{\chi}_{\bar{n}}(r^{+},r_{\perp})\frac{\gamma^{+}}{2}{|}X_{\bar{n}}{\rangle}{\langle}X_{\bar{n}}{|}\chi_{\bar{n}}(0){|}p(P){\rangle}\\\
\times\frac{1}{N_{c}}\sum_{X_{sc}}{\langle}0{|}U_{n}^{{\dagger}}W_{t}(r_{\perp}){|}X_{sc}{\rangle}{\langle}X_{sc}{|}W_{t}^{{\dagger}}U_{n}(0){|}0{\rangle}\\\
\times\frac{z}{2N_{c}}\text{tr}\sum_{X_{n}}{\langle}0{|}\frac{\gamma^{-}}{2}\chi_{n}(r^{-},r_{\perp}){|}X_{n}{\rangle}{\langle}X_{n}{|}\bar{\chi}_{n}(0){|}0{\rangle}{|}_{p_{\perp}^{X_{n}}=0}\;.$
(46)

The hard matching coefficient in general has two Lorentz structures, given the
two types of currents, vector and axial. For the case of photon with vector
current, we simply have $H^{\mu\nu}\sim g^{\mu\nu}_{\perp}$. We have also
multipole-expanded the final result. To proceed with the factorization theorem
in momentum space, we remove $r_{\perp}$ dependence from the various EFT
matrix elements by acting the corresponding fields on the final states. This
gives us

$\frac{d\sigma}{dxdydzd\boldsymbol{q}_{T}}(ep\to ef)=\sigma_{0}(x,Q)\times
H_{2}(Q)\int d^{4}re^{iq\cdot
r+i(\boldsymbol{p}^{X^{R}_{\bar{n}}}_{\perp}+\boldsymbol{p}_{\perp}^{S})\cdot\boldsymbol{r}_{\perp}}\\\
\frac{1}{N_{c}}\sum_{X_{s}}{\langle}0{|}S_{n}S_{\bar{n}}^{{\dagger}}(0){|}X_{s}{\rangle}{\langle}X_{s}{|}S_{\bar{n}}S_{n}^{{\dagger}}(0){|}0{\rangle}\\\
\times\sum_{X_{\bar{n}}}{\langle}p(P){|}\bar{\chi}_{\bar{n}}(r^{+},0_{\perp})\frac{\gamma^{+}}{2}{|}X_{\bar{n}}{\rangle}{\langle}X_{\bar{n}}{|}\chi_{\bar{n}}(0){|}p(P){\rangle}\\\
\times\frac{1}{N_{c}}\sum_{X_{sc}}{\langle}0{|}U_{n}^{{\dagger}}W_{t}(0){|}X_{sc}{\rangle}{\langle}X_{sc}{|}W_{t}^{{\dagger}}U_{n}(0){|}0{\rangle}\\\
\times\frac{z}{2N_{c}}\;\text{tr}\sum_{X_{n}}{\langle}0{|}\frac{\gamma^{-}}{2}\chi_{n}(r^{-},0_{\perp}){|}X_{n}{\rangle}{\langle}X_{n}{|}\bar{\chi}_{n}(0){|}0{\rangle}{|}_{p_{\perp}^{X_{n}}=0}\;,$
(47)

where

$p_{\perp}^{X_{\bar{n}}^{R}}{|}_{p_{\perp}^{\text{g-jet}}=0}=p_{\perp}^{X_{\bar{n}}}-P_{\perp}{|}_{p_{\perp}^{\text{g-jet}}=0}=p_{\perp}^{X_{\bar{n}}}{|}_{P_{\perp}=0}\Big{(}1+\mathcal{O}(\lambda)\Big{)}\;,$
(48)

is the difference in the transverse momentum of the recoiling initial state
collinear radiation and the proton with respect to the hadrons direction,
which up to power-corrections of order $\mathcal{O}(\lambda)$ is simply the
transverse momentum of the recoiling radiation with respect to the proton.
Performing the integral over $d^{4}r$ we get:

$\frac{d\sigma}{dxdydzd\boldsymbol{q}_{T}}(ep\to ef)=\sigma_{0}(x,Q)\times
H_{2}(Q)\delta^{(2)}(\boldsymbol{q}_{T}+\boldsymbol{p}^{X^{R}_{\bar{n}}}_{\perp}+\boldsymbol{p}_{\perp}^{X_{s}}+\boldsymbol{p}_{\perp}^{X_{sc}})\\\
\frac{1}{N_{c}}\sum_{X_{s}}{\langle}0{|}S_{n}S_{\bar{n}}^{{\dagger}}(0){|}X_{s}{\rangle}{\langle}X_{s}{|}S_{\bar{n}}S_{n}^{{\dagger}}(0){|}0{\rangle}\\\
\times\sum_{X_{\bar{n}}}{\langle}p(P){|}\bar{\chi}_{\bar{n}}(0)\frac{\gamma^{+}}{2}\delta(q^{-}-p_{X_{\bar{n}}}^{-}){|}X_{\bar{n}}{\rangle}{\langle}X_{\bar{n}}{|}\chi_{\bar{n}}(0){|}p(P){\rangle}\\\
\times\frac{1}{N_{c}}\sum_{X_{sc}}{\langle}0{|}U_{n}^{{\dagger}}W_{t}(0){|}X_{sc}{\rangle}{\langle}X_{sc}{|}W_{t}^{{\dagger}}U_{n}(0){|}0{\rangle}\\\
\times\frac{z}{2N_{c}}\;\text{tr}\sum_{X_{n}}{\langle}0{|}\frac{\gamma^{-}}{2}\chi_{n}(0)\delta(q^{+}-p_{X_{n}}^{+}){|}X_{n}{\rangle}{\langle}X_{n}{|}\bar{\chi}_{n}(0){|}0{\rangle}{|}_{p_{\perp}^{X_{n}}=0}\;.$
(49)

In order to simplify our result further we introduce “measurement” delta
functions for the soft and initial state matrix elements. This will allow us
to absorb the $p_{\perp}^{X_{i}}$ into the corresponding matrix elements and
use

${\boldsymbol{1}}_{i}=\sum_{X_{i}}{|}X_{i}{\rangle}{\langle}X_{i}{|}\;,$ (50)

to further simplify the form of EFT matrix elements. We also perform a type-I
RPI transformation in order to rewrite the proton matrix elements as function
of fields with respect to the initial state proton axis. We thus get

$\frac{d\sigma}{dxdydzd\boldsymbol{q}_{T}}(ep\to ef)=\sigma_{0}(x,Q)\times
H_{2}(Q)\int
d\boldsymbol{p}_{\perp}^{\;s}d\boldsymbol{p}_{\perp}^{\;sc}d\boldsymbol{p}_{\perp}^{\;c}\;\;\delta^{(2)}(\boldsymbol{q}_{T}+\boldsymbol{p}^{c}_{\perp}+\boldsymbol{p}_{\perp}^{s}+\boldsymbol{p}_{\perp}^{sc})\\\
\frac{1}{N_{c}}{\langle}0{|}T\Big{(}S_{n}S_{\bar{n}}^{{\dagger}}(0)\Big{)}\delta^{(2)}(\boldsymbol{p}_{\perp}^{s}-\boldsymbol{\mathcal{P}}_{\perp})\bar{T}\Big{(}S_{\bar{n}}S_{n}^{{\dagger}}(0)\Big{)}{|}0{\rangle}\\\
\times{\langle}p(P){|}\bar{\chi}_{\bar{n}}(0)\frac{\gamma^{+}}{2}\delta(q^{-}-\mathcal{P}^{-})\delta^{(2)}(\boldsymbol{p}_{\perp}^{c}-\boldsymbol{\mathcal{P}}_{\perp})\chi_{\bar{n}}(0){|}p(P){\rangle}{|}_{P_{\perp}=0}\\\
\times\frac{1}{N_{c}}{\langle}0{|}T\Big{(}U_{n}^{{\dagger}}W_{t}(0)\Big{)}\mathcal{M}_{\perp}^{SD}\bar{T}\Big{(}W_{t}^{{\dagger}}U_{n}(0)\Big{)}{|}0{\rangle}\\\
\times\frac{z}{2N_{c}}\;\text{tr}\sum_{X_{n}}{\langle}0{|}\frac{\gamma^{-}}{2}\chi_{n}(0)\delta(q^{+}-\mathcal{P}^{+}){|}X_{n}{\rangle}{\langle}X_{n}{|}\bar{\chi}_{n}(0){|}0{\rangle}{|}_{p_{\perp}^{X_{n}}=0}\;,$
(51)

where $\mathcal{M}_{\perp}^{SD}$ is the measurement function given in eq.
(99). Since we are considering only large radius jets with $R\gtrsim 1$ we may
trivially perform the integration of the energy fraction $z$ using $Q\simeq
p^{X_{n}}_{+}$ up to power corrections. Also performing change of integration
variables,

$dxdy=\frac{dxdQ^{2}}{xs}\;,$ (52)

we get eq. (32) with

$\mathcal{N}(x,Q)=\frac{\sigma_{0}(x,Q)}{xs}\;,$ (53)

and the matrix elements involved in the functions $S$, $B$, and $\mathcal{J}$
are given in the appendix. For the case of groomed jets with invariant mass
cutoff it is possible to refactorize the jet function. This is done in ref.
Frye:2016aiz  and thus we do not demonstrate it here. Then integrating over
$e\in(0,e_{\text{cut}})$ gives the dependence of the jet function in the
parameter $e_{\text{cut}}$. This is identical to the analysis in the previous
section on $e^{+}e^{-}$. This is our final result for the factorization
theorem in DIS.

###  3.3 Numerical results for DIS

In this section we use the factorization theorem in eq. (32) to obtain
numerical results for the TMD spectrum of groomed jets in DIS process. Our
analysis is done for two center-of-mass energies, EIC: $\sqrt{s}=100$ GeV and
HERA: $318$ GeV. For both energies we integrate over $y=Q^{2}/(xs)$ and
$Q=\sqrt{-q^{2}}$ in the regions $0.01<y<0.95$ and $40<Q<50$ GeV. For the
TMDPDFs we use the fits obtained from Drell-Yan data Bertone:2019nxa  with the
use of $\zeta$-prescription. In fig. 7 we show our results for NLL and NNLL
accuracies for the two center of mass choices, including theoretical
uncertainties. We estimate the theoretical scale variations as described in
sec. 2.4. The groomed jet parameters that we choose are the same as in the di-
lepton case: $\beta=0$, $z_{\text{cut}}=0.2$, and $e_{\text{cut}}=0.01$. As
before we find good convergence between the NLL and NNLL result. The absolute
value of theoretical scale variation is improvable with higher logarithmic
accuracy (NNLL-prime or perhaps N3LL), which needs the explicit calculation of
several jet hadronic matrix elements at two loops.

Figure 7: The NLL and NNLL TMD spectra for groomed jets in DIS for EIC (left:
$\sqrt{100}$ GeV) and HERA (right: $\sqrt{s}=318$ GeV) kinematics. The cross
section are integrated in $y=Q^{2}/(xs)$ and $Q=\sqrt{-q^{2}}$ (see details in
the main text).

We further investigate the size of the uncertainty due to the hadronic initial
state and the non-perturbative effects induced by TMD evolution. We do that by
varying the model parameters as constrained by the phenomenological analysis
in ref. Bertone:2019nxa  for our NNLL result. The results are shown in figure
8. We consider both variable and fixed $B_{\text{NP}}=2.5$ GeV-1 (for details
on the difference of the two schemes see Bertone:2019nxa ). We find that the
effects (for our kinematics) are particularly small, of the order of $\sim
5\%$, which is much smaller than the theoretical uncertainties. This suggests
that we need a better control over the theoretical uncertainties in order to
further constrain TMD distributions from groomed jets in DIS. As mentioned
earlier the uncertainty can be mitigated with higher logarithmic accuracy or
by choosing larger values of $e_{\text{cut}}$, still compatible with
factorization. This, will require to treat the region III shown in fig. 3. For
this reason it is interesting to investigate the range of values of
$e_{\text{cut}}$ for which the energetic wide angle radiation is avoided.

Figure 8: The NNLL cross-section including modeling of the initial hadronic
state effects fitted from Derll-Yan processes using two different scenes:
fixed and variable $B_{\text{NP}}$.

##  4 Hadronization effects

One of the goals of the paper is to study the non-perturbative effects
associated with TMD distributions, in this case the TMDPDF. Usually in any
experiment, there are multiple sources of non-perturbative corrections
associated with both the initial and final states. To have access to a
specific source of corrections, its therefore necessary to separate out the
pieces of interest from the uninteresting ones, which in this case constitute
the final state hadronization corrections. To access the TMD then, we must
already have a good extraction of the rest of the non-perturbative effects.
This is the reason why we consider distinct experiments in this paper. The
idea, as we shall demonstrate, is that the final state hadronization
corrections are exactly the same in the two experiments. The
$e^{+}e^{-}\to\text{2\;jets}$ case can be used to extract out all the final
state hadronization corrections, which can then be used for DIS.

For the $e^{+}e^{-}$ observable, the factorization takes the form in eq. (8),
we can then study the non-perturbative corrections for each collinear object
$\mathcal{J}^{\perp}_{i}$, which by symmetry, are the same for the two
objects. If we now look at the factorization for DIS, eq. (34), the key point
to note is that
$\mathcal{J}_{j}^{\perp}(e_{\text{cut}},Q,z_{\text{cut}},\boldsymbol{b};\mu,\zeta_{B})$
is the same object that appears in the the case of $e^{+}e^{-}$, while
$F_{i\leftarrow h}$ is just the TMDPDF. Thus it now becomes possible to
exclusively access the complete TMDPDF. We now wish to systematically list the
sources of the non-perturbative corrections associated with each factorized
function that appear in our cross section.

In order to use jets it is important to consider all the non-perturbative
effects for the case of our observables and in particular the ones coming from
the implementation of (groomed) jets. In fig. 1 we have shown that such
corrections are expected to be particularly small and we provide here a
discussion about their origin from a theory perspective. We have two
measurements on the jet: the jet mass, which is ultimately integrated over
some interval and acts as a normalization, and the transverse momentum
($p_{\perp}$) of the radiation that is groomed away. Since we are interested
in the shape of the $q_{T}$ spectrum, we will only consider the non-
perturbative effects in cross sections sensitive to it. As was explained in
sec. 2, we are working in the region II of EFT and we are going to discuss how
non-perturbative effect arise when we increase the value of $q_{T}$ (that is,
we discuss here the non-perturbative corrections in the small-$b$ limit, where
$b\equiv|\boldsymbol{b}|$). Our factorization theorem has four functions in
the IR, the collinear, the global soft, the collinear-soft, the soft-collinear
functions, see eq. (10-11), and all of them can potentially contribute to non-
perturbative power corrections. Even though the collinear and collinear-soft
functions do not contribute to $\boldsymbol{q}_{T}$ perturbatively, they can
still give a non-perturbative power correction to the ${\boldsymbol{q}}_{T}$
spectrum444There are also power corrections of similar magnitude in this
region due to the factorization of the $sc$ function from the $cs$, but they
are perturbative in nature and can be handled by making a smooth transition to
region III..

There are two types of non-perturbative corrections that we will consider
here. We call _shift_ non-perturbative effects the ones which are not altered
by the pass and fail procedure of the grooming conditions. An example is the
global soft function that is independent of the grooming procedure and it is
common to other TMD analysis. We refer to this kind of correction as _shift_
non-perturbative effects since, as we will see later, in the simplest case it
generates a shift in the TMD spectrum. The second correction instead is
related to the grooming procedure with $cs$ and $sc$ soft functions and the
jet shape function. In this case non-perturbative effects are driven by the so
called “non-perturbative particles” and it is obviously only possible when
perturbative modes are on the boundary of passing and failing soft-drop. We
refer to these contributions as _boundary_ non-perturbative effects.

###  4.1 _Shift_ non-perturbative correction

For the case of shift correction, we assume that the soft-drop condition
remains unaltered by any non-perturbative emissions. Now consider the
contribution to the shift correction by each function in turn.

The non-perturbative part of the global soft function defined in eq. (7) has
been studied in the literature in several frameworks Beneke:1995pq ;
Korchemsky:1994is ; Korchemsky:1997sy ; Beneke:1997sr ; Becher:2013iya ;
Scimemi:2016ffw . Up to ${\cal O}(b^{4})$ terms it can be written as

$\displaystyle\langle
0|T[S_{n}S^{\dagger}_{\bar{n}}(\boldsymbol{b})]\bar{T}[S_{\bar{n}}S_{n}^{\dagger}(0)]|0\rangle=\tilde{S}(b)+b^{2}\;\bar{C}^{(s)}_{i}(b)\langle
0|O^{i}|0\rangle\;,$ (54)

where $O^{i}$ is the complete set of local operators that have the same
quantum numbers as the soft function. Summation over $i$ is implied. Here
$\tilde{S}$ is the perturbative calculable part of the soft function and it
contains rapidity and UV divergences as well as the rest of other terms in the
equation. We can pull this out as a common factor to write

$\displaystyle\langle
0|T[S_{n}S^{\dagger}_{\bar{n}}(\boldsymbol{b})]\bar{T}[S_{\bar{n}}S_{n}^{\dagger}(0)]|0\rangle$
$\displaystyle=\tilde{S}(b)\Big{(}1+b^{2}\;C^{(s)}_{i}(b)\langle
0|O^{i}|0\rangle\Big{)}\;.$ (55)

To maintain the UV scale invariance of the cross section, we need that the
second term in the brackets be independent of UV divergences. However
additional rapidity divergences may be present in the non-perturbative matrix
element on the r.h.s. that cancel with the corresponding rapidity divergence
arising in the non-perturbative power corrections to the collinear or soft-
collinear functions. This is related to the origin of the non-perturbative
correction to the rapidity anomalous dimension and it is usually included also
in TMD analysis.

We can perform a similar analysis for the soft-collinear ($sc$) function. When
an $sc$ (perturbative) mode passes soft-drop, then it does not contribute to
$q_{T}$ since it becomes part of the groomed jet. But since it has a large +
component, it drives the groomed jet mass outside the region of measurement
and hence such events are dropped. Therefore, we only need to consider the
case when the $sc$ mode fails soft drop. In this case the non-perturbative
emission contributes to the $q_{T}$ measurement if it lies outside the groomed
jet. Given the angular scaling of this mode, which is much larger than the
collinear-soft ($cs$) and collinear modes that form the groomed jet, the phase
space region available is effectively unconstrained (this is also the reason
why we ignore any phase space constraints on the soft non-perturbative
emissions). Hence the correction in this case will also be a simple shift type
and is implemented in the same manner as in the case of the global soft
function. As before, we can pull out a common perturbative factor (that
includes the perturbative soft drop condition), and write

$\displaystyle\tilde{S}_{sc}^{\perp}(b,z_{\text{cut}}){|}_{\text{hadr.}}=\tilde{S}_{sc}^{\perp}(b,Qz_{\text{cut}})\Big{(}1+b^{2}C^{(sc)}_{i}(b,z_{\text{cut}})\;\langle
0|O^{i}|0\rangle\Big{)}.$ (56)

Notice that now all the $z_{\mathrm{cut}}$ dependence of the power correction
is included in the perturbative calculable coefficient
$C^{(sc)}(b,z_{\text{cut}})$, which multiplies the same non-perturbative power
correction present also in the global soft function case. The calculation of
$C^{(s)},\;C^{(sc)}$ is doable perturbatively, although this consideration
goes beyond the present work.

We can then combine all shift corrections that have an unconstrained phase
space for non-perturbative emissions together so that in $b$ space we have a
multiplicative correction to the perturbative cross section of the form

$\displaystyle
SS_{sc}^{\perp}{|}_{\text{hadr.}}=(1+b^{2}(\Omega_{s}+\Omega_{sc}))SS_{sc}^{\perp}{|}_{\text{pert.}}\;,$
(57)

where $\Omega_{s}$ is the same as the TMD case and $\Omega_{s}$ is a single
parameter to be fitted from $e^{+}e^{-}$ experiments. It is clear that, in the
event of non-trivial $C^{\\{(s),\,(sc)\\}}$, $\Omega_{s,\,sc}$ can have a mild
(logarithmic) dependence on $q_{T}$ so that this model will work well over a
limited range of $q_{T}$ which may be sufficient for most cases.

We now consider the shift corrections coming from the collinear-soft and the
collinear functions. Since these modes determine the region of the groomed
jet, we can consider two possible scenarios which give a non-trivial power
correction.

  1. 1.

Collinear-soft ($cs$) particles pass soft-drop:
If the $cs$ particles pass the soft-drop (for phase space see figure 9(a))
then any non-perturbative emission scaling as the $cs$ mode can contribute to
$q_{T}$ when it lies outside the groomed jet. In this case, we need to
calculate the catchment area of the groomed jet that is determined by the
angular distance of the $cs$ subject that passed soft-drop. As was pointed out
in Hoang:2019ceu , it is possible at NLL, using a coherent branching
formalism, to factorize a purely non-perturbative function from all the
calculable perturbative effects (including grooming). A detailed analysis of
these corrections will be presented in a future work.

Figure 9: (a) When the collinear-soft ($cs$) function passes soft drop, the
non-perturbative (NP) emissions, with the angular scaling of the $cs$ mode ,
with a virtuality $\Lambda_{QCD}$ must fall in the phase space shown by the
blue shaded area in order to contribute to $q_{T}$. (b) When the $cs$ function
fails soft drop, the NP emission with the angular scaling of the collinear
modes must not be clustered with the collinear sub-jet in order to contribute
to $q_{T}$.

  2. 2.

Collinear-soft particles fail soft-drop:
In this case collinear modes are the only ones that pass soft-drop (for phase
space see figure 9(b)), so that any non-perturbative mode scaling as $cs$ has
an unconstrained phase space, by the same logic as for the soft and the $sc$
functions, so that we get a simple shift correction of the same form as the
soft, $sc$ and TMD collinear functions.555Technically in this case the
perturbative value of $p_{\perp cs}$ would give a larger correction. However,
this correction can eventually be handled by transitioning to a new EFT in
which the $sc$ and $cs$ functions merge together. For now we will ignore them
and only keep track of the other non-perturbative corrections. There is
another possible interesting correction that will come from the collinear NP
emission that lies outside the catchment region that is now determined by the
collinear modes alone.

In this case there are two ways of approaching the problem. In one, we
consider separating out the non-perturbative corrections before factorizing
the $cs$ and collinear modes. The other way is to realize that in the case
where $cs$ fails soft-drop, the entire groomed jet mass measurement comes from
the jet function alone and using this condition we can define a catchment area
for the collinear non-perturbative emissions without explicitly accessing any
information from the $cs$ function, so that the factorization between the
collinear and $cs$ modes is maintained. In this case, we can do a diagrammatic
analysis, similar to Hoang:2019ceu , for the collinear function, to check if
it is possible to factorize the non-perturbative effects from the
perturbative. We leave this work for the future.

###  4.2 Boundary corrections

We now consider boundary corrections that leave the $q_{T}$ measurement
function unchanged but only require an expansion of the soft-drop condition in
$q^{-}/Q$. The functions that do not explicitly have a soft-drop condition can
then be ignored, which leaves us with only the $sc$ and $cs$ functions. We can
follow the same line of reasoning as in Hoang:2019ceu .

Figure 10: (a) The case where the sc subjet loses an NP emission (b) The case
when the sc subjet gains an NP emission

  1. 1.

$sc$ emissions
In this case we demand that either an addition or removal of the non-
perturbative emission cause the soft-collinear function to fail soft-drop.
Otherwise it will drive up the jet mass outside the measured range. If we
consider a non-perturbative emission $q^{\mu}$ along with a perturbative
momentum $p^{\mu}$, then we can expand out the soft-drop condition in the non-
perturbative momentum. We can write the complete measurement function as

$\displaystyle\Theta^{p\pm
q}=\Theta\left(\frac{p+q}{E_{J}}-z_{\text{cut}}\right)\delta^{2}(\boldsymbol{p}_{\perp
sc}-\boldsymbol{p}_{\perp}\mp\boldsymbol{q}_{T})\;,$ (58)

where $p$ is the momentum of the perturbative $sc$ sub-jet while $q^{\mu}$ is
the momentum of the non-perturbative emission. The $\pm$ signs indicate
whether the perturbative $cs$ subject gains or loses a non-perturbative
momentum after hadronization. In the case where the $sc$ sub-jet gains a non-
perturbative emission, the measurement expanded to leading order looks like

$\displaystyle\Theta^{p+q}\approx\Theta_{sd}^{p}\delta^{2}(\boldsymbol{p}_{\perp
sc}-\boldsymbol{p}_{\perp})+\frac{q^{-}}{E_{J}}\Theta^{\text{b.c.}}(\theta_{q},\theta_{p},\Delta\phi)\delta^{p}_{sd}\Big{[}\delta^{2}(\boldsymbol{p}_{\perp
sc}-\boldsymbol{p}_{\perp})\Big{]}\;,$ (59)

with

$\displaystyle\Theta^{p}_{sd}$
$\displaystyle\equiv\Theta\left(\frac{p}{E_{J}}-z_{\text{cut}}\right)\;,$
$\displaystyle\delta^{p}_{sd}$
$\displaystyle\equiv\delta\left(\frac{p}{E_{J}}-z_{\text{cut}}\right)\;.$ (60)

In this case, the non-perturbative emission $q^{\mu}$ gets clustered with the
$sc$ subject. Note that we have expanded $q_{i}$ from the
$\boldsymbol{p}_{\perp}$ measurement since we are working at leading order.
The phase-space constraint, $\Theta^{\text{b.c.}}$ (see figure 10(a)), , gives
the condition that ensures $q^{\mu}$ gets clustered with the $sc$ part.

The second case is when $q^{\mu}$ is emitted off $p^{\mu}$ but it is not
clustered with the $sc$ jet. The short distance condition now acts on $p-q$,
which can then be expanded out to give

$\displaystyle\Theta^{p-q}\approx\Theta_{sd}^{p}\delta^{2}(\boldsymbol{p}_{\perp
sc}-\boldsymbol{p}_{\perp})-\frac{q^{-}}{E_{J}}\bar{\Theta}^{\text{b.c.}}(\theta_{q},\theta_{p},\Delta\phi)\delta^{p}_{sd}\Big{[}\delta^{2}(\boldsymbol{p}_{\perp
sc}-\boldsymbol{p}_{\perp})\Big{]}\;,$ (61)

$\bar{\Theta}^{\text{b.c.}}$ (see figure 10(b)), is the phase space region for
$q^{\mu}$ so that it falls outside the sc subjet. We can see that the leading
power correction scales as $q^{-}/E_{J}$, which, given the angular scaling of
the $sc$ mode, scales as $q_{T}z_{cut}/Q$. Given a typical value of
$z_{cut}\sim 0.1$, this factor is then comparable to the $q_{T}^{2}/Q^{2}$
correction that we get from the shift terms.

  2. 2.

Soft -Collinear function
We expect that since perturbatively this function does not contribute to
$q_{T}$, the boundary correction should have no effect on the $q_{T}$
measurement.

We now have listed out all the possible NP corrections to the transverse
momentum measurement.

##  5 Conclusions

In this paper, we have presented the computation of the transverse momentum
de-correlation observable for fat jets groomed using the Soft-Drop algorithm.
We consider two scattering experiments: $e^{+}e^{-}\rightarrow$ di-jets and
semi-inclusive DIS. In the former, we measure the transverse momentum
imbalance between the two groomed jets. We impose a jet mass constraint on our
jets in order to ensure collimated jet configurations. Simulation using PYTHIA
show that grooming greatly reduces the impact of underlying events as well as
final state hadronization. We show that the factorization theorem for this
observable involves the universal soft function which also appears in the
traditional definition of TMDs. We propose that this observable can be used as
a probe of the non-perturbative rapidity anomalous dimension, which is a
universal parameter for TMD distributions. We prove within our EFT that the
cumulant jet mass constraint only adds to the overall normalization of the
perturbative cross section and hence does not impact the shape of the
transverse momentum distribution although it does contribute to the
uncertainty. We gather or compute all the ingredients necessary to evaluate
the cross section to NNLL accuracy and a numerical study for the cases of
interest. In the implementation we have used the artemide code web ;
Scimemi:2017etj ; Bertone:2019nxa  which contains the most recent extraction
TMDPDF at higher perturbative orders. As part of the numerical analysis we
have used the $\zeta$-prescription Scimemi:2018xaf  which allows a complete
disentanglement of non-perturbative effects of rapidity evolution from the
rest. An uncertainty analysis gives us an error band of approximately $\pm$ 10
$\%$. The main ingredient of this error is the perturbative uncertainty which
can be systematically improved. As shown in fig. 1 the hadronization
corrections at low $q_{T}$ are significantly smaller than the case of a
standard jet axis and it is therefore one of the major advantage of using
grooming. These effects are expected to be the same in $e^{+}e^{-}$ and SIDIS
because of the factorization of the cross section. In the case of $e^{+}e^{-}$
these corrections constitute all of non-perturbative effects and they are
associated with the final state shower. In order to do a meaningful extraction
of non-perturbative parameters in this case, it is therefore necessary to
improve the uncertainty from perturbative physics to be better than 5%. This
can be achieved by moving to a higher order in resummation accuracy (N3LL).
This is something we leave as a follow up to this paper.

In the SIDIS case we measure the transverse momentum imbalance between the
groomed jet and the recoiling lepton. Once again we demand a jet mass
measurement in order to ensure sensitivity to collinear physics only. A large
part of the contribution to this comes from the soft and collinear radiation
that lies outside the jet and, for low transverse momentum, probes the
complete TMDPDF. The cross section is again presented to NNLL accuracy and
involve much of the same ingredients as in the case of $e^{+}e^{-}\rightarrow$
dijets. A higher order perturbative calculation is expected to reduce
significatevely errors also in this case.

Concerning the hadronization effects we observe that grooming the jet allows
us to have a wide angle jet, which is preferred in low energy experiments,
while still being free from non-global logarithms, which are non-factorizable
and they are usually present in un-groomed jets. Nevertheless it is possible
to measure directly the hadronization effects due to grooming. The idea is to
parametrize and extract all of the non-perturbative effects from
$e^{+}e^{-}\rightarrow$ dijets and use them in SIDIS since they contain all
the same matrix elements (in addition to the TMDPDF) as explained in sec. 4.
This gives us a robust way to access the TMDPDF while maintaining control over
all other uninteresting non-perturbative effects.

## Acknowledgements

The authors would like to thank Aditya Pathak, Iain W. Stewart, and Wouter J.
Waalewijn for useful discussions. I.S. likes to acknowledge the support of Los
Alamos National Lab for his visit, during which part of the work was done.
D.G.R. and I.S. are supported by the Spanish MECD grant FPA2016-75654-C2-2-P.
This project has received funding from the European Union Horizon 2020
research and innovation program under grant agreement No 824093 (STRONG-2020).
D.G.R. acknowledges the support of the Universidad Complutense de Madrid
through the predoctoral grant CT17/17-CT18/17. Y.M. and V.V. are supported by
the U.S. Department of Energy through the Office of Science, Office of Nuclear
Physics under Contract DE-AC52-06NA25396 and by an Early Career Research
Award, through the LANL/LDRD Program. V.V. is also supported within the
framework of the TMD Topical Collaboration. L.Z. is supported by ERC grant
ERC- STG-2015-677323.

##  Appendix A Laplace and Fourier transformations

We define the Fourier transform,
$\mathcal{FT}[f](\boldsymbol{b})={f}(\boldsymbol{b})$ of a function,
$f(\boldsymbol{q}_{T})=\mathcal{FT}^{-1}[{f}](\boldsymbol{q}_{T})$ as follows,

${f}(\boldsymbol{b})=\int_{-\infty}^{+\infty}d\boldsymbol{q}_{T}\;f(\boldsymbol{q}_{T})\exp(-i\boldsymbol{b}\cdot\boldsymbol{q}_{T})\;,$
(62)

and the inverse transform

$f(\boldsymbol{q}_{T})=\int_{-\infty}^{+\infty}\frac{d\boldsymbol{b}\;}{(2\pi)^{2}}{f}(\boldsymbol{b})\exp(i\boldsymbol{b}\cdot\boldsymbol{q}_{T})\;.$
(63)

In order to get the Fourier transforms of the plus distributions that appear
in the factorization theorem we use,

$\frac{1}{(2\pi)\mu^{2}}\Big{(}\frac{\mu^{2}}{q_{T}^{2}}\Big{)}^{1+\alpha}=-\frac{1}{2\alpha}\delta^{(2)}(\boldsymbol{q}_{T})+\mathcal{L}_{0}(q_{T}^{2},\mu^{2})-\alpha\mathcal{L}_{1}(q_{T}^{2},\mu^{2})+\mathcal{O}(\alpha^{2})\;.$
(64)

Taking the Fourier transform of the left-hand-side (LHS) we get (see eq. (E.2)
of ref. Chiu:2012ir )

$\displaystyle\int_{-\infty}^{+\infty}\frac{d\boldsymbol{q}_{T}}{(2\pi)}\;\frac{1}{\mu^{2}}\Big{(}\frac{\mu^{2}}{q_{T}^{2}}\Big{)}^{1+\alpha}\exp(-i\boldsymbol{b}\cdot\boldsymbol{q}_{T})$
$\displaystyle=-\frac{e^{-2\alpha\gamma_{E}}}{2\alpha}\frac{\Gamma(1-\alpha)}{\Gamma(1+\alpha)}\Big{(}\frac{\mu}{\mu_{E}}\Big{)}^{2\alpha}$
$\displaystyle=-\frac{1}{2\alpha}+\ln\Big{(}\frac{\mu}{\mu_{E}}\Big{)}+\alpha\ln^{2}\Big{(}\frac{\mu}{\mu_{E}}\Big{)}+\mathcal{O}(\alpha^{2})\;,$
(65)

where $\mu_{E}=2\exp(-\gamma_{E})/b$ and $b\equiv|\boldsymbol{b}|$ and in the
second line we expanded in $\alpha$. Comparing this result with the RHS of eq.
(64) we get,

$\displaystyle\mathcal{FT}\Big{[}\delta^{(2)}(\boldsymbol{q}_{T})\Big{]}(\boldsymbol{b})$
$\displaystyle=1\;,$
$\displaystyle\mathcal{FT}\Big{[}\mathcal{L}_{0}(q_{T}^{2},\mu^{2})\Big{]}(\boldsymbol{b})$
$\displaystyle=\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}\;,$
$\displaystyle\mathcal{FT}\Big{[}\mathcal{L}_{1}(q_{T}^{2},\mu^{2})\Big{]}(\boldsymbol{b})$
$\displaystyle=\ln^{2}\Big{(}\frac{\mu_{E}}{\mu}\Big{)}\;.$ (66)

We define the convolution $f\otimes g$ with

$\Big{[}f\otimes g\Big{]}(\boldsymbol{q}_{T})=\int
d\boldsymbol{\ell}_{\perp}\;f(\boldsymbol{q}_{T}-\boldsymbol{\ell}_{\perp})g(\boldsymbol{\ell}_{\perp})\;,$
(67)

such that

$\mathcal{F}\Big{[}f\otimes
g\Big{]}(\boldsymbol{b})={f}(\boldsymbol{b})\times{g}(\boldsymbol{b})\;.$ (68)

Similarly for the distribution in the jet-thrust we often work in Laplace
space where the corresponding convolutions translate to products. For these
reason we define the Laplace transformation $\mathcal{LT}[f](s)=f(s)$ of jet-
trust distribution $f(e)=\mathcal{LT}^{-1}[f](e)$ as follows:

$f(s)=\int_{-\infty}^{\infty}de\;\exp(-se)f(e)\;,$ (69)

and the corresponding inverse transform

${f}(e)=\frac{1}{2\pi
i}\int_{\gamma-i\infty}^{\gamma+i\infty}ds\;\exp(se)f(s)\;.$ (70)

Similarly with the case of Fourier transform we use the following expansion to
identify the Laplace transform of plus distributions that are present in the
fixed order expansion of the jet and collinear-soft functions,

$\frac{1}{\xi}\Big{(}\frac{\xi}{e}\Big{)}^{1+\alpha}{|}_{e>0}=-\frac{1}{\alpha}\delta(e)+\mathcal{L}_{0}(e,\xi)-\alpha\mathcal{L}_{1}(e,\xi)+\mathcal{O}(\alpha^{2})\;,$
(71)

taking the Laplace transform of the LHS we get

$\int_{0}^{\infty}\frac{de}{\xi}\Big{(}\frac{\xi}{e}\Big{)}^{1+\alpha}\exp(-se)=s^{\alpha}\Gamma(-\alpha)=-\frac{1}{\alpha}-\ln(\xi\tilde{s})-\alpha\Big{(}\frac{1}{2}\ln^{2}(\xi\tilde{s})+\frac{\pi^{2}}{12}\Big{)}+\mathcal{O}(\alpha^{2})\;,$
(72)

where $\tilde{s}\equiv s\exp(\gamma_{E})$ and thus from comparing eq. (71) and
(72) we have

$\displaystyle\mathcal{LT}\Big{[}\delta(e)\Big{]}(s)$ $\displaystyle=1\;,$
$\displaystyle\mathcal{LT}\Big{[}\mathcal{L}_{0}(e,\xi)\Big{]}(s)$
$\displaystyle=-\ln(\xi\tilde{s})\;,$
$\displaystyle\mathcal{LT}\Big{[}\mathcal{L}_{1}(e,\xi)\Big{]}(s)$
$\displaystyle=\frac{1}{2}\ln^{2}(\xi\tilde{s})+\frac{\pi^{2}}{12}\;.$ (73)

##  Appendix B Operator definitions and one loop results

In this appendix we give the operator definitions of the factorization
elements and their NLO expansions. From those we determine the renormalization
functions, group equations, and corresponding anomalous dimensions. Many of
the results presented here are already known and found in literature.

###  B.1 Jet functions

The quark and gluon jet function definitions, one loop calculation, and the
corresponding Laplace transforms can be found in ref. Frye:2016aiz . Here we
summarize their results. The quark jet function is given by,

$J_{q}(e,Q)=\frac{(2\pi)^{3}}{N_{c}}\text{tr}\;{\langle}\frac{\not{\bar{n}}}{2}\chi_{n}(0)\delta(Q-\mathcal{P}^{-})\delta^{(2)}(\boldsymbol{\mathcal{P}}_{\perp})\delta(e-\mathcal{E})\bar{\chi}_{n}{\rangle}\;,$
(74)

and the gluon

$J_{q}(e,Q)=\frac{(2\pi)^{3}}{N_{c}}\text{tr}\;{\langle}\frac{\not{\bar{n}}}{2}\mathcal{B}^{\mu}_{n\perp}(0)\delta(Q-\mathcal{P}^{-})\delta^{(2)}(\boldsymbol{\mathcal{P}}_{\perp})\delta(e-\mathcal{E})\mathcal{B}_{n\perp\mu}{\rangle}\;,$
(75)

where $N_{c}$ is the number of colors and $\mathcal{B}_{n\perp}^{\mu}$ is the
gauge invariant gluon building block of the effective field theory,

$\mathcal{B}_{n\perp}^{\mu}=\frac{1}{g}[W_{n}^{{\dagger}}(\mathcal{P}_{\perp}^{\mu}+gA_{n,\perp}^{\mu})W_{n}]\;.$
(76)

As demonstrated earlier when working with the cumulant distribution (i.e.,
when integrating out to $e_{\text{cut}}$) it is useful to work in Laplace
space. The renormalized groomed jet function up to NLO contributions in
Laplace space is given by

$J_{i}(s,Q;\mu)=1+\frac{\alpha_{s}C_{i}}{2\pi}\Big{\\{}L_{J}^{2}+\bar{\gamma}_{i}L_{J}-\frac{\pi^{2}}{3}+c_{i}\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(77)

where for quark initiated jets we have

$\displaystyle C_{q}$ $\displaystyle=C_{F}=\frac{N_{c}^{2}-1}{2N_{c}}\;,$
$\displaystyle\bar{\gamma}_{q}$ $\displaystyle=\frac{3}{2}\;,$ $\displaystyle
c_{q}$ $\displaystyle=\frac{7}{2}\;,$ (78)

and for gluon initiated jets we have

$\displaystyle C_{g}$ $\displaystyle=C_{A}=N_{c}\;,$
$\displaystyle\bar{\gamma}_{g}$ $\displaystyle=\frac{\beta_{0}}{2C_{A}}\;,$
$\displaystyle c_{g}$
$\displaystyle=\frac{67}{18}-\frac{10}{9}\frac{n_{f}T_{R}}{C_{A}}\;,$ (79)

The logarithms, $L_{J}$ that appear in eq. (77) and the corresponding one loop
anomalous dimensions are

$\displaystyle L_{J}$
$\displaystyle=\ln\Big{(}\frac{\mu^{2}\tilde{s}}{Q^{2}}\Big{)}\;,$
$\displaystyle\gamma^{J}=\frac{\alpha_{s}C_{i}}{\pi}\Big{(}2L_{J}+\bar{\gamma}_{i}\Big{)}+\mathcal{O}(\alpha_{s}^{2})\;.$
(80)

The anomalous dimension is defined through the RG equation satisfied by
renormalized jet functions. In Laplace space this is

$\frac{d}{d\ln\mu}J_{i}(s,Q;\mu)=\gamma^{J}(s,Q;\mu)J_{i}(s,Q;\mu)\;.$ (81)

In momentum space the above equation is written as convolution (in the
invariant mass variable $e$), of the anomalous dimension and the renormalized
jet function.

###  B.2 Collinear-soft function

The operator definition of the invariant mass measurement collinear soft
function is given by

$S_{cs}(e,Qz_{cut})=\frac{1}{N_{R}}\text{tr}{\langle}T\left(U_{n}^{\dagger}W_{t}\right)\mathcal{M}^{SD}_{e}\bar{T}\left(W_{t}^{\dagger}U_{n}\right){\rangle}\;,$
(82)

where $\mathcal{M}^{SD}_{e}$ is the invariant measurement function,

$\mathcal{M}^{SD}_{e}=\delta\left(e-(1-\Theta_{SD})\;\mathcal{E}\right)\;.$
(83)

Here we dropped the jet flavor (quark/anti-quark or gluon) for simplicity of
notation and the normalization constant $N_{R}$ is simply the size of the
representation for SU($N_{c}$) of the $W_{t}$ and $U_{n}$ Wilson lines. For
quark jets (fundamental representation) we have $N_{R}=N_{c}$ and for gluon
jets (adjoint representation) we have $N_{R}=N_{c}^{2}-1$. At NLO the bare
collinear soft function is given by

$S_{cs,\text{bare}}(e,Qz_{\text{cut}})=\delta(e)+\frac{\alpha_{s}C_{i}}{\pi}\Big{\\{}-\frac{1}{\epsilon^{2}}\delta(e)+\frac{1}{\epsilon}\mathcal{L}_{0}(e,\xi)-\mathcal{L}_{1}(e,\xi)+\frac{\pi^{2}}{12}\delta(e)\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(84)

where

$\xi\equiv\frac{\mu^{2}}{Q^{2}z_{\text{cut}}}\;.$ (85)

Therefore we have for the renormalized function

$S_{cs}(e,Qz_{\text{cut}})=\delta(e)+\frac{\alpha_{s}C_{i}}{\pi}\Big{\\{}-\mathcal{L}_{1}(e,\xi)+\frac{\pi^{2}}{12}\delta(e)\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(86)

where

$S_{cs,\text{bare}}(e,Qz_{\text{cut}})=Z_{cs}\otimes
S_{cs}(e,Qz_{\text{cut}})\;,$ (87)

with

$Z_{cs}(e)=\delta(e)+\frac{\alpha_{s}C_{i}}{\pi}\Big{\\{}-\frac{1}{\epsilon^{2}}\delta(e)+\frac{1}{\epsilon}\mathcal{L}_{0}(e,\xi)\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;.$
(88)

In Laplace space for the renormalized collinear-soft function we get,

$S_{cs}(s,Qz_{\text{cut}};\mu)=1-\frac{\alpha_{s}C_{i}}{2\pi}L_{cs}^{2}+\mathcal{O}(\alpha_{s}^{2})\;,$
(89)

which satisfies the following RGE

$\frac{d}{d\ln\mu}S_{cs}(s,Qz_{\text{cut}};\mu)=\gamma^{cs}(s,\mu)S_{cs}(s,Qz_{\text{cut}};\mu)\;.$
(90)

The logarithm $L_{cs}$ and the corresponding anomalous dimension are

$\displaystyle L_{cs}$ $\displaystyle=\ln(\xi\tilde{s})\;,$
$\displaystyle\gamma^{cs}(s,\mu)$
$\displaystyle=-2\frac{\alpha_{s}C_{i}}{\pi}L_{cs}+\mathcal{O}(\alpha_{s}^{2})\;.$
(91)

###  B.3 Soft function

The soft function that appears in the factorization theorems in eq. (15) and
(32) is defined in eq. (7) and it has been calculated in several schemes at
higher orders in QCD, as reported in sec. 2.2. Here we report a one loop
expression using the analytic regulator in momentum space,

$S_{\text{bare}}=\delta^{(2)}(\boldsymbol{q}_{T})+\frac{\alpha_{s}(\mu)C_{i}}{\pi}\Big{\\{}\frac{4}{\eta}\Big{[}\mathcal{L}_{0}(q_{T}^{2},\mu^{2})-\frac{1}{2\epsilon}\delta^{(2)}(\boldsymbol{q}_{T})\Big{]}+\frac{1}{\epsilon}\Big{[}\frac{1}{\epsilon}-2\ln\Big{(}\frac{\nu}{\mu}\Big{)}\Big{]}\delta^{(2)}(\boldsymbol{q}_{T})\\\
+4\mathcal{L}_{0}(q_{T}^{2},\mu^{2})\ln\Big{(}\frac{\nu}{\mu}\Big{)}-2\mathcal{L}_{1}(q_{T}^{2},\mu^{2})-\frac{\pi^{2}}{12}\delta^{(2)}(\boldsymbol{q}_{T})\Big{\\}}+\mathcal{O}(\alpha_{s}^{2}).$
(92)

The renormalized soft function, $S$, is defined through

$S_{\text{bare}}=Z_{s}^{\perp}(\mu,\nu)\otimes S(\mu,\nu),$ (93)

and satisfies the following renormalization group equations

$\displaystyle\frac{d}{d\ln{\mu}}S(\mu,\nu)$
$\displaystyle=\gamma^{s}(\mu,\nu)S(\mu,\nu)\;,$
$\displaystyle\frac{d}{d\ln{\nu}}S(\mu,\nu)$
$\displaystyle=\gamma^{s}_{\nu}(\mu,\nu)\otimes S(\mu,\nu)\;.$ (94)

Therefore we find for the one-loop corresponding impact parameter space
quantities

$S(\mu,\nu)=1+\frac{\alpha_{s}(\mu)C_{i}}{\pi}\Big{\\{}4\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}\ln\Big{(}\frac{\nu}{\mu}\Big{)}-2\ln^{2}\Big{(}\frac{\mu_{E}}{\mu}\Big{)}-\frac{\pi^{2}}{12}\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(95)
$Z_{s}^{\perp}(\mu,\nu)=1+\frac{\alpha_{s}(\mu)C_{i}}{\pi}\Big{\\{}\frac{4}{\eta}\Big{[}\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}-\frac{1}{2\epsilon}\Big{]}+\frac{1}{\epsilon}\Big{[}\frac{1}{\epsilon}-2\ln\Big{(}\frac{\nu}{\mu}\Big{)}\Big{]}\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(96)

with

$\displaystyle\gamma^{s}(\mu,\nu)$
$\displaystyle=-4\frac{\alpha_{s}(\mu)C_{i}}{\pi}\ln\Big{(}\frac{\nu}{\mu}\Big{)}+\mathcal{O}(\alpha_{s}^{2})\;,$
$\displaystyle\gamma^{s}_{\nu}(\mu,\nu)$
$\displaystyle=4\frac{\alpha_{s}(\mu)C_{i}}{\pi}\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}+\mathcal{O}(\alpha_{s}^{2})\;.$
(97)

The rapidity and renormalization scales used to produce our result are fixed
using the $\zeta$-prescription Scimemi:2018xaf  adapted for this case. Later
in the appendix we give a description of how one can use the rapidity
regulated objects that have $\nu$ dependence to construct the subtracted
rapidity divergences free objects but yet keep trace of the rapidity logs
using the $\zeta$ parameter.

###  B.4 Soft-collinear function

Figure 11: The order $\mathcal{O}(\alpha_{s})$ diagrams that contribute to the
soft-collinear function.

The soft-collinear function is defined by the matrix element

$S_{sc}^{\perp}(Qz_{cut})=\frac{1}{N_{R}}\text{tr}{\langle}T\left(U_{n}^{\dagger}W_{t}\right)\mathcal{M}^{SD}_{\perp}\bar{T}\left(W_{t}^{\dagger}U_{n}\right){\rangle}\;,$
(98)

and the groomed jet measurement function, $\mathcal{M}^{SD}_{\perp}$ is given
in terms of the label momentum operator, $\mathcal{P}$,

$\mathcal{M}^{SD}_{\perp}=\Theta_{SD}\times\delta^{2}\left(\boldsymbol{q}_{T}-\Theta_{SD}\boldsymbol{\mathcal{P}}_{\perp}\right)\;,$
(99)

where $\Theta_{SD}$ denotes the soft drop groomer. The collinear-soft modes
only contribute to the invariant mass measurement if they pass soft-drop,
which is implemented by the $\Theta_{SD}$ term. The NLO calculation involves
one real and one virtual diagram shown in fig. 11. While the virtual diagram
is scaleless. The diagram with a real gluon needs to be integrated over the
phase-space of soft gluon. This then yields non-vanishing contribution from
when the soft gluon fails the grooming,

$S_{sc,\text{NLO}}^{\perp}(Qz_{cut})=4g^{2}C_{i}\tilde{\mu}^{2\epsilon}\;\nu^{\eta}\;\int\frac{d^{d}k}{(2\pi)^{d-1}}\frac{\delta(k^{2})\;\delta^{(2)}(\boldsymbol{q}_{T}-\boldsymbol{k}_{\perp})}{k^{+}\;(k^{-})^{1+\eta}}\theta(Qz_{\text{cut}}-k^{-})\;.$
(100)

Performing the integrals we find for the bare quantity

$S_{sc,\text{bare}}^{\perp}(Qz_{\text{cut}})=\delta^{(2)}(\boldsymbol{q}_{T})+\frac{\alpha_{s}C_{i}}{\pi}\Big{\\{}-\frac{2}{\eta}\Big{[}\mathcal{L}_{0}(q_{T}^{2},\mu^{2})-\frac{1}{2\epsilon}\delta^{(2)}(\boldsymbol{q}_{T})\Big{]}+\frac{1}{\epsilon}\ln\Big{(}\frac{\nu}{Qz_{\text{cut}}}\Big{)}\delta^{(2)}(\boldsymbol{q}_{T})\\\
-2\ln\Big{(}\frac{\nu}{Qz_{\text{cut}}}\Big{)}\mathcal{L}_{0}(q_{T}^{2},\mu^{2})\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(101)

and for the renormalized quantity, $S_{sc,}^{\perp}(Qz_{\text{cut}};\mu,\nu)$
we have

$S_{sc,\text{bare}}^{\perp}(Qz_{\text{cut}})=Z_{sc}^{\perp}(\mu,\nu)\otimes
S_{sc}^{\perp}(Qz_{\text{cut}};\mu,\nu)\;,$ (102)

and satisfies the following renormalization group equations

$\displaystyle\frac{d}{d\ln{\nu}}S_{sc}^{\perp}(\mu,\nu)$
$\displaystyle=\gamma^{sc}_{\nu}(\mu,\nu)\otimes S_{sc}^{\perp}(\mu,\nu)\;,$
$\displaystyle\frac{d}{d\ln{\mu}}S_{sc}^{\perp}(\mu,\nu)$
$\displaystyle=\gamma^{sc}(\mu,\nu)S_{sc}^{\perp}(\mu,\nu)\;,$ (103)

where the $Qz_{\text{cut}}$ dependence is suppressed to improve readability.
In $\overline{\text{MS}}$ scheme the corresponding Fourier transform can be
obtained using eq. (A):

$\tilde{S}_{sc}^{\perp}(Qz_{\text{cut}};\mu,\nu)=1+\frac{\alpha_{s}C_{i}}{\pi}\Big{\\{}-2\ln\Big{(}\frac{\nu}{Qz_{\text{cut}}}\Big{)}\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(104)
$\tilde{Z}_{sc}^{\perp}(\mu,\nu)=1+\frac{\alpha_{s}C_{i}}{\pi}\Big{\\{}-\frac{2}{\eta}\Big{[}\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}-\frac{1}{2\epsilon}\Big{]}+\frac{1}{\epsilon}\ln\Big{(}\frac{\nu}{Qz_{\text{cut}}}\Big{)}\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;,$
(105)

and thus for the one-one-loop anomalous dimensions we get

$\displaystyle\gamma_{\nu}^{sc}(\mu,\nu)$
$\displaystyle=-2\frac{\alpha_{s}(\mu)C_{i}}{\pi}\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}+\mathcal{O}(\alpha_{s}^{2})\;$
$\displaystyle\gamma^{sc}(\mu,\nu)=2\frac{\alpha_{s}(\mu)C_{F}}{\pi}\ln\Big{(}\frac{\nu}{Qz_{\text{cut}}}\Big{)}+\mathcal{O}(\alpha_{s}^{2}).$
(106)

##  Appendix C Solution of renormalization group evolution equations

In this appendix we discuss the solutions of both virtuality and rapidity
renormalization group equations written in eq. (23). All factorization
elements (hard, soft, soft-collinear, collinear-soft, and jet) satisfy
renormalization group equations, but only transverse momentum dependent
quantities have rapidity RGE.

###  C.1 Renormalization group evolution

The solution to the RGE in eq. (23) is

$\displaystyle G(\mu)$
$\displaystyle=\mathcal{U}_{G}(\mu,\mu_{0})G(\mu_{0})\,,$
$\displaystyle\mathcal{U}_{G}(\mu,\mu_{0})=\exp\left(K_{G}(\mu,\mu_{0})\right)\left(\frac{\mu_{0}}{m_{G}}\right)^{2\;\omega_{G}(\mu,\mu_{0})},$
(107)

with

$\displaystyle K_{G}(\mu,\mu_{0})$
$\displaystyle=2\int_{\alpha(\mu_{0})}^{\alpha(\mu)}\frac{d\alpha}{\beta[\alpha]}\Gamma^{G}[\alpha]\int_{\alpha(\mu_{0})}^{\alpha}\frac{d\alpha^{\prime}}{\beta[\alpha^{\prime}]}+\int_{\alpha(\mu_{0})}^{\alpha(\mu)}\frac{d\alpha}{\beta[\alpha]}\Delta\gamma^{G}[\alpha],$
(108) $\displaystyle\omega_{G}(\mu,\mu_{0})$
$\displaystyle=\int_{\alpha(\mu_{0})}^{\alpha(\mu)}\frac{d\alpha}{\beta[\alpha]}\Gamma^{G}[\alpha].$
(109)

Since in this work we are interested only in the NLL and NLL’ result we may
keep only the first two terms in the perturbative expansion of the cusp part
(i.e., $\Gamma_{0}^{G}$, $\Gamma_{0}^{\text{cusp}}$, and
$\Gamma_{1}^{\text{cusp}}$) and only the first term form the non-cusp part
($\gamma^{G}_{0}$). Performing this expansion we get,

$\displaystyle K_{G}(\mu,\mu_{0})$
$\displaystyle=-\frac{\gamma^{G}_{0}}{2\beta_{0}}\ln
r-\frac{2\pi\Gamma^{G}_{0}}{(\beta_{0})^{2}}\Big{[}\frac{r-1-r\ln
r}{\alpha_{s}(\mu)}+\left(\frac{\Gamma_{1}^{\text{cusp}}}{\Gamma_{0}^{\text{cusp}}}-\frac{\beta_{1}}{\beta_{0}}\right)\frac{1-r+\ln
r}{4\pi}+\frac{\beta_{1}}{8\pi\beta_{0}}\ln^{2}r\Big{]},$ (110)
$\displaystyle\omega_{G}(\mu,\mu_{0})$
$\displaystyle=-\frac{\Gamma^{G}_{0}}{2\beta_{0}}\Big{[}\ln
r+\left(\frac{\Gamma^{1}_{\text{cusp}}}{\Gamma^{0}_{\text{cusp}}}-\frac{\beta_{1}}{\beta_{0}}\right)\frac{\alpha_{s}(\mu_{0})}{4\pi}(r-1)\Big{]},$
(111)

where $r=\alpha(\mu)/\alpha(\mu_{0})$ and $\beta_{n}$ are the coefficients of
the QCD $\beta$-function,

$\beta[\alpha_{s}]=\mu\frac{d\alpha_{s}}{d\mu}=-2\alpha_{s}\sum_{n=0}^{\infty}\left(\frac{\alpha_{s}}{4\pi}\right)^{1+n}\beta_{n}\;.$
(112)

The expressions for all ingredients necessary to perform the evolution of any
function that appears in the factorization theorems we considered in this
paper are given in tab. 1.

Function | $\Gamma^{G}_{0}$ | $\gamma^{G}_{0}$ | $m_{G}$
---|---|---|---
$H_{ij}$ | $-4(C_{i}+C_{\bar{j}})$ | $-4\bar{\gamma}_{i}(C_{i}+C_{j})$ | $Q$
$S_{cs}$ | $-8C_{i}$ | 0 | $Q\sqrt{z_{\text{cut}}/\tilde{s}}$
$J_{i}$ | $8C_{i}$ | $4\bar{\gamma}_{i}C_{i}$ | $Q/\sqrt{\tilde{s}}$
$B_{i/h}$ | 0 | $4\bar{\gamma}_{i}C_{i}+\gamma_{sc}^{0}$ | 0
$S$ | $4(C_{i}+C_{j})$ | 0 | $\nu_{s}$
$S_{sc}^{\perp}$ | 0 | $\gamma^{sc}_{0}$ | n.a
Table 1: Anomalous dimensions coefficients for up to NLL accuracy:
$\bar{\gamma}_{q}=3/2$, $\bar{\gamma}_{g}=\beta_{0}/(2C_{A})$, and
$\gamma^{sc}_{0}=2\alpha_{s}(\mu)C_{F}/\pi\ln(\nu/Qz_{\text{cut}})$.

The coefficients for the expansion of the cusp anomalous dimension are

$\displaystyle\Gamma^{\text{cusp}}_{0}=$ $\displaystyle 4C_{F}\;,$
$\displaystyle\Gamma^{\text{cusp}}_{1}=$ $\displaystyle
4C_{F}\Big{[}\Big{(}\frac{67}{9}-\frac{\pi^{2}}{3}\Big{)}C_{A}-\frac{20}{9}n_{f}T_{R}\Big{]}\;,$
$\displaystyle\Gamma^{\text{cusp}}_{2}=$ $\displaystyle
4C_{F}\Big{[}\Big{(}\frac{245}{6}-\frac{134}{27}\pi^{2}+\frac{11}{45}\pi^{4}+\frac{22}{3}\zeta_{3}\Big{)}C_{A}^{2}+\Big{(}-\frac{209}{108}+\frac{5}{27}\pi^{2}-\frac{7}{3}\zeta_{3}\Big{)}8C_{A}n_{f}T_{R}$
$\displaystyle+\Big{(}16\zeta_{3}-\frac{15}{3}C_{F}n_{f}T_{R}-\frac{64}{27}T_{R}^{2}n_{f}^{2}\Big{)}\Big{]}\;.$
(113)

The two loop non-cusp anomalous dimensions we need to NNLL RGEs are given by
ref. Frye:2016aiz

$\displaystyle\frac{1}{2}\gamma_{1}^{s}+\gamma_{1}^{sc}$
$\displaystyle=\frac{C_{i}}{2}\Big{[}34.01C_{F}+\Big{(}\frac{1616}{27}-56\zeta_{3}-9.31\Big{)}C_{A}-\Big{(}\frac{448}{27}+14.04\Big{)}n_{f}T_{R}-\frac{2}{3}\pi^{2}\beta_{0}\Big{]}\;,$
$\displaystyle\gamma^{cs}_{1}=$ $\displaystyle
C_{i}\Big{[}-17.00C_{F}+\Big{(}-55.20+\frac{22}{9}\pi^{2}+56\zeta_{3}\Big{)}C_{A}+\Big{(}23.61-\frac{8}{9}\pi^{2}\Big{)}n_{f}T_{R}\Big{]}$
$\displaystyle\gamma^{q}_{1}=$ $\displaystyle
C_{F}\Big{[}\Big{(}3-4\pi^{2}+48\zeta_{3}\Big{)}C_{F}+\Big{(}\frac{1769}{27}+\frac{22}{9}\pi^{2}-80\zeta_{3}\Big{)}C_{A}$
$\displaystyle+\Big{(}-\frac{484}{27}-\frac{8}{9}\pi^{2}\Big{)}n_{f}T_{R}\Big{]}\;,$
$\displaystyle\gamma_{1}^{s}+\gamma_{1}^{sc}+\gamma_{1}^{B}=$ $\displaystyle
C_{i}\Big{[}\Big{(}20-4\pi^{2}+48\zeta_{3}\Big{)}C_{F}+\Big{(}60.87+\frac{22}{9}\pi^{2}-80\zeta_{3}\Big{)}C_{A}$
$\displaystyle+\Big{(}-24.94-\frac{8}{9}\pi^{2}\Big{)}n_{f}T_{R}\Big{]}\;.$
(114)

###  C.2 The connection between $\zeta$-parameter and rapidity regulator

In the standard EFT approach one used the rapidity renormalization group (RRG)
equations in order to resum large logarithms at the level of individual
rapidity regulated terms Chiu:2011qc ; Chiu:2012ir . A more recent approach
for performing the resummation of large logarithms in the TMD evolution it was
introduced in ref. Scimemi:2018xaf . The approach is referred to as the
$\zeta$-prescription. Here we rewrite the fixed order results using the
rapidity regulator in the past sections in the form appropriate for
implementing the $\zeta$-prescription. In the framework of ref.
Scimemi:2018xaf  one works with the rapidity divergent free quantity,

$S^{sub}_{sc}(b;\mu,\zeta)\equiv\sqrt{S^{\perp}_{2}({\boldsymbol{b}};\mu,\nu_{s})}\;S_{sc}^{\perp}({\boldsymbol{b}},Qz_{\text{cut}};\mu,\nu_{sc})\;,$
(115)

where we have explicitly show the dependence on the rapidity regulator
parameters $\nu_{s}$ and $\nu_{sc}$. In the RRG approach this combination does
not acquire rapidity evolution thus here in order to establish the rapidity
evolution we fix the rapidity scales at two different values. Particularly we
evaluate the soft-collinear rapidity scale at its canonical value,
$\nu_{sc}=Qz_{\text{cut}}$, and we allow for the corresponding soft scale to
float through a parameter $\zeta$: $\nu_{s}=\sqrt{\zeta}$.666Note that this is
not a unique choice of scales since any choice for which
$\nu_{s}/\nu_{sc}=\sqrt{\zeta}/(Qz_{\text{cut}})$ will give the same result.
With this choice of scales we have,

$S^{sub}_{sc}(b;\mu,\zeta)=1+\frac{\alpha_{s}(\mu)C_{i}}{2\pi}\Big{\\{}2\ln\Big{(}\frac{\mu_{E}}{\mu}\Big{)}\ln\Big{(}\frac{\zeta}{\mu^{2}}\Big{)}-2\ln^{2}\Big{(}\frac{\mu_{E}}{\mu}\Big{)}-\frac{\pi^{2}}{12}\Big{\\}}+\mathcal{O}(\alpha_{s}^{2})\;.$
(116)

And according to the notation of eq. (20) and (21) satisfies the following
equations

$\displaystyle\mu^{2}\frac{d}{d\mu^{2}}S^{sub}_{sc}(b;\mu,\zeta)$
$\displaystyle=\frac{1}{2}\gamma^{sc/sub}(\mu,\zeta)S^{sub}_{sc}(b;\mu,\zeta)\;,$
$\displaystyle\zeta\frac{d}{d\zeta}S^{sub}_{sc}(b;\mu,\zeta)$
$\displaystyle=-\mathcal{D}(\mu)S^{sub}_{sc}(b;\mu,\zeta)\;.$ (117)

Its easy to show that the anomalous dimensions $\gamma^{sc/sub}$ and
$\mathcal{D}$ are related to the RG and RRG anomalous dimensions of the global
soft and soft-collinear function as follows,

$\gamma^{sc/sub}(\mu,\zeta)=\frac{1}{2}\gamma^{s}+\gamma^{sc}=\Gamma^{\text{cusp}}[\alpha_{s}]\ln\Big{(}\frac{\mu^{2}}{\zeta}\Big{)}+\frac{1}{2}\Delta\gamma^{s}[\alpha_{s}]+\Delta\gamma^{sc}[\alpha_{s}]\;,$
(118)

and

$\mathcal{D}(\mu)=\Gamma^{\text{cusp}}[\alpha_{s}]\ln\Big{(}\frac{\mu}{\mu_{E}}\Big{)}-\frac{1}{4}\Delta\gamma_{\nu}^{s}[\alpha_{s}]\;,$
(119)

where

$\displaystyle\Delta\gamma_{\nu}^{s}=-\Big{(}\frac{\alpha_{s}(\mu)}{4\pi}\Big{)}^{2}C_{i}\Big{[}\Big{(}\frac{128}{9}-56\zeta_{3}\Big{)}C_{A}+\frac{112}{9}\beta_{0}\Big{]}+\mathcal{O}(\alpha_{s}^{3}).$
(120)

It is easy to confirm by looking the above equations that the anomalous
dimensions $\gamma^{sc/sub}$ and $\mathcal{D}$ satisfy the following
differential equations,

$\displaystyle\frac{d}{d\ln\zeta}\gamma^{sc/sub}(\mu,\zeta)$
$\displaystyle=-\Gamma_{\text{cusp}}\;,$
$\displaystyle\frac{d}{d\ln\mu}\mathcal{D}(\mu)$
$\displaystyle=+\Gamma_{\text{cusp}}\;.$ (121)

Also comparing against the notation of eq. (23) we see that the non-cusp part,
$\Delta\gamma^{sc/sub}$, of the anomalous dimension $\gamma^{sc/sub}$ is a
linear combination of the corresponding non-cusp pieces of the global soft and
soft-collinear functions. Particularly:

$\Delta\gamma^{sc/sub}(\mu)=\Big{(}\frac{1}{2}\Delta\gamma^{s}[\alpha_{s}(\mu)]+\Delta\gamma^{sc}[\alpha_{s}(\mu)]\Big{)}\;,$
(122)

and this statement is true to all orders in perturbative expansion.

###  C.3 $\zeta$-prescription

The implementation of the $\zeta$-prescription leads to the definition of
optimal TMDs. We sketch here the procedure to obtain optimal TMDs referring to
the original work Scimemi:2018xaf  for further details. The anomalous
dimensions $\gamma_{F}(\mu,\zeta)$ and $\mathcal{D}(\mu,b)$ governing the
evolution can be thought as two components of a vector field in the plane
($\ln\mu^{2},\,\ln\zeta$). The integrability condition, e.g. eq. (121), states
that such field is irrotational, i.e. locally conservative. This allows to
define a scalar potential and guarantees that the evolution between two points
in the ($\ln\mu^{2},\,\ln\zeta$) space is independent of the path; in
particular, no evolution occurs along equipotential lines. However, the
perturbative expansion breaks the validity of such statement and in fact it
was shown that numerical predictions largely depend on the choice of path.
This limit is overcome by the _improved_ $\gamma$ _solution_ , that reinstates
path-invariance by supplementing $\gamma_{F}$ with formally higher-order
terms. If we let $F$ be a generic TMD, then the evolution kernel $R$,
implicitly defined as

$\displaystyle
F(x,\boldsymbol{b},\mu_{f},\zeta_{f})=R(\boldsymbol{b};\mu_{f},\zeta_{f};\mu_{i},\zeta_{i})F(x,\boldsymbol{b},\mu_{i},\zeta_{i})\,,$
(123)

within the improved $\gamma$ solution yields

$\displaystyle
R(\boldsymbol{b};\mu_{f},\zeta_{f};\mu_{i},\zeta_{i})=\exp\bigg{\\{}\mathcal{D}(\mu_{f},\boldsymbol{b})\ln\Big{(}\frac{\mu_{f}^{2}}{\zeta_{f}}\Big{)}-\mathcal{D}(\mu_{i},\boldsymbol{b})\ln\Big{(}\frac{\mu_{i}^{2}}{\zeta_{i}}\Big{)}-\int_{\mu_{i}}^{\mu_{f}}\frac{d\mu}{\mu}\big{[}2\mathcal{D}(\mu,\boldsymbol{b})+\gamma_{V}(\mu)\big{]}\bigg{\\}}\,,$
(124)

where $\gamma_{V}$ is the noncusp anomalous dimension.

Path independence allows one to apply the $\zeta$-_prescription_ , the key
point of the method. The idea is setting the initial rapidity scale
$\zeta_{i}=\zeta_{\mu_{i}}$ as a function of $\mu_{i}$ such that the scale-
dependence of the initial TMDs vanishes independent of $\mu_{i}$. At one loop,
this simply reads

$\displaystyle\zeta_{\mu}=e^{-\frac{\gamma_{V}}{\Gamma_{\mathrm{cusp}}}}\mu^{2}\,,$
(125)

and the corrections to higher loops are evaluated in Scimemi:2018xaf .

Figure 12: Sketch of the geometry of the $(\zeta,\mu)$ plane where the double
scale evolution takes place. The anomalous dimensions determine a conservative
field (grey arrows) and the evolution is null among equipotential lines
(shades of red). The intersection of two special equipotential lines (bright
red) determines a saddle point; the zeta-prescription corresponds to running
the evolution from this point, after reinstating path invariance. The
equations for the special equipotential lines in the figure correspond to the
one-loop result.

The relation between $\zeta$ and $\mu$ draws a line in the
$(\ln\zeta,\ln\mu^{2})$ plane (fig. 12). Since by requirement the TMDs are
constant along it, this must be an equipotential line, which is well defined
only if path-independence is restored. The remarkable fact with the $\zeta$
prescription is that, contrarily from standard evolution, the cancellation of
large rapidity logarithms affecting the un-evolved TMDs is an internal
mechanism. The rapidity evolution is still responsible for cancelling the
large logarithms in the hard function, but the scale uncertainty of the
evolution is now entirely decoupled from the definition of the TMDs (and in
particular, from the non-perturbative model that enters their definition).

The definition of _optimal TMDs_ requires one more specification, which
concerns the choice of initial scale $\mu_{i}$ (and consequently
$\zeta_{\mu_{i}}$), and follows from TMD factorization. Considering TMD PDFs
for definiteness, we have up to nonperturbative corrections

$\displaystyle F_{a\leftarrow
h}(x,\boldsymbol{b},\mu,\zeta_{\mu})=\sum_{b}\int_{x}^{1}\frac{dy}{y}\,\mathcal{C}_{a\leftarrow
b}(\tfrac{x}{y},\boldsymbol{b},\mu,\zeta_{\mu},\mu_{\mathrm{OPE}})f_{b\leftarrow
h}(y,\mu_{\mathrm{OPE}})\Big{[}1+\mathcal{O}\big{(}b^{2}\Lambda_{\mathrm{QCD}}^{2}\big{)}\Big{]}\,,$
(126)

where $f_{b\leftarrow h}$ are the collinear PDFs, $\mathcal{C}_{a\leftarrow
b}$ are transverse momentum matching coefficients known at two loop from ref.
Echevarria:2016scs . The matching is performed at the scale
$\mu_{\mathrm{OPE}}$. The choice of $\mu_{\mathrm{OPE}}$ is in general
constrained by $\mu_{i}$, as they need to lie on the same half-plane with
respect to the saddle point. This undesired feature is eliminated by choosing
$\mu_{i}=\mu_{\mathrm{saddle}}$.

###  C.4 Modeling of TMDPDF

At small transverse momenta non-perturbative effects inside a TMD become
dominant. A non-perturbative model valid for optimal TMDs was recently
extracted in ref. Bertone:2019nxa  by fitting combined data from Drell-Yan and
Z-boson production. Since our groomed jet functions have the same rapidity
evolution as the standard TMDs we will use the same model. First, for large
values of $b$ the initial scale, $\mu_{\mathrm{saddle}}$ enters the non-
perturbative region. We correct for this by adopting the definition

$\displaystyle\mu_{i}=\frac{2e^{-\gamma_{E}}}{b}+2\mbox{GeV}\,,$ (127)

which effectively imposes a higher cutoff on $b$. Second, the rapidity
anomalous dimension is modified as follows,

$\mathcal{D}(\mu,\boldsymbol{b})=\mathcal{D}_{\mathrm{res}}(\mu,\boldsymbol{b}^{*})+c_{0}b\,b^{*}\,,$
(128)

where $c_{0}$ is a constant, the resummed anomalous dimension can be found at
three loop in Refs. Echevarria:2012pw ; Scimemi:2018xaf , and the b-star
prescription is

$\displaystyle\boldsymbol{b}^{*}=\boldsymbol{b}\Big{(}1+\frac{b^{2}}{B_{\mathrm{NP}}^{2}}\Big{)}^{-\frac{1}{2}}\,.$
(129)

The constants $c_{0}$ and $B_{\mathrm{NP}}$ specify the nonpertubative model
in the case of $e^{+}e^{-}\to\text{2\;jets}$. For SIDIS, additional input is
required when building the TMD PDFs. Non-perturbative corrections to the
factorization formula are modeled with a multiplicative, flavor-independent
function $f_{\mathrm{NP}}$,

$\displaystyle F_{a\leftarrow
h}(x,\boldsymbol{b},\mu,\zeta_{\mu})=f_{\mathrm{NP}}(\boldsymbol{b},x)\sum_{b}\int_{x}^{1}\frac{dy}{y}\,\mathcal{C}_{a\leftarrow
b}(\tfrac{x}{y},\boldsymbol{b},\mu,\zeta_{\mu},\mu_{\mathrm{OPE}})f_{b\leftarrow
h}(y,\mu_{\mathrm{OPE}})\,,$ (130)

whose explicit expression reads

$\displaystyle
f_{\mathrm{NP}}(\boldsymbol{b},x)=\exp\bigg{\\{}-\frac{\lambda_{1}(1-x)+\lambda_{2}x+\lambda_{3}x(1-x)}{\sqrt{1+\lambda_{4}x^{\lambda_{5}}b^{2}}}b^{2}\bigg{\\}}\,,$
(131)

generalizing the common choices of gaussian or exponential functions. The five
parameters $\lambda_{i}$, together with $c_{0}$ and $B_{\mathrm{NP}}$, are
listed in table 4 of ref. Bertone:2019nxa  and were fitted within two
different schemes: the first one treats them all as free parameters, while in
the second one $B_{\mathrm{NP}}$ is fixed to 2.5 $\mbox{GeV}^{-1}$. The set of
PDF used is NNPDF3.1 Ball:2017nwa .

##  Appendix D IRC safety of the observable

It is known that the energy difference between groomed and ungroomed jet is an
IRC unsafe quantity Larkoski:2014wba ; Larkoski:2015lea . For the lepton-
antilepton annihilation process it is trivial to show even at the leading non
vanishing order, $\mathcal{O}(\alpha_{s})$. In the standard jet cross section
the collinear and soft divergence from the additional real gluon cancel
against the IR divergences form virtual contribution. This is possible because
in both the collinear and soft limits ($p_{\text{gluon}}^{0}\to 0$ and
$p^{\mu}_{\text{gluon}}\parallel p^{\mu}_{q/\bar{q}}$ ) the energy of the
(ungroomed) jet, $E_{\text{jet}}$, equals half of the center of mass energy,
$Q$, i.e. $z_{\text{jet}}\equiv 2E_{\text{jet}}/Q=1$. In the case of groomed
jets the phase space condition $p^{\mu}_{\text{gluon}}\parallel
p^{\mu}_{q/\bar{q}}$ does not guarantee $z_{\text{g-jet}}=1$ where
$z_{\text{g-jet}}$ is the energy fraction fo the groomed jet,
$z_{\text{g-jet}}\equiv 2E_{\text{g-jet}}/Q$. This is demonstrated in fig. 13
where configuration (a) corresponds to collinear gluon emission that passes
soft-drop ($p^{\mu}_{\text{gluon}}\parallel p^{\mu}_{q/\bar{q}}$ and
$p_{\text{gluon}}^{0}>z_{\text{cut}}Q$), (b) soft gluon emission
($p_{\text{gluon}}^{0}\to 0$), (c) collinear gluon emission that fails soft-
drop ($p^{\mu}_{\text{gluon}}\parallel p^{\mu}_{q/\bar{q}}$ and
$p_{\text{gluon}}^{0}<z_{\text{cut}}Q$). While divergences from the phase
space configurations (a) and (b) can cancel against the virtual divergences,
the ones in configuration (c) cannot.

Figure 13: IR divergent configurations for soft-drop groomed jets at
$\mathcal{O}(\alpha_{s})$.

This problem is usually solved when a jet substructure measurement (e.g. jet
thrust $e$) is included. In this case the configuration (c) will only
contribute to the $e=0$ bin and thus does not constitute a problem of IRC
safety, if we constrain the result for finite values of the jet thrust. In the
observable we propose, the jet thrust measurement does not help since we
require integrating over the range $e\in(0,e_{\text{cut}})$, see for example
eq. (13). However, the transverse momentum ${\boldsymbol{q}}_{T}$ does since
configuration (c) yields only ${\boldsymbol{q}}_{T}=0$ and we are interested
only in finite values of the transverse momentum. Therefore in our proposed
observable the ${\boldsymbol{q}}_{T}$ measurement regulates the IRC
divergences in a similar manner to the (differential) jet-shape measurements.

Furthermore is important to notice that the quantity
${\boldsymbol{p}}_{T\;\text{g-jet}}/z_{\text{g-jet}}$ does not directly relate
to groomed jet energy, but rather to the groomed jet direction, which is
collinear safe to all orders in $\alpha_{s}$,

$\frac{{\boldsymbol{p}}_{T\;\text{g-jet}}}{z_{\text{g-jet}}}=Q\sin(\theta_{\text{g-jet}}){\boldsymbol{n}}_{T\;\text{g-jet}}$
(132)

where $\theta_{\text{g-jet}}$ is the angle between the groomed jet axis and
the reference axis ${\boldsymbol{n}}_{T\;\text{g-jet}}$ is the direction of
the groomed jet with respect to the reference axis. Since the direction of the
groomed jet is IR safe then the observable ${\boldsymbol{q}}_{T}$ is as well.

## References

  * (1) J. Collins, _Foundations of perturbative QCD_ , _Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol._ 32 (2011) 1–624.
  * (2) M. G. Echevarria, A. Idilbi and I. Scimemi, _Factorization Theorem For Drell-Yan At Low $q_{T}$ And Transverse Momentum Distributions On-The-Light-Cone_, _JHEP_ 07 (2012) 002, [1111.4996].
  * (3) M. G. Echevarria, A. Idilbi and I. Scimemi, _Unified treatment of the QCD evolution of all (un-)polarized transverse momentum dependent functions: Collins function as a study case_ , _Phys. Rev._ D90 (2014) 014003, [1402.0869].
  * (4) J. R. Gaunt, _Glauber Gluons and Multiple Parton Interactions_ , _JHEP_ 07 (2014) 110, [1405.2080].
  * (5) A. Vladimirov, _Structure of rapidity divergences in multi-parton scattering soft factors_ , _JHEP_ 04 (2018) 045, [1707.07606].
  * (6) D. Neill, I. Scimemi and W. J. Waalewijn, _Jet axes and universal transverse-momentum-dependent fragmentation_ , _JHEP_ 04 (2017) 020, [1612.04817].
  * (7) D. Gutierrez-Reyes, I. Scimemi, W. J. Waalewijn and L. Zoppi, _Transverse momentum dependent distributions with jets_ , _Phys. Rev. Lett._ 121 (2018) 162001, [1807.07573].
  * (8) D. Gutierrez-Reyes, I. Scimemi, W. J. Waalewijn and L. Zoppi, _Transverse momentum dependent distributions in $e^{+}e^{-}$ and semi-inclusive deep-inelastic scattering using jets_, 1904.04259.
  * (9) D. Bertolini, T. Chan and J. Thaler, _Jet Observables Without Jet Algorithms_ , _JHEP_ 1404 (2014) 013, [1310.7584].
  * (10) Y. Makris, D. Neill and V. Vaidya, _Probing Transverse-Momentum Dependent Evolution With Groomed Jets_ , 1712.07653.
  * (11) Y. Makris and V. Vaidya, _Transverse Momentum Spectra at Threshold for Groomed Heavy Quark Jets_ , _JHEP_ 10 (2018) 019, [1807.09805].
  * (12) A. H. Hoang, S. Mantry, A. Pathak and I. W. Stewart, _Nonperturbative Corrections to Soft Drop Jet Mass_ , 1906.11843.
  * (13) T. Sjostrand, S. Mrenna and P. Z. Skands, _PYTHIA 6.4 Physics and Manual_ , _JHEP_ 05 (2006) 026, [hep-ph/0603175].
  * (14) T. Sjostrand, S. Mrenna and P. Z. Skands, _A Brief Introduction to PYTHIA 8.1_ , _Comput. Phys. Commun._ 178 (2008) 852–867, [0710.3820].
  * (15) A. J. Larkoski, S. Marzani, G. Soyez and J. Thaler, _Soft Drop_ , _JHEP_ 1405 (2014) 146, [1402.2657].
  * (16) S. D. Ellis and D. E. Soper, _Successive combination jet algorithm for hadron collisions_ , _Phys.Rev._ D48 (1993) 3160–3166, [hep-ph/9305266].
  * (17) S. Catani, Y. L. Dokshitzer, M. Seymour and B. Webber, _Longitudinally invariant $K_{t}$ clustering algorithms for hadron hadron collisions_, _Nucl.Phys._ B406 (1993) 187–224.
  * (18) Y. L. Dokshitzer, G. Leder, S. Moretti and B. Webber, _Better jet clustering algorithms_ , _JHEP_ 9708 (1997) 001, [hep-ph/9707323].
  * (19) M. Wobisch and T. Wengler, _Hadronization corrections to jet cross-sections in deep inelastic scattering_ , in _Monte Carlo generators for HERA physics. Proceedings, Workshop, Hamburg, Germany, 1998-1999_ , pp. 270–279, 1998.  hep-ph/9907280.
  * (20) M. Wobisch, _Measurement and QCD analysis of jet cross-sections in deep inelastic positron proton collisions at $\sqrt{s}$ = 300 GeV_.  PhD thesis, Aachen, Tech. Hochsch., 2000.
  * (21) M. Cacciari, G. P. Salam and G. Soyez, _The Anti-k(t) jet clustering algorithm_ , _JHEP_ 0804 (2008) 063, [0802.1189].
  * (22) G. Salam, _$E_{t}^{\infty}$ Scheme_, _Unpublished_ .
  * (23) A. J. Larkoski, D. Neill and J. Thaler, _Jet Shapes with the Broadening Axis_ , _JHEP_ 1404 (2014) 017, [1401.2158].
  * (24) S. Marzani, L. Schunk and G. Soyez, _A study of jet mass distributions with grooming_ , _JHEP_ 07 (2017) 132, [1704.02210].
  * (25) A. J. Larkoski, S. Marzani and J. Thaler, _Sudakov Safety in Perturbative QCD_ , _Phys.Rev._ D91 (2015) 111501, [1502.01719].
  * (26) J.-Y. Chiu, A. Jain, D. Neill and I. Z. Rothstein, _A Formalism for the Systematic Treatment of Rapidity Logarithms in Quantum Field Theory_ , _JHEP_ 1205 (2012) 084, [1202.0814].
  * (27) M. G. Echevarria, I. Scimemi and A. Vladimirov, _Universal transverse momentum dependent soft function at NNLO_ , _Phys. Rev._ D93 (2016) 054004, [1511.05590].
  * (28) Y. Li and H. X. Zhu, _Bootstrapping Rapidity Anomalous Dimensions for Transverse-Momentum Resummation_ , _Phys. Rev. Lett._ 118 (2017) 022004, [1604.01404].
  * (29) A. A. Vladimirov, _Correspondence between Soft and Rapidity Anomalous Dimensions_ , _Phys. Rev. Lett._ 118 (2017) 062001, [1610.05791].
  * (30) I. Scimemi and A. Vladimirov, _Power corrections and renormalons in Transverse Momentum Distributions_ , _JHEP_ 03 (2017) 002, [1609.06047].
  * (31) M. G. Echevarria, I. Scimemi and A. Vladimirov, _Unpolarized Transverse Momentum Dependent Parton Distribution and Fragmentation Functions at next-to-next-to-leading order_ , _JHEP_ 09 (2016) 004, [1604.07869].
  * (32) M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J._ C72 (2012) 1896, [1111.6097].
  * (33) I. Scimemi and A. Vladimirov, _Systematic analysis of double-scale evolution_ , _JHEP_ 08 (2018) 003, [1803.11089].
  * (34) T. Gehrmann, T. Luebbert and L. L. Yang, _Calculation of the transverse parton distribution functions at next-to-next-to-leading order_ , _JHEP_ 06 (2014) 155, [1403.6451].
  * (35) U. D’Alesio, M. G. Echevarria, S. Melis and I. Scimemi, _Non-perturbative QCD effects in $q_{T}$ spectra of Drell-Yan and Z-boson production_, _JHEP_ 11 (2014) 098, [1407.3311].
  * (36) M. G. Echevarria, A. Idilbi, A. Schäfer and I. Scimemi, _Model-Independent Evolution of Transverse Momentum Dependent Distribution Functions (TMDs) at NNLL_ , _Eur. Phys. J._ C73 (2013) 2636, [1208.1281].
  * (37) J.-y. Chiu, A. Jain, D. Neill and I. Z. Rothstein, _The Rapidity Renormalization Group_ , _Phys. Rev. Lett._ 108 (2012) 151601, [1104.0881].
  * (38) S. Moch, J. A. M. Vermaseren and A. Vogt, _The Three loop splitting functions in QCD: The Nonsinglet case_ , _Nucl. Phys._ B688 (2004) 101–134, [hep-ph/0403192].
  * (39) S. Moch, J. A. M. Vermaseren and A. Vogt, _Three-loop results for quark and gluon form-factors_ , _Phys. Lett._ B625 (2005) 245–252, [hep-ph/0508055].
  * (40) P. A. Baikov, K. G. Chetyrkin, A. V. Smirnov, V. A. Smirnov and M. Steinhauser, _Quark and gluon form factors to three loops_ , _Phys. Rev. Lett._ 102 (2009) 212002, [0902.3519].
  * (41) A. Vogt, F. Herzog, S. Moch, B. Ruijl, T. Ueda and J. A. M. Vermaseren, _Anomalous dimensions and splitting functions beyond the next-to-next-to-leading order_ , _PoS_ LL2018 (2018) 050, [1808.08981].
  * (42) C. W. Bauer, F. J. Tackmann, J. R. Walsh and S. Zuberi, _Factorization and Resummation for Dijet Invariant Mass Spectra_ , _Phys.Rev._ D85 (2012) 074006, [1106.6047].
  * (43) C. Frye, A. J. Larkoski, M. D. Schwartz and K. Yan, _Factorization for groomed jet substructure beyond the next-to-leading logarithm_ , _JHEP_ 07 (2016) 064, [1603.09338].
  * (44) D. Kang, C. Lee and I. W. Stewart, _Using 1-Jettiness to Measure 2 Jets in DIS 3 Ways_ , _Phys. Rev._ D88 (2013) 054004, [1303.6952].
  * (45) V. Bertone, I. Scimemi and A. Vladimirov, _Extraction of unpolarized quark transverse momentum dependent parton distributions from Drell-Yan/Z-boson production_ , 1902.08474.
  * (46) M. Beneke and V. M. Braun, _Power corrections and renormalons in Drell-Yan production_ , _Nucl. Phys._ B454 (1995) 253–290, [hep-ph/9506452].
  * (47) G. P. Korchemsky and G. F. Sterman, _Nonperturbative corrections in resummed cross-sections_ , _Nucl. Phys._ B437 (1995) 415–432, [hep-ph/9411211].
  * (48) G. P. Korchemsky, G. Oderda and G. F. Sterman, _Power corrections and nonlocal operators_ , _AIP Conf. Proc._ 407 (1997) 988, [hep-ph/9708346].
  * (49) M. Beneke, V. M. Braun and L. Magnea, _Phenomenology of power corrections in fragmentation processes in e+ e- annihilation_ , _Nucl. Phys._ B497 (1997) 297–333, [hep-ph/9701309].
  * (50) T. Becher and G. Bell, _Enhanced nonperturbative effects through the collinear anomaly_ , _Phys. Rev. Lett._ 112 (2014) 182002, [1312.5327].
  * (51) “artemide web-page, https://teorica.fis.ucm.es/artemide/
artemide repository, https://github.com/vladimirovalexey/artemide-public.”

  * (52) I. Scimemi and A. Vladimirov, _Analysis of vector boson production within TMD factorization_ , _Eur. Phys. J._ C78 (2018) 89, [1706.01473].
  * (53) NNPDF collaboration, R. D. Ball et al., _Parton distributions from high-precision collider data_ , _Eur. Phys. J._ C77 (2017) 663, [1706.00428].