Coarse-grained binning in Drell-Yan transverse momentum spectra

Measurements of transverse momentum spectra of electroweak bosons via Drell-Yan (DY) lepton-pair production [1] are at the core of many aspects of the physics program at the Large Hadron Collider (LHC), ranging from precision determinations of Standard Model parameters in the strong and electroweak sectors, to non-perturbative features of hadron structure and transverse momentum dependent (TMD) parton distribution functions [2].

Recent applications of DY transverse momentum measurements include the determination of the strong coupling from perturbative Quantum Chromodynamics (QCD) predictions matched with resummation [3, 4]; the extraction of TMD parton distributions both in analytic-resummation [5, 6, 7, 8] and parton-branching [9] approaches; the modeling of TMD contributions in soft-collinear effective theory [10]; the determination of intrinsic-$k_{T}$ parameters in the tuning of Monte Carlo event generators [11].

A common thread in these applications is that high-precision measurements of the low transverse-momentum region ($\Lambda_{\rm{QCD}}{\raisebox{-2.58334pt}{\hbox to0.0pt{$\,\sim\,$\hss}}\raisebox{1.72218pt}{$\,<\,$}}p_{T}(ll)\ll m(ll)$, where $p_{T}$ and $m$ are the dilepton’s transverse momentum and invariant mass) with fine binning in $p_{T}$ will improve our ability to unravel QCD dynamics, involving multiple soft-gluon emissions as well as the non-perturbative intrinsic transverse motion of partons. See for instance the role of fine-binned $p_{T}(ll)$ measurements in the context of analytic-resummation studies [12] and parton-branching studies [13].

However, recent measurements of the DY process by the ATLAS [14] and CMS [15] collaborations indicate that such fine-binned treatments require an extremely delicate control of systematic uncertainties. Experimentally, the determination of the low-$p_{T}(ll)$ fine-binned structure is challenging. Under these circumstances, it becomes important to assess the capabilities of methodologies which do not require the fine-binned measurement of $p_{T}(ll)$ distributions, and rather rely on coarse-grained binning.

In this paper we explore one such methodology, by studying the intrinsic $k_{T}$ determination in TMD parton distributions from the measurement of $pp\to Z/\gamma\to l^{+}l^{-}$ in the region $0<p_{T}(ll)<p_{T,{\rm{max}}}$ (with $p_{T,{\rm{max}}}$ small compared to $m(ll)$) based on coarse-grained partitions of this region. We subdivide the $p_{T}$ region into two bins with separation momentum $p_{s}$, $p_{T}<p_{s}$ and $p_{T}>p_{s}$, and investigate the sensitivity to the intrinsic $k_{T}$ distribution from measuring the ratio between the low $p_{T}(ll)$ cross section (predominantly sensitive to TMD dynamics and resummed soft-parton radiation) and the high $p_{T}(ll)$ cross section (predominantly sensitive to fixed-order hard-parton radiation) as a function of varying $p_{s}$. We study the role of bin-to-bin migration effects on the $p_{T}$ ratio. Although the high $p_{T}(ll)$ cross section is not sensitive to intrinsic $k_{T}$, it acts as a reference in the relative ratio between the strength of the intrinsic $k_{T}$ and fixed-order contribution. We find that the intrinsic $k_{T}$ can be determined by measuring its overall strength through the $p_{T}$-ratio instead of measuring the low $p_{T}(ll)$ differential cross section. This conclusion is relevant, from a practical viewpoint, in order to achieve a good determination of intrinsic $k_{T}$ from analyses of experimental data with reduced systematic uncertainties. We illustrate this methodology by presenting an extraction of intrinsic $k_{T}$ from DY $p_{T}(ll)$ measurements by the CMS collaboration [16].

To perform this study, we use the parton branching (PB) approach [17, 18] to TMD evolution. This approach provides the basis for including TMD distributions in parton shower Monte Carlo calculations, and has been shown to give a consistent description of DY $p_{T}$ spectra [9] as well as of the multiple-jet structure associated with DY production [19]. At inclusive level, it can be related to the evolution of collinear parton distribution functions (PDF), and correctly describes deep-inelastic scattering structure functions [20]. Given the wide applicability of the PB TMD approach, this provides a useful framework to assess the performance of the $p_{T}$-ratio methodology.

The paper is organized as follows. In Sect. 2, we briefly review the PB TMD approach. In Sect. 3, we propose the $p_{T}$-ratio as a new observable and methodology for studies of TMD dynamics. The sensitivity of the $p_{T}$-ratio to the intrinsic $k_{T}$ distribution is investigated with a pseudo-data sample, and is compared to that from the fine-binning structure of the low-$p_{T}(ll)$ differential cross section. In Sect. 4, we apply the methodology to the $p_{T}(ll)$ distribution measured by the CMS collaboration [16] at 13 TeV across a wide range in DY mass from 50 GeV to 1 TeV. We give conclusions in Sect. 5.

In this section we briefly describe the main features of the PB approach [17, 18, 20] which will be used for the analysis in this work. We first recall the PB evolution equations; then we discuss the intrinsic-$k_{T}$ distribution, the soft-gluon resolution scale, the treatment of the strong coupling, and the application of the method to the computation of DY $p_{T}$ distributions.

We consider the TMD parton distribution of parton flavor $a$, $A_{a}(x,{\bf k},\mu^{2})$, as a function of the longitudinal momentum fraction $x$, the transverse momentum ${\bf k}$ and the evolution scale $\mu$. According to Refs. [18, 20], the distributions $A_{a}(x,{\bf k},\mu^{2})$ fulfill evolution equations of the form

where $\Delta_{a}(\mu^{2},\mu_{0}^{2})$ is the Sudakov form factor and ${\cal E}_{ab}$ are the evolution kernels, which are specified in Ref. [18] as functionals of the Sudakov form factors $\Delta_{a}$, of the real-emission splitting functions $P_{ab}^{R}$, and of phase space constraints collectively denoted by $\Theta$ in Eq. (2.1). The functions that appear in the evolution kernels are perturbatively computable as power series expansions in the strong coupling $\alpha_{s}$. The explicit expressions of these expansions for all flavor channels are given to two-loop order in Ref. [18].

The evolution in Eq. (2.1) is expressed in terms of two branching variables, $z$ and ${\boldsymbol{\mu}}^{\prime}$: $z$ is the longitudinal momentum transfer at the branching, controlling the rapidity of the parton emitted along the branching chain; $\mu^{\prime}=\sqrt{{\boldsymbol{\mu}}^{\prime 2}}$ is the mass scale at which the branching occurs, and is related to the branching’s kinematic variables according to the ordering condition. It is worth noting that the double evolution in rapidity and mass for TMD distributions can also be formulated in a CSS [21, 22], rather than PB, approach — see e.g. [23, 12, 24]. A study of the relationship of the PB Sudakov evolution kernel with the CSS formalism is performed in [25, 26]. Eq. (2.1) is obtained for the case of angular ordering [27, 28, 29], which gives $\mu^{\prime}=q_{\perp}/(1-z)$, where $q_{\perp}$ is the emitted parton’s transverse momentum. The angular-ordered branching is motivated by the treatment of the endpoint region in the TMD case [18, 30].

The initial evolution scale in Eq. (2.1) is denoted by $\mu_{0}$ ($\mu_{0}>\Lambda_{\rm{QCD}}$), and is usually taken to be of order 1 GeV. The distribution ${A}_{a}(x,{\bf k},\mu_{0}^{2})$ at scale $\mu_{0}$ in the first term on the right hand side of Eq. (2.1) provides the intrinsic $k_{T}$ distribution. This is a nonperturbative boundary condition to the evolution equation, and may be determined from comparisons of theory predictions to experimental data. In the calculations of this paper, for simplicity we will parameterize ${{\cal A}}_{a}(x,{\bf k},\mu^{2}_{0})$, following previous applications of the PB TMD method [20, 9], as

with the width of the Gaussian distribution given by $\sigma=q_{s}/\sqrt{2}$, independent of parton flavor and $x$, where $q_{s}$ is the intrinsic-$k_{T}$ parameter.

The evolution kernels ${\cal E}_{ab}$ and Sudakov form factors $\Delta_{a}$ in Eq. (2.1) contain kinematic constraints which embody the phase space of the branchings along the parton cascade. These can be described, using the “unitarity” picture of QCD evolution [31], by separating resolvable and non-resolvable branchings in terms of a resolution scale to classify soft-gluon emissions [17]. (See e.g. [32] for a study of the interplay of resolution scale with transverse momentum recoils in parton showers.) Applications of the PB TMD method can be carried out either by taking the soft-gluon resolution scale to be a fixed constant value close to the kinematic limit $z=1$ or by allowing for a running, ${\mu}^{\prime}$-dependent resolution (see, e.g., discussions in [29, 33, 9]). In the computations of this paper, we will take fixed resolution scale. This is the same as what is done in the studies of Refs. [20, 13, 9], which we will use as benchmarks.

The scale at which the strong coupling $\alpha_{s}$ is to be evaluated in Eq. (2.1) is a function of the branching variable. Two scenarios are commonly studied in PB TMD applications: i) $\alpha_{s}=\alpha_{s}({\mu}^{\prime 2})$; ii) $\alpha_{s}=\alpha_{s}(q_{\perp}^{2})=\alpha_{s}({\mu}^{\prime 2}(1-z)^{2})$. Case i) corresponds to DGLAP evolution [34, 35, 36] , while case ii) corresponds to angular ordering, e.g. CMW evolution [27, 28]. In Ref. [20], fits to precision deep inelastic scattering HERA data [37] are performed for both scenarios i) and ii), using the fitting platform xFitter [38, 39]. It is found that fits with good $\chi^{2}$ values can be achieved in either case. Correspondingly, PB-NLO-2018 Set1 (with the DGLAP-type $\alpha_{s}({\bf q}^{\prime 2})$) and PB-NLO-2018 Set2 (with the angular-ordered CMW-type $\alpha_{s}(q_{T}^{2})$) are obtained.

On the other hand, it is found that PB-NLO-2018 Set2 provides a much better description, compared to PB-NLO-2018 Set1, of measured $Z/\gamma$ $p_{T}$ spectra at the LHC [13] and in low-energy experiments [40], and of di-jet azimuthal correlations near the back-to-back region at the LHC [41, 42]. Further, it is shown that a good description is also obtained for DY + jets final-state distributions [43, 44, 19, 45]. We will therefore base the analysis of this work on PB-NLO-2018 Set2. As in [20, 9], the strong coupling is modeled according to a “pre-confinement” picture [46, 47] as $\alpha_{s}=\alpha_{s}(\max(q^{2}_{c},{\bf q}_{\perp}^{2}))$, where $q_{c}$ is a semi-hard scale, taken to be $q_{c}=1$ GeV. The PB TMD sets are available from the TMDlib library [48, 49].

The calculation of DY production cross sections in the PB TMD method proceeds as described in Ref. [13]. NLO hard-scattering matrix elements are obtained from the MadGraph5_aMC@NLO [50] (hereafter, MCatNLO) event generator and matched with TMD parton distributions and showers obtained from PB evolution [20, 18, 17] and implemented in the Cascade Monte Carlo event generator [51, 52], using the subtractive matching procedure proposed in [13] and further analyzed in [45]. In particular, the Herwig6 [53, 54] subtraction terms are used in MCatNLO, as they are based on the same angular ordering conditions as the PB TMD distributions [45]. Final state parton showers are generated from Pythia6.428 [55], including photon radiation.

Ref. [9] uses the approach described above to compute theoretical results for DY transverse momentum distributions, and to make a determination of the intrinsic-$k_{T}$ parameter $q_{s}$ in Eq. (2) by performing fits of the results to DY experimental measurements of DY $p_{T}(ll)$ [16, 56, 57, 58, 59, 60, 61, 62, 63]. The study [9] includes a detailed treatment of statistical, correlated and uncorrelated uncertainties. The CMS data [16] at center of mass energy $\sqrt{s}=$ 13 TeV cover DY invariant masses from 50 GeV to 1 TeV. The other data cover energies from 13 TeV down to 38 GeV. With the fitted $q_{s}$ values [9], PB TMD distributions are thus obtained. Further discussions of the PB TMD distributions [9] are given in [64, 65, 66, 67].

In the next section we will use Monte Carlo predictions from the PB TMD method to explore the feasibility of extracting intrinsic-$k_{T}$ distributions from experimental measurements with coarse-grained binning in the low $p_{T}(ll)$ region.

In this section we present the $p_{T}$-ratio observable and methodology for studies of TMD dynamics. We will focus on the process $pp\to Z/\gamma\to l^{+}l^{-}$ in the transverse momentum region $0<p_{T}(ll)<p_{T,{\rm{max}}}$, where $p_{T,{\rm{max}}}$ is taken to be small compared to the invariant mass $m(ll)$.

We construct PB TMD predictions for the DY $p_{T}(ll)$ distribution, as described in the previous section, from MCatNLO + Cascade Monte Carlo, with initial-scale TMD distributions given in Eq. (2). To explore the intrinsic-$k_{T}$ parameter space, we generate 7 template samples with different $q_{s}$ values using the TMD grid files available on the TMDlib [48] website: these are $q_{s}$ = 0.5 GeV (using the PB TMD set PB-NLO-HERAI+II-2018-set2 [20]), $q_{s}$ = 0.59, 0.74, 0.89 GeV (using the set PB-NLO-HERAI+II-2023-set2-qs=0.74 [9]), and $q_{s}$ = 0.96, 1.04, 1.12 GeV (using the set PB-NLO-HERAI+II-2023-set2-qs=1.04 [9]). We study the influence of the intrinsic-$k_{T}$ width by examining the templates firstly in the full phase space (later on we will consider an experimental fiducial phase space). Results are shown in Fig. 1 over the range $p_{T}(ll)<20$ GeV.

The $p_{T}(ll)$ shape in the rising part of the spectrum is altered by $q_{s}$, with the peak position shifting to higher $p_{T}$ when $q_{s}$ increases. However, the trend is mild, suggesting that the essential information on the intrinsic $k_{T}$ width may be extracted from a well-defined ratio between the low and high $p_{T}$ regions. To this end, we introduce a momentum parameter $p_{s}$ to separate the regions $p_{T}<p_{s}$ and $p_{T}>p_{s}$, and construct the 2-bin $p_{T}$-ratio between the lower and higher $p_{T}$ regions. The $p_{T}$-ratio is defined from the integral of event numbers in the relatively low ($p_{L}$) and high ($p_{H}$) region, i.e., from the two-binned $p_{T}(ll)$,

The uncertainty of the $p_{T}$-ratio is obtained from propagating statistical uncertainties in the low and high $p_{T}$ bins, and is given by

The 2-bin distribution and corresponding $p_{T}$-ratio are plotted, for the case of a given value $p_{s}=3$ GeV, in Fig. 2. The fall-off of the $p_{T}$-ratio with increasing $q_{s}$ indicates that one may be able to extract TMD parameters from this ratio instead of the complete fine-binned $p_{T}$ shape, which results into a significant advantage from the standpoint of controlling experimental systematic uncertainties.

The separation momentum $p_{s}$ will play a critical role in the $q_{s}$ extraction from the $p_{T}$-ratio. TMD sensitivity arises primarily from the low-$p_{T}$ region. If $p_{s}$ is too low inside this region, a proportion of soft-gluon splitting transverse momenta will be integrated into the upper region, causing a loss of sensitivity to TMD parameter determination. If, on the other hand, $p_{s}$ is higher than the peak, the information on low $p_{T}$ will tend to become undetectable.

We next perform a sensitivity test on both the fine-binned $p_{T}$ and the $p_{T}$-ratio. We test the former with different bin width from 0.5 GeV to 3 GeV, to examine the potential of the fine-binned $p_{T}$ shape in extracting $q_{s}$. We test the latter, on the other hand, focusing on the sensitivity change due to different choice of $p_{s}$. We shift $p_{s}$ from 1.5 GeV to 4 GeV to find an optimal definition of the $p_{T}$-ratio, as well as to verify whether it has statistically the same sensitivity as the fine-binned $p_{T}$ shape. The $p_{T}(ll)$ greater than 10 GeV is expected to be ancillary in TMD performance studies. In this test, we will consider $p_{T,{\rm{max}}}=$ 10 GeV, 20 GeV, in a similar spirit to previous studies of intrinsic-$k_{T}$ distributions (see e.g. Refs. [8, 9, 12]).

In this test, we choose the sample with $q_{s}=0.89$ as the pseudo-data, and extract $q_{s}$ with the PB TMD replicas in the phase space with final-state lepton selections $p_{T}(l^{\pm})>25$ GeV and $|\eta(l^{\pm})|<2.4$, in which leptons are dressed with photons reconstructed within $\Delta R=\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}<0.1$. This is close to the fiducial phase space in real experiments at ATLAS and CMS. Sensitivity is represented by the fitting uncertainty of $q_{s}$. Each sample in PB TMD template contains 100 million events. The fitting is performed with least square method, in which the estimator $\chi^{2}$ is defined as

with $m_{i}$ and $\mu_{i}$ being the observed and predicted data point, respectively, and $C_{ij}$ the covariance matrix. We consider uncorrelated statistical uncertainties from the pseudo-data and template as the only sources contributing to $C_{ij}$, making it a diagonal matrix. The fitting uncertainties in both cases are depicted in Fig. 3. The result indicates that the $p_{T}$-ratio is as sensitive as the fine-binned $p_{T}$ shape.

In real experiments, momentum resolution causes bin-to-bin migration. Since the migration effect is a dominant systematic uncertainty, it is essential to verify that such effect on the $p_{T}$-ratio is not destructive. To study this, we apply a 3% resolution on the four-momentum of the dressed leptons as

where $g$ follows a Gaussian distribution $\mathcal{G}(0,0.03)$. The migration effect is illustrated in Fig. 4. It moves the peak position to a higher $p_{T}$ value, which is a similar trend as the $q_{s}$ effect. Also, owing to this the optimal separation choice $p_{s}$ is different from the non-migration scenario. If $p_{s}$ is defined correctly in this case, we expect the behavior to be the same as it is at truth level. We also note that in the bottom plot in Fig. 4 the relative ratio for $p_{T}(ll)<1$ GeV is smaller at the reco-level. This should result in a mild decrease in sensitivity.

We now repeat the sensitivity test after the 3% energy resolution is applied to the dressed leptons. The results are shown in Fig. 5. The uncertainties are slightly larger than those in Fig. 3. Also, by comparing the fitting uncertainties in the cases of fine-binned $p_{T}(ll)$ and $p_{T}$-ratio, we see that the decrease in sensitivity due to the migration effect is similar. We are led to conclude that the $p_{T}$-ratio still has a good sensitivity to the intrinsic $k_{T}$, so that this methodology can be applied directly in the experiments for TMD performance studies.

Having performed a test with pseudo-data in the previous section, in this section we move on to consider a real experimental data analysis.

DY $p_{T}(ll)$ spectra are measured at the LHC with bin widths equal to or greater than 1 GeV. Measuring the fine-binned $p_{T}(ll)$ structure, on the other hand, is challenging. Sizeable systematic uncertainties have been reported in recent DY differential cross section measurements [14, 15], with luminosity of approximately 36 $fb^{-1}$ data collected in 2016, especially for the $p_{T}(ll)<2$ GeV bins. This is due to the uncertainties mostly induced from lepton efficiency and momentum resolution. Moreover, since the bin-to-bin correlation of the unfolding method uncertainties and momentum migration uncertainties are negative, increasing the bin numbers in such regions inevitably increases these systematic uncertainties.

Given these difficulties in achieving finer binning in the very low $p_{T}(ll)$ range, we next test the $p_{T}$-ratio methodology of the previous section, based on less fine binning in the low-$p_{T}$ fiducial region, by applying it to real experimental data. We emphasize that the purpose of this study is not to carry out a precision extraction of $q_{s}$, but to perform a feasibility study for the 2-bin $p_{T}$-ratio in the context of TMD parameters determination.

The CMS collaboration has measured the $p_{T}(ll)$ distribution in a wide mass range from 50 GeV to 1 TeV, with a bin width of 1 GeV in the DY invariant mass bin around the $Z$-boson mass [16]. This measurement has been used for the determination of the $q_{s}$ dependence on DY mass in Ref. [9]. We here perform an extraction of $q_{s}$ from this measurement using the $p_{T}$-ratio, and compare the results with those in Ref. [9]. Besides testing the feasibility of the $p_{T}$-ratio proposal with real data, we also aim to check that the $p_{T}$-ratio is not biased by the high-$p_{T}$ region, where contributions from higher jet multiplicities are important [43, 19].

The selections of the leptons are

• •

  $|\eta|<2.4$ for both lepton,

• •

  $p_{T}>25$ GeV for leading lepton,

• •

  $p_{T}>20$ GeV for subleading lepton.

We use dressed-level leptons, defined with $\Delta R<$ 0.1. The scale uncertainties are treated as fully correlated, and computed by varying the renormalization and factorization scales by a factor of 2 separately, excluding the two extreme cases. Since the model-induced statistical uncertainty is negligible comparing to the experimental uncertainty, the final covariance matrix $C_{ij}$ comprises total covariance matrix from experiment and scale uncertainty:

The description of the CMS data by the theoretical predictions is shown in Fig. 6. In our study, the $p_{T}(ll)$ range is the same as in Ref. [9] – 6 GeV for the first $m_{ll}$ bin, 7 GeV for the second bin, and 8 GeV for the rest. This limitation prevents the ratio from being biased by the difference between prediction and data in the higher $p_{T}$ region. The choice of $p_{s}$ is largely restricted by the measurement, which affects the final sensitivity, as is suggested in the sensitivity study in Sect. 3. We stress that, if the $p_{T}$-ratio is used for TMD parameter extractions, the choice of the $p_{s}$ value will need careful investigation. For our test in this paper, $p_{s}$ is set to 2 GeV for all $m_{DY}$ bins.

Final results are reported in Fig. 7. For comparison, the statistical uncertainties and scale uncertainties from data are considered, and compared to those in the analysis [9]. The uncertainties labeled with data at 68% confidence level are obtained from the $q_{s}$ gap defined by $\chi^{2}=\chi^{2}_{min}+1$, and should be compared with the uncertainties labeled with data in Ref. [9]. Another source of uncertainties comes from the choice of the $p_{T}$ range. In Ref. [9], this is estimated by varying the numbers of bins in the fit. In this study, we obtain it in the same way. For the last two $m_{DY}$ bins, it is not estimated for the lack of statistics. The result in every $m_{DY}$ region is consistent with the previous result [9].

In this paper we propose the $p_{T}(ll)$-ratio, which is defined as the ratio of low and high transverse momentum region of DY lepton pairs, as a sensitive observable and methodology for the extraction of a TMD parameter, the intrinsic $k_{T}$ Gaussian width $q_{s}$. We firstly review the basic idea of the parton branching method for TMD evolution, and the calculation of $p_{T}(ll)$ in the DY process. We then point out that the sensitivity of the $p_{T}(ll)$ shape to the intrinsic $k_{T}$ distribution actually lies in the shift of the spectrum from low $p_{T}$ to high $p_{T}$. This observation leads to the proposal of the $p_{T}(ll)$-ratio. We numerically test the sensitivity of the $p_{T}$-ratio to $q_{s}$ by comparing its statistical uncertainty with the same procedure on fine-binning $p_{T}$ structure. The result shows that the $p_{T}$-ratio has comparable sensitivity to the fine-binned $p_{T}$ shape.

The precision in the most sensitive $p_{T}(ll)$ regions in real measurements at the LHC is largely restricted by experimental systematic uncertainties. These make the fine-binned $p_{T}$ structure difficult to access experimentally, and point to the need to explore approaches based on coarse-grained binning. The relative ratio of low and high $p_{T}(ll)$ regions avoids the fine binning. In this paper, we illustrate its effectiveness by extracting $q_{s}$ from the $p_{T}$-ratio in a wide invariant-mass range of DY lepton pairs measured by the CMS collaboration recently. Since the CMS $p_{T}$ is binned into 2 GeV (or 1 GeV), we can only conduct this test by rebinning the distributions and covariance matrices, which means we cannot arbitrarily choose the separation momentum $p_{s}$. Nevertheless, the result is consistent with previous results obtained from the low-$p_{T}(ll)$ distribution of DY lepton pairs as a function of invariant mass, and has the same level of precision.

From this perspective, this methodology can be tested in future studies by using more complex TMD parameterization forms, or investigating extreme conditions of very small or very large $q_{s}$ values, to make sure the extraction is unbiased. It could be used in forthcoming experiments for TMD performance studies, since it avoids the need to control large systematics in the low $p_{T}(ll)$ region. Unlike the extraction from a $p_{T}$ distribution that has been already binned, the $p_{T}$-ratio requires a dedicated study of the separation momentum $p_{s}$, affecting the sensitivity to TMD parameters, for example along the lines of the pseudo-data test carried out in this work.

We are grateful to Hannes Jung for discussion and for assistance with the use of PB TMD templates. We thank Laurent Favart, Jibo He, Hengne Li, Louis Moureaux and Hang Yin for useful conversations. This work was supported by the National Natural Science Foundation of China under Grants No.11721505, No.12061141005, and No.12105275, and supported by the “USTC Research Funds of the Double First-Class Initiative”.