aainstitutetext: INPAC, Key Laboratory for Particle Astrophysics and Cosmology
(MOE), Shanghai Key Laboratory for Particle Physics and Cosmology, School of
Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240,
Chinabbinstitutetext: Department of Physics, University of Maryland, College
Park, MD 20742

# Hybrid Renormalization for Quasi Distribution Amplitudes of A Light Baryon

Chao Han b,* Yushan Su††*Corresponding author a Wei Wang a Jia-Lu Zhang
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We develop a hybrid scheme to renormalize quasi distribution amplitudes of a
light baryon on the lattice, which combines the self-renormalization and ratio
scheme. By employing self-renormalization, the UV divergences and linear
divergence at large spatial separations in quasi distribution amplitudes are
removed without introducing extra nonperturbative effects, while making a
ratio with respect to the zero-momentum matrix element can properly remove the
UV divergences in small spatial separations. As a specific application,
distribution amplitudes of the $\Lambda$ baryon made of $uds$ are
investigated, and the requisite equal-time correlators, which define quasi
distribution amplitudes in coordinate space, are perturbatively calculated up
to the next-to-leading order in strong coupling constant $\alpha_{s}$. These
perturbative equal-time correlators are used to convert lattice QCD matrix
elements to the continuum space during the renormalization process.
Subsequently, quasi distribution amplitudes are matched onto lightcone
distribution amplitudes by integrating out hard modes and the corresponding
hard kernels are derived up to next-to-leading order in $\alpha_{s}$ including
the hybrid counterterms. These results are valuable in the lattice-based
investigation of the lightcone distribution amplitudes of a light baryon from
the first principles of QCD.

##  1 Introduction

Lightcone distribution amplitudes (LCDAs) of light baryons are the fundamental
non-perturbative inputs in QCD factorization for exclusive processes with a
large momentum transfer Shih:1998pb . An example of this type is weak decays
of bottom baryons which are valuable to extract the CKM matrix element
$|V_{ub}|$ LHCb:2015eia  and to probe new physics beyond the standard model
through flavor changing neutral current process LHCb:2015tgy ; LHCb:2021byf .
In addition, the knowledge of LCDAs is also crucial for understanding the
internal structure of baryons. The LCDAs characterize the distributions of the
longitudinal momentum among the quarks and gluons within the dominant leading
Fock state of a baryon, which are complementary to parton distribution
functions that encode the probability density distribution of parton momenta
in hadrons. Thus the LCDAs of light baryons have been widely investigated by
theoretical techniques such as QCD sum rules Chernyak:1984bm ; King:1986wi ;
Braun:1999te ; Anikin:2013aka  and Lattice QCD Gockeler:2008xv ; QCDSF:2008qtn
; Bali:2015ykx ; RQCD:2019hps .

Despite its pivotal importance, LCDAs of light mesons and baryons have been
less studied compared to parton distribution functions (PDFs). A major
difficulty is that in an exclusive process it is very likely that more than
one LCDAs enter into one physical observable in a rather complicated way
through convolution integrals. Thus it makes an experimental determination of
LCDAs extremely difficult. In addition, the leading-twist baryon LCDA
$\Phi\left(y_{1},y_{2},\mu\right)$ with $y_{1},y_{2}$ being the momentum
fractions of two involved quarks and the momentum fraction for the third quark
satisfying $y_{3}=1-y_{1}-y_{2}$ describes the momentum distributions of the
three quarks, and by definition is a two-dimensional distribution function,
which is even more complicated than the meson LCDAs. Our limited knowledge of
baryon LCDAs mostly relies on non-perturbative methods, many of which are
model-dependent and inevitably introduce uncontrollable uncertainties. In
recent years lattice QCD has been applied to determine the normalization
constants and the ﬁrst moments of the distribution amplitudes for the lowest-
lying baryon octet QCDSF:2008qtn ; Bali:2015ykx ; RQCD:2019hps . In
particular, a latest lattice QCD study has used a large number of $n_{f}=2+1$
ensembles with physical pion (and kaon) masses and ﬁve diﬀerent lattice
spacings RQCD:2019hps . After making the extrapolation to the continuum and
inﬁnite volume limits, they obtained results for the first two moments of
LCDAs. Despite these progresses, a complete description of LCDAs can not be
constructed from these few moments and the LCDA is far from being deciphered.

In a previous publication Deng:2023csv , a direct method to extract the shape
distribution of LCDA of a light baryon is proposed through the simulation of
equal-time correlation functions, named as quasi distribution amplitude (DA),
under the framework of large momentum effective theory (LaMET) Ji:2013dva ;
Ji:2014gla  (please see Refs. Cichy:2018mum ; Zhao:2018fyu ; Ji:2020ect  for
reviews on the development and successful applications in LaMET). The quasi-DA
$\tilde{\Phi}(x_{1},x_{2},P^{z},\mu)$ with $P^{z}$ being the hadron momentum
on the $z$ direction and $x_{1},x_{2}$ being the momentum fractions is a
calculable quantity on the lattice. Since the quasi-DA and LCDA have the same
infrared (IR) structure, the LCDA of a light baryon can be obtained by
performing a boost on the hadron momentum to infinity, which is captured by
the matching formula

$\tilde{\Phi}\left(x_{1},x_{2},P^{z},\mu\right)=\int
dy_{1}dy_{2}\mathcal{C}\left(x_{1},x_{2},y_{1},y_{2},P^{z},\mu\right)\Phi\left(y_{1},y_{2},\mu\right)+\mathcal{O}\left(\frac{\Lambda_{\rm
QCD}}{x_{1}P^{z}},\frac{\Lambda_{\rm QCD}}{x_{2}P^{z}},\frac{\Lambda_{\rm
QCD}}{\left(1-x_{1}-x_{2}\right)P^{z}}\right),$ (1)

where $\mathcal{C}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)$ is a hard kernel to
compensate for the ultraviolet (UV) differences between these two
distributions. The hard kernel has been calculated up to one-loop accuracy in
the $\overline{\rm MS}$ scheme Deng:2023csv . The $\mu$ in the quasi-DA comes
from the renormalization of the logarithmic divergences, while the scale $\mu$
in the LCDA is the factorization scale to split the collinear and hard modes.
Thus the hard kernel $\mathcal{C}$ contains both renormalization and
factorization scales, which are chosen to be the same for convenience. To
remove the remnant ultraviolet divergence in the quasi-DA, a regularization
invariant momentum subtraction method (RI/MOM) Martinelli:1994ty  was adopted
to renormalize the quasi-DA and the corresponding one-loop counterterm was
obtained Deng:2023csv .

Despite the theoretical advantages of the RI/MOM scheme, discrepancies emerge
when implementing this scheme to renormalize the quasi-PDFs and quasi-DAs on
the lattice. Residual linear divergence continues to manifest Zhang:2020rsx
and additional uncontrollable infrared effects are avoidably introduced
LatticePartonCollaborationLPC:2021xdx . Moreover, the discretization effects
arising in the lattice calculations also need to be considered carefully.
Given all these difficulties, it remains a challenge to perform a practical
lattice calculation for baryon quasi-DA.

In order to address these issues, we develop a hybrid renormalization method
following the same spirit of Ref. Ji:2020brr . The concept of the hybrid
renormalization scheme was first proposed in Ref. Ji:2020brr  in which the
quasi-parton distribution functions were considered. This scheme has a wide
range of applications Hua:2020gnw ; Gao:2021dbh ; LatticeParton:2022zqc ;
Chou:2022drv ; Hua:2022wop ; Gao:2022iex ; LatticeParton:2022xsd ; Su:2022fiu
; Ji:2022ezo ; Gao:2022ytj ; Gao:2022uhg ; Ji:2022thb ; Zhang:2023tnc ;
Holligan:2023rex ; Zhang:2023bxs ; Gao:2023lny  since its proposal. In this
work, the hybrid scheme for the quasi-DA of a light baryon is developed
utilizing the perturbative coordinate-space quasi-DA calculations, in which
divergences appearing in both long-distance and short-distance spatial
separations can be eliminated properly. By taking a ratio with respect to the
zero-momentum matrix element at short distances, one can eliminate the UV
divergences in lattice matrix elements, and part of the discretization effects
are also expected to be canceled. Correspondingly, the short distance UV
logarithms in the $\overline{\rm MS}$ scheme can also be eliminated properly.
Through the self-renormalization at large distances, one can eliminate the UV
divergences without introducing additional uncontrollable non-perturbative
effects. Since the baryon DA is a two-dimensional distribution, there are
multi regions involving both short distance and large distance simultaneously,
which will be treated separately in this scheme. As a result, this method
ensures that lattice matrix elements approach the continuum limit in a more
appropriate manner and allows for a realistic determination of LCDA.

The rest of this paper is organized as follows: In Sec. 2, we give a brief
overview of the lightcone distribution amplitudes and quasi-DA of a light
baryon. The detailed calculations of the one-loop spatial correlation are
presented in Appendix A. In Sec. 3, the hybrid renormalization scheme is
developed and the matching kernel is presented. Some detailed results are
collected in Appendix B and C. A summary is provided in the last section.

##  2 Lightcone distribution amplitudes and quasi distribution amplitudes for
a light baryon

In this section, we introduce the requisite notations and conventions required
for subsequent discussions. In particular, we will give the definition of
LCDAs, and quasi-DAs and collect the results for the one-loop matching in the
${\overline{\rm MS}}$ scheme.

###  2.1 LCDAs

We start with the LCDAs, which are defined as the hadron-to-vacuum matrix
elements of non-local operators consisting of quarks and gluon which live on
the light cone. In the case of a light baryon, the three-quark matrix element
can be constructed as Braun:1999te

$\left\langle
0\left|\varepsilon^{ijk}u_{\alpha}^{i^{\prime}}\left(z_{1}\right)U_{i^{\prime}i}\left(z_{1},z_{0}\right)d_{\beta}^{j^{\prime}}\left(z_{2}\right)U_{j^{\prime}j}\left(z_{2},z_{0}\right)s_{\gamma}^{k^{\prime}}\left(z_{3}\right)U_{k^{\prime}k}\left(z_{3},z_{0}\right)\right|\Lambda(P,\lambda)\right\rangle,$
(2)

where $\left.\left.\right|\Lambda(P,\lambda)\right\rangle$ stands for the
$\Lambda$ baryon state with the momentum $P$, $P^{2}=0$ and the helicity
$\lambda$. $\alpha$, $\beta$ and $\gamma$ are Dirac indices. $i^{(\prime)}$,
$j^{(\prime)}$ and $k^{(\prime)}$ denote color charges. In this paper, two
light-cone unit vectors are defined as $n^{\mu}=(1,0,0,-1)/\sqrt{2}$ and
$\bar{n}^{\mu}=(1,0,0,1)/\sqrt{2}$. The momentum of the baryon is along the
$\bar{n}$ direction, $P^{\mu}=P^{+}\bar{n}^{\mu}=(P^{z},0,0,P^{z})$. The
coordinates are set in the $n$ direction, $z_{i}^{\mu}=z_{i}n^{\mu}$. The
Wilson lines $U(x,y)$

$U(x,y)=\mathcal{P}\exp\left[ig\int_{0}^{1}\mathrm{~{}d}t(x-y)_{\mu}A^{\mu}(tx+(1-t)y)\right]$
(3)

are inserted to preserve the gauge invariance. In the definition, $z_{0}$ can
be chosen freely due to the gauge invariance, and in the following we will use
$z_{0}=0$ for simplicity. Also for brevity Wilson lines, color indexes, and
helicity will not be written out explicitly below.

Based on Lorentz invariance, and the spin and parity requirement, the matrix
element can be decomposed in terms of three functions, $V(z_{i}P\cdot n)$,
$A(z_{i}P\cdot n)$, and $T(z_{i}P\cdot n)$ to the leading twist

$\displaystyle\left\langle
0\left|u_{\alpha}\left(z_{1}\right)d_{\beta}\left(z_{2}\right)s_{\gamma}\left(z_{3}\right)\right|\Lambda(P)\right\rangle$
(4)
$\displaystyle={f_{N}}\left\\{(P\\!\\!\\!\\!/C)_{\alpha\beta}\left(\gamma_{5}u_{\Lambda}\right)_{\gamma}V\left(z_{i}P\cdot
n\right)+\left(P\\!\\!\\!\\!/\gamma_{5}C\right)_{\alpha\beta}(u_{\Lambda})_{\gamma}A\left(z_{i}P\cdot
n\right)+\left(i\sigma_{\mu\nu}P^{\nu}C\right)_{\alpha\beta}\left(\gamma_{\mu}\gamma_{5}u_{\Lambda}\right)_{\gamma}T\left(z_{i}P\cdot
n\right)\right\\},$

where C signifies the charge conjugation. $u_{\Lambda}$ stands for the
$\Lambda$ baryon spinor. Equivalently, the three leading twist functions can
be projected by inserting a specific gamma matrix $\Gamma$ into the $u$ and
$d$ quark fields. In the following discussion, we will take $A(z_{i}P\cdot n)$
as an example while the other matrix elements can be similarly analyzed. Then
we have

$\displaystyle M_{L}(z_{1},z_{2},z_{3},P^{+},\mu)=\left\langle
0\left|u^{T}\left(z_{1}\right){\Gamma}d\left(z_{2}\right)s\left(z_{3}\right)\right|\Lambda(P)\right\rangle_{R},$
(5)
$\displaystyle{\Phi_{L}}\left(x_{1},x_{2},\mu\right)f_{\Lambda}(\mu)P^{+}u_{\Lambda}(P)=\int_{-\infty}^{+\infty}\frac{d\,P^{+}z_{1}}{2\pi}\frac{d\,P^{+}z_{2}}{2\pi}e^{ix_{1}P^{+}z_{1}+ix_{2}P^{+}z_{2}}M_{L}(z_{1},z_{2},0,P^{+},\mu),$

where $T$ means transpose and ${\Gamma}=C\gamma_{5}\not{n}$. $R$ stands for
renormalization. $x_{i}$s label the longitudinal momentum fractions carried by
the three quarks and $0\leq x_{i}\leq 1$. The $\mu$ denotes the
renormalization scale which will be converted to the factorization scale when
the factorization of quasi-DA is established. $f_{\Lambda}(\mu)$ is the
$\Lambda$ baryon decay constant defined as follows

$f_{\Lambda}(\mu)P^{+}u_{\Lambda}(P)=M_{L}(0,0,0,P^{+},\mu).$ (6)

It should be noted that we have defined the LCDA
${\Phi_{L}}\left(x_{1},x_{2},\mu\right)$ by separating the baryon decay
constant $f_{\Lambda}(\mu)$, which has a different convention with the recent
LQCD calculation RQCD:2019hps . Note that $f_{\Lambda}(\mu)$ depends on the
renormalization scale $\mu$ since the local operator here is not a conserved
current. The LCDA ${\Phi_{L}}\left(x_{1},x_{2},\mu\right)$ in Eq. (5) is
dimensionless and normalized.

###  2.2 Quasi-DAs

Figure 1: One loop corrections for the equal-time matrix element of the
$\Lambda$ baryon.

We consider a spatial correlator

$M(z_{1},z_{2},z_{3},P^{z},\mu)=\left\langle
0\left|u^{T}\left(z_{1}\right)\widetilde{\Gamma}d\left(z_{2}\right)s\left(z_{3}\right)\right|\Lambda(P)\right\rangle_{R}$
(7)

to define the quasi-DA

$\tilde{\Phi}\left(x_{1},x_{2},P^{z},\mu\right)\tilde{f}_{\Lambda}(\mu)P^{z}u_{\Lambda}(P)=(P^{z})^{2}\int_{-\infty}^{+\infty}\frac{d\,z_{1}}{2\pi}\frac{d\,z_{2}}{2\pi}e^{ix_{1}P^{z}z_{1}+ix_{2}P^{z}z_{2}}M(z_{1},z_{2},0,P^{z},\mu),$
(8)

where $\tilde{\Gamma}=C\gamma_{5}n\\!\\!\\!/_{z}$. For the quasi-DAs, the
coordinates are set as $z_{i}^{\mu}=z_{i}n_{z}^{\mu}$, where
$n_{z}^{\mu}=(0,0,0,1)$. The baryon decay constant $\tilde{f}_{\Lambda}(\mu)$
here can be similarly defined as Eq. (5):

$\displaystyle\tilde{f}_{\Lambda}(\mu)=\frac{M(0,0,0,P^{z},\mu)}{P^{z}u_{\Lambda}(P)}.$
(9)

To obtain the quasi-DAs from the lattice simulations, one needs to consider
their short and long-distance properties separately. From this viewpoint, it
is advantageous to work in the coordinate space. To explicitly demonstrate the
ultraviolet and infrared structures, we perform a perturbative calculation of
quasi-DAs by sandwiching the operators between the vacuum state $\left\langle
0|\right.$ and the lowest-order Fock state $\left.|uds\right\rangle$:

${M}_{p}(z_{1},z_{2},z_{3}=0,P^{z},\mu)=\left\langle
0\left|u^{T}\left(z_{1}\right)\widetilde{\Gamma}d\left(z_{2}\right)s\left(0\right)\right|u(x_{1}P)d(x_{2}P)s(x_{3}P)\right\rangle_{R},$
(10)

where the lower index “$p$” in $M_{p}$ denotes the perturbative calculation.
The next-to-leading order Feynman diagrams are shown in Fig. 1, and the
calculation details are collected in Appendix A. The final results for the
spatial correlator up to one loop are given as

$\displaystyle{M}_{p}(z_{1},z_{2},0,P^{z},\mu)=\left\\{1+\frac{\alpha_{s}C_{F}}{\pi}\left(\frac{1}{2}L_{1}^{\text{UV}}+\frac{1}{2}L_{2}^{\text{UV}}+\frac{1}{2}L_{12}^{\text{UV}}+\frac{3}{2}\right)\right\\}{M}_{0}\left(z_{1},z_{2},0,P^{z},\mu\right)$
$\displaystyle-\frac{\alpha_{s}C_{F}}{8\pi}\int_{0}^{1}d\eta_{1}\int_{0}^{1-\eta_{1}}d\eta_{2}$
$\displaystyle\times\left\\{\left(L_{1}^{\text{IR}}-1+\frac{1}{\epsilon_{\mathrm{IR}}}\right){M}_{0}\left(\left(1-\eta_{1}\right)z_{1},z_{2},\eta_{2}z_{1},P^{z},\mu\right)+\left(L_{2}^{\text{IR}}-1+\frac{1}{\epsilon_{\mathrm{IR}}}\right){M}_{0}\left(z_{1},\left(1-\eta_{1}\right)z_{2},\eta_{2}z_{2},P^{z},\mu\right)\right.$
$\displaystyle\left.+2\left(L_{12}^{\text{IR}}-3+\frac{1}{\epsilon_{\mathrm{IR}}}\right){M}_{0}\left(\left(1-\eta_{1}\right)z_{1}+\eta_{1}z_{2},\left(1-\eta_{2}\right)z_{2}+\eta_{2}z_{1},0,P^{z},\mu\right)\right\\}$
$\displaystyle-\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta\times\left\\{{M}_{0}\left((1-\eta)z_{1}+\eta
z_{2},z_{2},0,P^{z},\mu\right)\left\\{\left(L_{12}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}\right.$
$\displaystyle+{M}_{0}\left(z_{1},(1-\eta)z_{2}+\eta
z_{1},0,P^{z},\mu\right)\left\\{\left(L_{12}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}$
$\displaystyle+{M}_{0}\left((1-\eta)z_{1},z_{2},0,P^{z},\mu\right)\left\\{\left(L_{1}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}$
$\displaystyle+{M}_{0}\left(z_{1},(1-\eta)z_{2},0,P^{z},\mu\right)\left\\{\left(L_{2}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}$
$\displaystyle-{M}_{0}\left(z_{1},z_{2},\eta
z_{1},P^{z},\mu\right)\left\\{\left(L_{1}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}$
$\displaystyle\left.-{M}_{0}\left(z_{1},z_{2},\eta
z_{2},P^{z},\mu\right)\left\\{\left(L_{2}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}\right\\},$
(11)

where $M_{0}$ stands for tree-level matrix element:

${M}_{0}(z_{1},z_{2},0,P^{z},\mu)=\sqrt{2}P^{z}e^{ix_{1}P^{z}z_{1}+ix_{2}P^{z}z_{2}}u_{s}(x_{3}P),$
(12)

where $x_{3}=1-x_{1}-x_{2}$, and
$\left[u_{u}\left(x_{1}P\right)\right]^{T}\tilde{\Gamma}u_{d}\left(x_{2}P\right)=\frac{1}{2}\operatorname{tr}\left[P\\!\\!\\!/C\gamma^{5}\tilde{\Gamma}\right]$
is employed. $u_{u/d/s}(P)$ denotes the spinor of u, d, or s quark with
momentum $P$. The plus function is defined as

$\displaystyle\displaystyle\int_{0}^{1}du\left[G(u)\right]_{+}F(u)=\displaystyle\int_{0}^{1}duG(u)[F(u)-F(0)],$
(13)

and some abbreviations are used in the above:

$\displaystyle L_{1}^{\text{IR,
UV}}=\ln\left(\displaystyle\frac{1}{4}\mu_{\text{IR,
UV}}^{2}z_{1}^{2}e^{2\gamma_{E}}\right),\;\;\;L_{2}^{\text{IR,UV}}=\ln\left(\displaystyle\frac{1}{4}\mu_{\text{IR,UV}}^{2}z_{2}^{2}e^{2\gamma_{E}}\right),$
(14) $\displaystyle
L_{12}^{\text{IR,UV}}=\ln\left(\displaystyle\frac{1}{4}\mu^{2}_{\text{IR,UV}}(z_{1}-z_{2})^{2}e^{2\gamma_{E}}\right).$
(15)

We have checked that these results are consistent with the calculation in the
momentum space Deng:2023csv . Moreover, one can see the UV and IR behaviors
clearly in the coordinate space, which is convenient for the renormalization
scheme to be established below.

Furthermore, one can obtain the zero-momentum matrix element in the coordinate
space by letting $P^{z}=0$ and performing the normalization with the local
matrix element

$\displaystyle\hat{M}_{p}\left(z_{1},z_{2},z_{3}=0,P^{z}=0,\mu\right)=\frac{{M}_{p}(z_{1},z_{2},0,0,\mu)}{{M}_{p}(0,0,0,0,\mu)}=1+\frac{\alpha_{s}C_{F}}{2\pi}\left[\frac{7}{8}L_{1}^{\text{UV}}+\frac{7}{8}L_{2}^{\text{UV}}+\frac{3}{4}L_{12}^{\text{UV}}+4\right],$
(16)

where the local matrix element
$\displaystyle{M}_{p}(0,0,0,0,\mu)=\left(1-\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{\epsilon_{\rm
IR}}\right){M}_{0}\left(0,0,0,0,\mu\right)$, see Eq. (80). The perturbative
zero momentum matrix element in Eq. (16) will be used in the hybrid
renormalization method.

###  2.3 Matching of quasi-DAs in the ${\overline{\rm MS}}$ scheme

In the large $P^{z}$ limit, the quasi and lightcone distribution amplitudes
can be related through the QCD factorization. After separating the hard and
collinear contributions, one can factorize the quasi-DAs in terms of the LCDAs
and a hard kernel which can be perturbatively calculated. In the
$\overline{\text{MS}}$ scheme, the one-loop hard kernel has been obtained
Deng:2023csv

$\displaystyle\mathcal{C}_{\rm\overline{MS}}\left(x_{1},x_{2},y_{1},y_{2},P^{z},\mu\right)$
$\displaystyle=\delta\left(x_{1}-y_{1}\right)\delta\left(x_{2}-y_{2}\right)+\frac{\alpha_{s}C_{F}}{8\pi}\times\left[C_{2}\left(x_{1},x_{2},y_{1},y_{2},P^{z},\mu\right)\delta\left(x_{2}-y_{2}\right)\right.$
(17)
$\displaystyle\left.+C_{3}\left(x_{1},x_{2},y_{1},y_{2},P^{z},\mu\right)\delta\left(x_{3}-y_{3}\right)+\left\\{x_{1}\leftrightarrow
x_{2},y_{1}\leftrightarrow y_{2}\right\\}\right]_{\oplus},$

where $\oplus$ denotes a double plus function defined as

${\left[g\left(x_{1},x_{2},y_{1},y_{2}\right)\right]_{\oplus}=}g\left(x_{1},x_{2},y_{1},y_{2}\right)-\delta\left(x_{1}-y_{1}\right)\delta\left(x_{2}-y_{2}\right)\int
dx_{1}^{\prime}dx_{2}^{\prime}g\left(x_{1}^{\prime},x_{2}^{\prime},y_{1},y_{2}\right).$
(18)

$x_{3}=1-x_{1}-x_{2}$, $y_{3}=1-y_{1}-y_{2}$, and

$\displaystyle C_{2}\left(x_{1},x_{2},y_{1},y_{2},P^{z},\mu\right)=$ (19)
$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\frac{\left(x_{1}+y_{1}\right)\left(x_{3}+y_{3}\right)\ln\frac{y_{1}-x_{1}}{-x_{1}}}{y_{1}\left(y_{1}-x_{1}\right)y_{3}}-\frac{x_{3}\left(x_{1}+y_{1}+2y_{3}\right)\ln\frac{x_{3}}{-x_{1}}}{\left(y_{1}-x_{1}\right)y_{3}\left(y_{1}+y_{3}\right)},x_{1}<0\\\
\displaystyle\frac{\left(x_{1}-3y_{1}-2y_{3}\right)x_{1}}{y_{1}\left(x_{3}-y_{3}\right)\left(y_{1}+y_{3}\right)}-\frac{\left[\left(x_{3}-y_{3}\right)^{2}-2x_{3}y_{1}\right]\ln\frac{x_{3}-y_{3}}{x_{3}}}{y_{1}\left(x_{3}-y_{3}\right)y_{3}}+\frac{2x_{1}\ln\frac{4x_{1}\left(x_{3}-y_{3}\right)P_{z}^{2}}{\mu^{2}}}{y_{1}\left(x_{3}-y_{3}\right)}+\frac{x_{1}\ln\frac{4x_{1}x_{3}P_{z}^{2}}{\mu^{2}}}{y_{1}\left(y_{1}+y_{3}\right)},0<x_{1}<y_{1}\\\
\displaystyle\frac{\left(x_{3}-2y_{1}-3y_{3}\right)x_{3}}{y_{3}\left(x_{1}-y_{1}\right)\left(y_{1}+y_{3}\right)}-\frac{\left[\left(x_{1}-y_{1}\right)^{2}-2x_{1}y_{3}\right]\ln\frac{x_{1}-y_{1}}{x_{1}}}{\left(x_{1}-y_{1}\right)y_{1}y_{3}}+\frac{2x_{3}\ln\frac{4x_{3}\left(x_{1}-y_{1}\right)P_{z}^{2}}{\mu^{2}}}{\left(x_{1}-y_{1}\right)y_{3}}+\frac{x_{3}\ln\frac{4x_{1}x_{3}P_{z}^{2}}{\mu^{2}}}{y_{3}\left(y_{1}+y_{3}\right)},y_{1}<x_{1}<y_{1}+y_{3}\\\
\displaystyle\frac{\left(x_{1}+y_{1}\right)\left(x_{3}+y_{3}\right)\ln\frac{y_{3}-x_{3}}{-x_{3}}}{y_{1}y_{3}\left(y_{3}-x_{3}\right)}-\frac{x_{1}\left(x_{3}+2y_{1}+y_{3}\right)\ln\frac{x_{1}}{-x_{3}}}{y_{1}\left(y_{3}-x_{3}\right)\left(y_{1}+y_{3}\right)},x_{1}>y_{1}+y_{3}.\end{array}\right.$
(24) $\displaystyle C_{3}\left(x_{1},x_{2},y_{1},y_{2},P^{z},\mu\right)=$ (25)
$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\frac{\left(x_{1}x_{2}+y_{1}y_{2}\right)\ln\frac{x_{2}-y_{2}}{x_{2}}}{y_{1}\left(x_{2}-y_{2}\right)y_{2}}-\frac{x_{1}\left(x_{2}+y_{1}\right)\ln\frac{-x_{1}}{x_{2}}}{y_{1}\left(x_{2}-y_{2}\right)\left(y_{1}+y_{2}\right)},x_{1}<0\\\
\displaystyle\frac{1}{x_{1}-y_{1}}+\frac{2x_{1}+x_{2}}{y_{1}\left(y_{1}+y_{2}\right)}+\frac{\left[\left(x_{1}+y_{2}\right)y_{1}-x_{1}^{2}\right]\ln\frac{x_{2}-y_{2}}{x_{2}}}{y_{1}\left(x_{2}-y_{2}\right)y_{2}}+\frac{x_{1}\ln\frac{4x_{1}\left(x_{2}-y_{2}\right)P_{z}^{2}}{\mu^{2}}}{y_{1}\left(x_{2}-y_{2}\right)}+\frac{x_{1}\ln\frac{4x_{1}x_{2}P_{z}^{2}}{\mu^{2}}}{y_{1}\left(y_{1}+y_{2}\right)},0<x_{1}<y_{1}\\\
\displaystyle\frac{1}{x_{2}-y_{2}}+\frac{x_{1}+2x_{2}}{y_{2}\left(y_{1}+y_{2}\right)}+\frac{\left[\left(x_{2}+y_{1}\right)y_{2}-x_{2}^{2}\right]\ln\frac{x_{1}-y_{1}}{x_{1}}}{\left(x_{1}-y_{1}\right)y_{1}y_{2}}+\frac{x_{2}\ln\frac{4x_{2}\left(x_{1}-y_{1}\right)P_{z}^{2}}{\mu^{2}}}{\left(x_{1}-y_{1}\right)y_{2}}+\frac{x_{2}\ln\frac{4x_{1}x_{2}P_{z}^{2}}{\mu^{2}}}{y_{2}\left(y_{1}+y_{2}\right)},y_{1}<x_{1}<y_{1}+y_{2}\\\
\displaystyle\frac{\left(x_{1}x_{2}+y_{1}y_{2}\right)\ln\frac{x_{1}-y_{1}}{x_{1}}}{y_{1}\left(x_{1}-y_{1}\right)y_{2}}-\frac{x_{2}\left(x_{1}+y_{2}\right)\ln\frac{-x_{2}}{x_{1}}}{y_{2}\left(x_{1}-y_{1}\right)\left(y_{1}+y_{2}\right)},x_{1}>y_{1}+y_{2}.\end{array}\right.$
(30)

It should be noticed that the above matching of quasi-DAs in the
${\overline{\rm MS}}$ scheme is problematic which can be understood as
follows. As shown in Eq. (18), the double plus function contains an integral
over the momentum fractions $x_{1}^{\prime}$ and $x_{2}^{\prime}$ which arises
from the so-called virtual corrections to DAs Deng:2023csv . In the
$x_{1}^{\prime},x_{2}^{\prime}\to\infty$ limit the integrals have the
asymptotic form $\int dx^{\prime}/x^{\prime}$ and are then divergent. This
divergence corresponds to the logarithmic UV divergence and should be
renormalized. However, in deriving the above matching kernel, one has only
calculated the so-called real diagrams to quasi-DAs and assumed that virtual
contributions are the integration of real diagrams and thus combining them
gives the double plus function. Keeping the remnant UV divergence leads to an
improper definition of double plus function.

In Ref. Izubuchi:2018srq , a properly defined plus function is introduced for
quasi parton distribution functions through the subtraction at the infinite
momentum, while in Ref. Deng:2023csv , the RI/MOM renormalization scheme has
been employed to subtract the asymptotic contributions for the quasi-DAs so
that the plus function becomes well-defined. However it has been shown that
implementing this scheme in the analysis of PDFs introduces unavoidably
additional nonperturbative effects LatticePartonCollaborationLPC:2021xdx  and
thus the final results contain uncontrollable systematic uncertainties.

##  3 Hybrid renormalization scheme for baryon quasi-DA

In this section, we will present the hybrid scheme to address the problems in
the matching as stated in the previous section. The aim of this scheme is to
remove the UV divergences in quasi-DAs without introducing nonperturbative
effects, and meanwhile, this scheme can be properly implemented on the
lattice.

Conceptually, the UV divergences correspond to the small-distance behavior,
while nonperturbative contributions are due to the long-distance behavior. So
from the viewpoint of separating these contributions, it is advantageous to
work in the coordinate space. In coordinate space, the results of quasi-DAs
and the hard kernel can be obtained by making a Fourier transformation of the
results in momentum space. We have also directly calculated the perturbative
contributions to the quasi-DAs in Appendix A, and checked the consistency.

At short distances, e.g. $z_{1}\rightarrow 0$, $z_{2}\rightarrow 0$ or
$|z_{1}-z_{2}|\rightarrow 0$ in Eq. (2.2) with $z_{3}=0$, perturbative quasi-
DAs behave as the logarithmic terms like $\ln(z_{1}^{2})$, $\ln(z_{2}^{2})$
and $\ln((z_{1}-z_{2})^{2})$. As stated in the above section, the $\ln
z_{i}^{2}$ terms in the real diagram of the perturbative matrix elements
prevent them from performing $z\to 0$. This is also reflected in our virtual
contribution definition in which corresponding logarithm UV divergence had
been kept. Subsequently, the double plus function is also not well-defined. By
dividing the perturbative matrix element by a proper zero momentum matrix
element, which is called the ratio scheme Orginos:2017kos ; Radyushkin:2017cyf
; Radyushkin:2017lvu , all these issues arising in the perturbative matrix
elements calculation can be resolved. By construction, in the ratio scheme,
the $\ln z_{i}^{2}$ terms in the numerator and the denominator are equal, then
will be canceled. The corresponding virtual contributions will no longer be
divergent. Then the double plus function can be guaranteed to be well-defined.
Since now the perturbative matrix element is fine for $z\to 0$, the
corresponding lattice matrix elements can also be constructed and used for
matching. It should be noted that the ratio scheme is only available in the
perturbative region. If one performs it over the long-distance region, the IR
structure may be changed. These require that only zero momentum matrix
elements in the short-distance region can be chosen as denominators and the
renormalization for the matrix elements in the long-distance region needs
additional consideration.

At large distances, e.g. $z_{1}\sim 1/\Lambda_{\rm QCD}$, $z_{2}\sim
1/\Lambda_{\rm QCD}$ or $|z_{1}-z_{2}|\sim 1/\Lambda_{\rm QCD}$ in Eq. (7), a
proper renormalization scheme should not introduce extra non-perturbative
effects. In this work, we use the self-renormalization scheme as advocated in
Ref. LatticePartonCollaborationLPC:2021xdx . In the self-renormalization
scheme, a renormalization factor $Z_{R}$, which includes all typical
divergences and discretion errors, is defined. One can use this
renormalization factor to convert the lattice matrix elements to continuous
matrix elements without bringing in any non-perturbative effects. Besides the
UV divergence inherited by dimensional regularization, the linear divergence
from the lattice simulation is eliminated as well. By definition, the
calculation of the renormalization factors requires the input of the UV
divergences and some parameters. Specifically, the UV divergences are fitted
from zero momentum lattice matrix elements at small lattice spacings, and
these parameters can be obtained through matching the renormalized lattice
matrix element to the continuum perturbative matrix element in the
perturbative region. As a result, the long-distance regions which merely
involve non-perturbative scales in the renormalization can also be handled.

A subtlety in applying the hybrid renormalization scheme to quasi-DAs of a
light baryon is that there are multi regions involving both perturbative and
non-perturbative scales, e.g. ($z_{1}\rightarrow 0$ and $z_{2}\sim
1/\Lambda_{\rm QCD}$), ($z_{2}\rightarrow 0$ and $z_{1}\sim 1/\Lambda_{\rm
QCD}$), or ($|z_{1}-z_{2}|\rightarrow 0$, $z_{1}\sim 1/\Lambda_{\rm QCD}$ and
$z_{2}\sim 1/\Lambda_{\rm QCD}$), see the blue bands in Fig. 4 except for the
one around $z_{1}\sim-z_{2}$. However, if the logarithmic divergences
regarding different scales ($z_{1}$, $z_{2}$ and $|z_{1}-z_{2}|$) can be
factorized separately, the UV divergences regarding those scales can be
multiplicatively renormalized separately. If so, the ratio scheme can be taken
for the perturbative scale and the self renormalization can be performed for
the non-perturbative scale. From the one-loop result in Eq. (16), one can see
that the logarithmic divergences indeed factorize, which allows to perform the
hybrid renormalization. An all-order proof is left for future studies.

In this section, we will give the normalization and self renormalization in
the beginning, which provide the building blocks for hybrid renormalization.
Then we provide a method for hybrid renormalization. The matching kernel in
the hybrid scheme is presented at the end.

###  3.1 Normalization

One starts from a quasi-DA matrix element on lattice
$M\left(z_{1},z_{2},z_{3}=0,P^{z},a\right)$, which is defined in Eq. (7). On
the lattice, the matrix element is regularized with the lattice spacing $a$.
To ensure the normalization, one has

$\hat{M}\left(z_{1},z_{2},0,P^{z},a\right)=\frac{M\left(z_{1},z_{2},0,P^{z},a\right)}{M\left(0,0,0,P^{z},a\right)},$
(31)

where part of the discretization effects are eliminated in the lattice matrix
elements through the normalization.

###  3.2 Self renormalization

In this subsection, we discuss the extraction of the renormalization factor
$Z_{R}(z_{1},z_{2},a,\mu)$ from the zero momentum matrix element
$\hat{M}\left(z_{1},z_{2},0,0,a\right)$ (Eq. (31) if $P^{z}=0$). This
renormalization factor will be applied to the large distances in hybrid
renormalization.

According to Ref. LatticePartonCollaborationLPC:2021xdx , the renormalization
factor is an asymptotic expansion with respect to $a$, with both power
dependence and logarithmic dependence

$Z_{R}(z_{1},z_{2},a,\mu)=\exp\Big{[}\left(\frac{k}{a\ln[a\Lambda_{\rm
QCD}]}-m_{0}\right)\tilde{z}+\frac{\gamma_{0}}{b_{0}}\ln\bigg{[}\frac{\ln[1/(a\Lambda_{\rm
QCD})]}{\ln[\mu/\Lambda_{\rm\overline{MS}}]}\bigg{]}+\ln\left[1+\frac{d}{\ln(a\Lambda_{\rm
QCD})}\right]+f(z_{1},z_{2})a\Big{]}\ ,$ (32)

where $\displaystyle\left(\frac{k}{a\ln[a\Lambda_{\rm
QCD}]}-m_{0}\right)\tilde{z}$ is the linear divergence Chen:2016fxx ;
Ji:2017oey ; Ishikawa:2017faj ; Green:2017xeu ; Ji:2020brr  and the mass
renormalization parameter Ji:1995tm ; Beneke:1998ui ; Bauer:2011ws ;
Bali:2013pla ; Zhang:2023bxs ,
$\displaystyle\frac{\gamma_{0}}{b_{0}}\ln\bigg{[}\frac{\ln[1/(a\Lambda_{\rm
QCD})]}{\ln[\mu/\Lambda_{\rm\overline{MS}}]}\bigg{]}+\ln\left[1+\frac{d}{\ln(a\Lambda_{\rm
QCD})}\right]$ contains the log divergence. $f(z_{1},z_{2})a$ (or
$f(z_{1},z_{2})a^{2}$) is the discretization effect. $\tilde{z}$ is the
effective length for the linear divergence, which is defined as follows

$\displaystyle\tilde{z}=\left\\{\begin{array}[]{ll}|z_{1}-z_{2}|,&\quad
z_{1}z_{2}<0\\\ {\rm max}\left(|z_{1}|,|z_{2}|\right),&\quad z_{1}z_{2}\geq
0.\end{array}\right.$ (35)

Figure 2: The schematic diagram of the Wilson links for $z_{1}z_{2}<0$.
$z_{1}$ and $z_{2}$ denote the locations of the heads of the Wilson links. The
tails of the Wilson links are located at the origin, where the color indices
are contracted with the Levi-Civita symbol $\epsilon$ to guarantee the gauge
invariance.

For $z_{1}z_{2}<0$, the effective length is the total length of the Wilson
link, see Fig. 2. For $z_{1}z_{2}>0$, there is an overlap region between the
two Wilson links. However, one can simply show (e.g. $z_{1}>z_{2}>0$)

$\displaystyle
U(z_{1},0)_{i^{\prime}i}U(z_{2},0)_{j^{\prime}j}\epsilon_{ijk}=U(z_{1},z_{2})_{i^{\prime}i^{\prime\prime}}U(z_{2},0)_{i^{\prime\prime}i}U(z_{2},0)_{j^{\prime}j}\epsilon_{ijk}$
$\displaystyle=U(z_{1},z_{2})_{i^{\prime}i^{\prime\prime}}U(z_{2},0)_{i^{\prime\prime}i}U(z_{2},0)_{j^{\prime}j}\delta_{kk^{\prime\prime}}\epsilon_{ijk^{\prime\prime}}$
$\displaystyle=U(z_{1},z_{2})_{i^{\prime}i^{\prime\prime}}U(z_{2},0)_{i^{\prime\prime}i}U(z_{2},0)_{j^{\prime}j}U^{{\dagger}}(z_{2},0)_{kk^{\prime}}U(z_{2},0)_{k^{\prime}k^{\prime\prime}}\epsilon_{ijk^{\prime\prime}}$
$\displaystyle=U(z_{1},z_{2})_{i^{\prime}i^{\prime\prime}}U^{{\dagger}}(z_{2},0)_{kk^{\prime}}{\rm
det}[U(z_{2},0)]\epsilon_{i^{\prime\prime}j^{\prime}k^{\prime}}=U(z_{1},z_{2})_{i^{\prime}i^{\prime\prime}}U^{{\dagger}}(z_{2},0)_{kk^{\prime}}\epsilon_{i^{\prime\prime}j^{\prime}k^{\prime}},$
(36)

where the unitary and special properties of the Wilson links as SU(3) group
elements are used in the second line and last line respectively. A schematic
diagram for the above relation is shown in Fig. 3. So the effective length is
$z_{1}$ for $z_{1}>z_{2}>0$. For a general case of $z_{1}z_{2}\geq 0$, the
effective length should be the maximum value between $|z_{1}|$ and $|z_{2}|$.

Figure 3: The schematic diagram of the Wilson links for $z_{1}z_{2}>0$. The
notations are the same as Fig. 2. Considering the color contraction
$\epsilon$, the two overlapped Wilson links in the same direction are
equivalent to a Wilson link in the opposite direction.

The leading order QCD beta function $\displaystyle
b_{0}=\frac{11C_{A}-2n_{f}}{6\pi}$ which satisfies
$\displaystyle\frac{d\alpha_{s}}{d\ln[\mu]}=-b_{0}\alpha_{s}^{2}$. The
perturbative zero momentum matrix element in $\overline{\rm MS}$ scheme Eq.
(16) satisfies the following renormalization group equation

$\frac{d\ln[\hat{M}_{p}\left(z_{1},z_{2},0,0,\mu\right)]}{d\ln[\mu]}=\gamma,$
(37)

where $\gamma=\gamma_{0}\alpha_{s}+...$ The leading anomalous dimension
$\displaystyle\gamma_{0}=\frac{C_{F}}{2\pi}\left(5-\frac{7}{4}\delta_{z_{1},0}-\frac{7}{4}\delta_{z_{2},0}-\frac{3}{2}\delta_{z_{1}-z_{2},0}\right)$
is scheme independent which can be applied to the renormalization factor Eq.
(32) of lattice matrix elements. It involves the quark-link interaction as
well as the evolution effect of the local operator. There are subtraction
terms since the UV fluctuation is frozen on lattice for a distance to be zero.
$\Lambda_{\rm\overline{MS}}$ is the RG invariant scale for the LO running
coupling, which is 0.142 GeV for $n_{f}=3$, 0.119 GeV for $n_{f}=4$ and 0.087
GeV for $n_{f}=5$, determined based on the method in Karbstein:2018mzo .

The parameters $k$, $\Lambda_{\rm QCD}$, $f(z_{1},z_{2})$, $m_{0}$ and $d$ are
extracted through the fit. The fit procedure is
LatticePartonCollaborationLPC:2021xdx :
1) fit the $a$ dependence in $\hat{M}\left(z_{1},z_{2},0,0,a\right)$ to
extract the global parameters $k$ and $\Lambda_{\rm QCD}$ as well as the
discretization effect $f(z_{1},z_{2})$:

$\displaystyle\hat{M}\left(z_{1},z_{2},0,0,a\right)$
$\displaystyle=\exp\Big{[}\frac{k}{a\ln[a\Lambda_{\rm
QCD}]}\tilde{z}+g(z_{1},z_{2},d)+\frac{\gamma_{0}}{b_{0}}\ln\bigg{[}\frac{\ln[1/(a\Lambda_{\rm
QCD})]}{\ln[\mu/\Lambda_{\rm\overline{MS}}]}\bigg{]}+\ln\left[1+\frac{d}{\ln(a\Lambda_{\rm
QCD})}\right]+f(z_{1},z_{2})a\Big{]}\ ,$ (38)

where $g(z_{1},z_{2},d)$ contains the non-perturbative intrinsic $z_{1},z_{2}$
dependences, which is also extracted through the fit. $g(z_{1},z_{2},d)$
depends on the choice of the global parameter $d$ during the fit;
2) extract $m_{0}$ and $d$ through requiring the renormalized matrix element
to be equal to the perturbative matrix element at short distances
($a<z_{1},z_{2}\ll 1/\Lambda_{\rm QCD}$)

$\frac{\hat{M}\left(z_{1},z_{2},0,0,a\right)}{Z_{R}(z_{1},z_{2},a,\mu)}=\exp\Big{[}g(z_{1},z_{2},d)+m_{0}\tilde{z}\Big{]}=\hat{M}_{p}\left(z_{1},z_{2},0,0,\mu\right).$
(39)

The perturbative zero momentum matrix element is obtained from Eqs. (2.2) and
(16),

$\displaystyle\hat{M}_{p}\left(z_{1},z_{2},0,0,\mu\right)=1+\frac{\alpha_{s}C_{F}}{2\pi}\left[\frac{7}{8}\ln\left(\frac{z_{1}^{2}\mu^{2}e^{2\gamma_{E}}}{4}\right)+\frac{7}{8}\ln\left(\frac{z_{2}^{2}\mu^{2}e^{2\gamma_{E}}}{4}\right)+\frac{3}{4}\ln\left(\frac{(z_{1}-z_{2})^{2}\mu^{2}e^{2\gamma_{E}}}{4}\right)+4\right],$
(40)

where the IR poles are canceled through the normalization with the local
matrix element. The details of the perturbative calculation are presented in
Appendix A.

Thus one can define the renormalized lattice matrix element in $\overline{\rm
MS}$ scheme for the whole range as the ratio of the normalized lattice matrix
element Eq. (31) to the renormalization factor Eq. (32)

$\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)=\frac{\hat{M}\left(z_{1},z_{2},0,P^{z},a\right)}{Z_{R}(z_{1},z_{2},a,\mu)},$
(41)

where the renormalization factor $Z_{R}$, though extracted from the zero
momentum matrix element, can be applied to the large momentum matrix element
since the renormalization is independent of the external states.

###  3.3 Hybrid renormalization

A hybrid renormalization method is presented in this subsection, based on the
building blocks provided in the previous subsections, such as the normalized
lattice matrix element Eq. (31), the renormalization factor Eq. (32) and the
renormalized lattice matrix element in $\overline{\rm MS}$ scheme Eq. (41).

Figure 4: A schematic diagram of renormalization. First, both large and zero
momentum matrix elements are converted to $\overline{\rm MS}$ scheme. Then,
additional ratios are taken with the $\overline{\rm MS}$ scheme matrix
elements as follows. The large momentum matrix element in the purple (hard)
region is divided by the zero momentum matrix element in the purple region
correspondingly. The large momentum matrix in the French blue (hard-soft)
region is divided by the zero momentum matrix element on the blue-purple
boundary. The large momentum matrix in the white (soft) region is divided by
the zero momentum matrix element at the white-purple intersection point.

In practical calculations, it is difficult to directly perform the matching in
$\overline{\rm MS}$ scheme since there are inconsistencies between the lattice
matrix element and continuum scheme at short distances. The lattice matrix
elements are finite as $z_{1}\rightarrow 0$ or $z_{2}\rightarrow 0$ while the
perturbative matrix elements, with the logarithmic terms $\ln(z_{1}^{2})$,
$\ln(z_{2}^{2})$ and $\ln((z_{1}-z_{2})^{2})$, are divergent as
$z_{1}\rightarrow 0$ or $z_{2}\rightarrow 0$. Moreover, those logarithmic
terms correspond to slowly decaying terms
$\sim\displaystyle\frac{1}{|x^{\prime}_{1}|}$ and
$\displaystyle\frac{1}{|x^{\prime}_{2}|}$ in the matching kernel, which
increase the difficulties of numerical calculations in preserving the
normalization, though not impossible. The hybrid scheme can be viewed as a
modification of $\overline{\rm MS}$ scheme at short distance. Through the
ratio at short distance, part of the discretization effects are canceled in
the lattice matrix elements and the singular log terms are canceled in the
perturbative matrix elements. So the lattice matrix elements become more
consistent with the continuum scheme under the hybrid scheme, where it is
easier to preserve the normalization.

Thus there are several principles in designing the hybrid scheme
LatticePartonCollaborationLPC:2021xdx :
1) Eliminate all the singular logarithmic terms in the perturbative matrix
elements including $\ln(z_{1}^{2})$, $\ln(z_{2}^{2})$ and
$\ln((z_{1}-z_{2})^{2})$ through the ratio;
2) Avoid introducing extra effects that are not perturbatively controllable at
large distances;
3) Keep the renormalized matrix element continuous and the method as simple as
possible.

A hybrid renormalization method that satisfies the above principles is
presented in the following, and a schematic diagram is shown in Fig. 4. For
the sake of convenience, besides $|z_{1}|$, $|z_{2}|$ and $|z_{1}-z_{2}|$, we
also treat $|z_{1}+z_{2}|$ as an argument. The purple region denotes that all
arguments are perturbative scales. The white regions denote that all are non-
perturbative scales. The blue regions denote that only one of the four scales
is perturbative.

  * •

If the scales $z_{1}$, $z_{2}$ and $z_{1}-z_{2}$ are all perturbative, e.g.
($|z_{1}|<z_{s}$ and $|z_{2}|<z_{s}$) or ($|z_{1}|<z_{s}$ and
$z_{s}<|z_{2}|<2z_{s}$) or ($z_{s}<|z_{1}|<2z_{s}$ and $|z_{2}|<z_{s}$), which
corresponds to the purple region in Fig. 4, one introduces the ratio scheme on
the normalized lattice matrix elements (Eq. (31))

$\displaystyle\frac{\hat{M}\left(z_{1},z_{2},0,P^{z},a\right)}{\hat{M}\left(z_{1},z_{2},0,0,a\right)}\left(\theta(z_{s}-|z_{1}|)\theta(z_{s}-|z_{2}|)+\theta(z_{s}-|z_{2}|)\theta(|z_{1}|-z_{s})\theta(2z_{s}-|z_{1}|)\right.$
(42)
$\displaystyle\left.+\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-z_{s})\theta(2z_{s}-|z_{2}|)\right)$
$\displaystyle=\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,0,\mu\right)}\left(\theta(2z_{s}-|z_{1}|)\theta(z_{s}-|z_{2}|)+\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-z_{s})\theta(2z_{s}-|z_{2}|)\right),$

where $z_{s}$ is the hybrid cutoff, which satisfies $a\ll 2z_{s}\ll
1/\Lambda_{\rm QCD}$. So no extra non-perturbative effects are introduced. The
ratio can be written with the renormalized lattice matrix element in
$\overline{\rm MS}$ scheme Eq. (41) since the renormalization factor $Z_{R}$
is independent of momentum $P^{z}$. There are several advantages to taking the
ratio scheme at short distances. First, part of the discretization effects are
canceled since $a\ll z_{s}$. Second, the normalization can be guaranteed in
both renormalized lattice matrix elements and perturbative matching. In the
perturbative matrix element, we will see the cancellation of $\ln(z_{1}^{2})$,
$\ln(z_{2}^{2})$ and $\ln((z_{1}-z_{2})^{2})$ through the ratio, which leads
to a current conserved matching kernel after Fourier transformation.

  * •

If $z_{1}$ is perturbative while $z_{2}$ and $z_{1}-z_{2}$ are not, e.g.
$|z_{1}|<z_{s}$ and $|z_{2}|>2z_{s}$, which is the blue vertical region in
Fig. 4, one needs to introduce the ratio scheme for $z_{1}$ and $\overline{\rm
MS}$ scheme (Eq. (41)) for $z_{2}$,

$\displaystyle\frac{\hat{M}\left(z_{1},z_{2},0,P^{z},a\right)Z_{R}(z_{1},{\rm
sign}(z_{2})2z_{s},a,\mu)}{\hat{M}\left(z_{1},{\rm
sign}(z_{2})2z_{s},0,0,a\right)Z_{R}(z_{1},z_{2},a,\mu)}\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-2z_{s})$
(43) $\displaystyle=\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{1},{\rm
sign}(z_{2})2z_{s},0,0,\mu\right)}\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-2z_{s}),$

where $Z_{R}(z_{1},z_{2},a,\mu)$ is the renormalization factor to convert the
lattice matrix element to $\overline{\rm MS}$ scheme as we have defined in Eq.
(32). The zero momentum matrix element in the denominator
$\hat{M}\left(z_{1},{\rm sign}(z_{2})2z_{s},0,0,a\right)$ will be crucial in
canceling the $\ln(z_{1}^{2})$ in the perturbative matrix element when we
deduce the matching kernel. No extra non-perturbative effect is introduced in
the denominator where the $z_{2}$ dependence is truncated at ${\rm
sign}(z_{2})2z_{s}$.

  * •

If $z_{2}$ is perturbative while $z_{1}$ and $z_{1}-z_{2}$ are not, e.g.
$|z_{1}|>2z_{s}$ and $|z_{2}|<z_{s}$, which is the blue horizontal region in
Fig. 4, one follows the similar strategy,

$\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left({\rm
sign}(z_{1})2z_{s},z_{2},0,0,\mu\right)}\theta(|z_{1}|-2z_{s})\theta(z_{s}-|z_{2}|).$
(44)

  * •

If both $z_{1}$ and $z_{2}$ are non-perturbative while $z_{1}-z_{2}$ is
perturbative, e.g. $|z_{1}|>z_{s}$, $|z_{2}|>z_{s}$ and $|z_{1}-z_{2}|<z_{s}$,
which is the blue diagonal region around $z_{1}\sim z_{2}$ in Fig. 4, one
takes

$\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{1}^{*},z_{2}^{*},0,0,\mu\right)}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}-z_{2}|),$
(45)

where $z_{1}^{*}=z_{s}+(z_{1}-z_{2})\theta(z_{1}-z_{2})$ and
$z_{2}^{*}=z_{s}+(z_{2}-z_{1})\theta(z_{2}-z_{1})$. The coordinate choices
$z_{1}^{*}$ and $z_{2}^{*}$ apply for $z_{1}<0$ and $z_{2}<0$ as well because
($z_{1}\rightarrow-z_{2},z_{2}\rightarrow-z_{1}$) is a symmetry for zero
momentum matrix element. The zero momentum matrix element in the denominator
will be crucial in canceling the $\ln((z_{1}-z_{2})^{2})$ in the perturbative
matrix element.

  * •

For $|z_{1}|>z_{s}$, $|z_{2}|>z_{s}$ and $|z_{1}+z_{2}|<z_{s}$, which is the
blue diagonal region around $z_{1}\sim-z_{2}$ in Fig. 4, one can introduce the
ratio as well

$\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{1}^{**},z_{2}^{**},0,0,\mu\right)}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}+z_{2}|),$
(46)

where $z_{1}^{**}=z_{s}+(z_{1}+z_{2})\theta(z_{1}+z_{2})$ and
$z_{2}^{**}=-z_{s}+(z_{2}+z_{1})\theta(-z_{2}-z_{1})$. This step is for
continuity and simplicity.

  * •

Finally, for $|z_{1}|>z_{s}$, $|z_{2}|>z_{s}$, $|z_{1}-z_{2}|>z_{s}$ and
$|z_{1}+z_{2}|>z_{s}$, we apply the $\overline{\rm MS}$ scheme (Eq. (41)),

$\displaystyle\frac{\hat{M}\left(z_{1},z_{2},0,P^{z},a\right)}{Z_{R}(z_{1},z_{2},a,\mu)\hat{M}_{\overline{\rm
MS}}\left({\rm sign}(z_{1})z_{s},{\rm
sign}(z_{2})2z_{s},0,0,\mu\right)}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(|z_{1}-z_{2}|-z_{s})\theta(|z_{1}+z_{2}|-z_{s})$
(47) $\displaystyle=\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left({\rm sign}(z_{1})z_{s},{\rm
sign}(z_{2})2z_{s},0,0,\mu\right)}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(|z_{1}-z_{2}|-z_{s})\theta(|z_{1}+z_{2}|-z_{s}).$

The isospin and parity symmetries allow different choices for the denominator:

$\hat{M}_{\overline{\rm
MS}}\left(z_{s},2z_{s},0,0,\mu\right)=\hat{M}_{\overline{\rm
MS}}\left(2z_{s},z_{s},0,0,\mu\right)=\hat{M}_{\overline{\rm
MS}}\left(-z_{s},-2z_{s},0,0,\mu\right)=\hat{M}_{\overline{\rm
MS}}\left(-2z_{s},-z_{s},0,0,\mu\right),$ $\hat{M}_{\overline{\rm
MS}}\left(-z_{s},2z_{s},0,0,\mu\right)=\hat{M}_{\overline{\rm
MS}}\left(2z_{s},-z_{s},0,0,\mu\right)=\hat{M}_{\overline{\rm
MS}}\left(z_{s},-2z_{s},0,0,\mu\right)=\hat{M}_{\overline{\rm
MS}}\left(-2z_{s},z_{s},0,0,\mu\right).$

To conclude, the hybrid renormalized matrix element is

$\displaystyle\hat{M}_{H}(z_{1},z_{2},0,P^{z})=\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,0,\mu\right)}\left(\theta(2z_{s}-|z_{1}|)\theta(z_{s}-|z_{2}|)+\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-z_{s})\theta(2z_{s}-|z_{2}|)\right)$
(48) $\displaystyle+\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{1},{\rm
sign}(z_{2})2z_{s},0,0,\mu\right)}\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-2z_{s})+\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left({\rm
sign}(z_{1})2z_{s},z_{2},0,0,\mu\right)}\theta(|z_{1}|-2z_{s})\theta(z_{s}-|z_{2}|)$
$\displaystyle+\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{s}+(z_{1}-z_{2})\theta(z_{1}-z_{2}),z_{s}+(z_{2}-z_{1})\theta(z_{2}-z_{1}),0,0,\mu\right)}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}-z_{2}|)$
$\displaystyle+\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left(z_{s}+(z_{1}+z_{2})\theta(z_{1}+z_{2}),-z_{s}+(z_{2}+z_{1})\theta(-z_{2}-z_{1}),0,0,\mu\right)}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}+z_{2}|)$
$\displaystyle+\frac{\hat{M}_{\overline{\rm
MS}}\left(z_{1},z_{2},0,P^{z},\mu\right)}{\hat{M}_{\overline{\rm
MS}}\left({\rm sign}(z_{1})z_{s},{\rm
sign}(z_{2})2z_{s},0,0,\mu\right)}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(|z_{1}-z_{2}|-z_{s})\theta(|z_{1}+z_{2}|-z_{s}),$

where $\hat{M}_{\overline{\rm MS}}$ is the renormalized lattice matrix element
in $\overline{\rm MS}$ scheme defined in Eq. (41).

###  3.4 Matching in hybrid scheme

The hybrid scheme quasi baryon DA is obtained through the Fourier
transformation

$\tilde{\Phi}_{H}(x_{1},x_{2},P^{z})=(P^{z})^{2}\int_{-\infty}^{+\infty}\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{ix_{1}P^{z}z_{1}+ix_{2}P^{z}z_{2}}\hat{M}_{H}\left(z_{1},z_{2},0,P^{z}\right).$
(49)

The factorization formula is

$\tilde{\Phi}_{H}(x_{1},x_{2},P^{z})=\int
dy_{1}dy_{2}\mathcal{C}_{H}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)\Phi_{L}(y_{1},y_{2},\mu),$
(50)

where $\Phi_{L}(y_{1},y_{2},\mu)$ is the baryon LCDA in Eq. (5).
$\mathcal{C}_{H}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)$ is the matching kernel in
hybrid scheme,

$\mathcal{C}_{H}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)=\mathcal{C}_{\overline{\rm
MS}}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)-\delta\mathcal{C}_{H}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu),$
(51)

where $\mathcal{C}_{\overline{\rm MS}}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)$ is
the $\overline{\rm MS}$ scheme matching kernel in Eq. (17).
$\delta\mathcal{C}_{H}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)$ is the hybrid
counterterm and the NLO result can be obtained through the Fourier
transformation Izubuchi:2018srq

$\displaystyle\delta\mathcal{C}_{H}^{(1)}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)=(P^{z})^{2}\int\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{i(x_{1}-y_{1})P^{z}z_{1}+i(x_{2}-y_{2})P^{z}z_{2}}\delta\hat{M}_{H}^{(1)}\left(z_{1},z_{2},\mu\right).$
(52)

$\delta\hat{M}_{H}^{(1)}\left(z_{1},z_{2},\mu\right)$ is the perturbative
correction in the hybrid scheme at NLO:

$\displaystyle\delta\hat{M}_{H}^{(1)}\left(z_{1},z_{2},\mu\right)=\hat{M}_{p}^{(1)}\left(z_{1},z_{2},0,0,\mu\right)\left(\theta(2z_{s}-|z_{1}|)\theta(z_{s}-|z_{2}|)+\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-z_{s})\theta(2z_{s}-|z_{2}|)\right)$
$\displaystyle+\hat{M}_{p}^{(1)}\left(z_{1},{\rm
sign}(z_{2})2z_{s},0,0,\mu\right)\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-2z_{s})+\hat{M}_{p}^{(1)}\left({\rm
sign}(z_{1})2z_{s},z_{2},0,0,\mu\right)\theta(|z_{1}|-2z_{s})\theta(z_{s}-|z_{2}|)$
$\displaystyle+\hat{M}_{p}^{(1)}\left(z_{s}+(z_{1}-z_{2})\theta(z_{1}-z_{2}),z_{s}+(z_{2}-z_{1})\theta(z_{2}-z_{1}),0,0,\mu\right)\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}-z_{2}|)$
$\displaystyle+\hat{M}_{p}^{(1)}\left(z_{s}+(z_{1}+z_{2})\theta(z_{1}+z_{2}),-z_{s}+(z_{2}+z_{1})\theta(-z_{2}-z_{1}),0,0,\mu\right)\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}+z_{2}|)$
$\displaystyle+\hat{M}_{p}^{(1)}\left({\rm sign}(z_{1})z_{s},{\rm
sign}(z_{2})2z_{s},0,0,\mu\right)\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(|z_{1}-z_{2}|-z_{s})\theta(|z_{1}+z_{2}|-z_{s}).$
(53)

$\hat{M}_{p}^{(1)}\left(z_{1},z_{2},0,0,\mu\right)$ is the NLO result of Eq.
(16),

$\hat{M}_{p}^{(1)}\left(z_{1},z_{2},0,0,\mu\right)=\frac{\alpha_{s}C_{F}}{2\pi}\left[\frac{7}{8}\ln\left(\frac{z_{1}^{2}\mu^{2}e^{2\gamma_{E}}}{4}\right)+\frac{7}{8}\ln\left(\frac{z_{2}^{2}\mu^{2}e^{2\gamma_{E}}}{4}\right)+\frac{3}{4}\ln\left(\frac{(z_{1}-z_{2})^{2}\mu^{2}e^{2\gamma_{E}}}{4}\right)+4\right].$

One performs the Fourier transformation on different regions and obtains the
hybrid counterterm:

$\displaystyle\delta\mathcal{C}_{H}^{(1)}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)=(P^{z})^{2}\frac{\alpha_{s}C_{F}}{2\pi}\Bigg{[}I_{\rm
H}[(x_{1}-y_{1})P^{z},(x_{2}-y_{2})P^{z}]+I_{\rm
HSI}[(x_{1}-y_{1})P^{z},(x_{2}-y_{2})P^{z}]$ $\displaystyle+I_{\rm
HSII}[(x_{1}-y_{1})P^{z},(x_{2}-y_{2})P^{z}]+I_{\rm
HSIII}[(x_{1}-y_{1})P^{z},(x_{2}-y_{2})P^{z}]+I_{\rm
HSIV}[(x_{1}-y_{1})P^{z},(x_{2}-y_{2})P^{z}]$ $\displaystyle+I_{\rm
S}[(x_{1}-y_{1})P^{z},(x_{2}-y_{2})P^{z}]+\delta[(x_{1}-y_{1})P^{z}]\delta[(x_{2}-y_{2})P^{z}]\left(\frac{5}{2}\ln\left(\frac{\mu^{2}e^{2\gamma_{E}}}{4}\right)+4\right)\Bigg{]},$
(54)

where all integrated formulas are collected in Appendix B.

Since the normalization with respect to the local matrix element is performed
and the ratio is taken with respect to the zero momentum matrix element at
short distances, the matching in the hybrid scheme preserves the normalization
requirement. That means, if the light cone DA is normalized $\int
dy_{1}dy_{2}\Phi_{L}(y_{1},y_{2},\mu)=1$, the quasi-DA obtained through
matching in the hybrid scheme is also normalized $\int
dx_{1}dx_{2}\tilde{\Phi}_{H}(x_{1},x_{2},P^{z})=1$. As shown in Appendix C, if
the normalization is preserved, the hybrid matching kernel at one loop is a
pure double plus function defined as Eq. (18),

$\mathcal{C}_{H}^{(1)}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)=\left[\mathcal{C}^{(1)}_{\overline{\rm
MS}}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)-\delta\mathcal{C}^{(1)}_{H}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)\right]_{\oplus}.$
(55)

A necessary condition for the pure double plus function is that there are no
linearly decaying terms $\sim\displaystyle\frac{1}{|x^{\prime}_{1}|}$ and
$\displaystyle\frac{1}{|x^{\prime}_{2}|}$ for $|x^{\prime}_{1}|\gg 1$ and
$|x^{\prime}_{2}|\gg 1$ respectively. One can check those linearly decaying
terms in $\mathcal{C}^{(1)}_{\overline{\rm
MS}}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)$ are explicitly cancelled by the
linearly decaying terms in
$\delta\mathcal{C}^{(1)}_{H}(x_{1},x_{2},y_{1},y_{2},P^{z},\mu)$, which is
equivalent to the cancellation of the log terms $\ln(z_{1}^{2})$,
$\ln(z_{2}^{2})$ and $\ln((z_{1}-z_{2})^{2})$ through the ratio at short
distance in coordinate space. For readers’ convenience, a Mathematica notebook
on the hybrid scheme matching kernel is attached to the source package on
arXiv.

##  4 Summary

To summarize, this paper continues the work on a direct method to extract the
shape distribution of LCDA of a light baryon through the simulation of equal-
time correlation functions under the framework of large momentum effective
theory. To clear the obstacles on the road to renormalize the quasi-DAs, we
have developed a hybrid renormalization scheme. By combining the self-
renormalization at large spatial separations and the ratio scheme at short
spatial separations, the hybrid renormalization scheme removes the UV
divergences without introducing extra nonperturbative effects. The
corresponding equal-time correlation functions have been calculated in the
coordinate space to the next-to-leading order, and the matching kernel between
the quasi-DAs and LCDAs has been derived.

In the calculation, the $\Lambda$ baryon has been taken as an example to
demonstrate the scheme, but such a scheme can be straightforwardly generalized
to other octet and decuplet baryons. This approach offers a practical
methodology for computing LCDA which can be carried out in lattice
calculations, and a preliminary analysis using this approach is now being
conducted by members of Larton Parton collaboration.

###### Acknowledgements.

We would like to thank Jun Zeng, Zhifu Deng, Minhuan Chu, Jun Hua, and Qi-An
Zhang for their insightful comments and invaluable discussions. This work is
supported in part by the Natural Science Foundation of China under Grants
No.12125503, No. U2032102, No. 12061131006, and No. 12335003. Y.S. is
partially supported by the U.S. Department of Energy, Office of Science,
Office of Nuclear Physics, contract no. DE-AC02-06CH11357. Part of the
computations in this paper was run on the Siyuan-1 cluster supported by the
Center for High Performance Computing at Shanghai Jiao Tong University.

##  Appendix A One-loop calculation of the spatial correlator in coordinate
space

In this section, we will present one-loop results for the spatial correlator
in dimensional regularization with $\overline{\text{MS}}$ renormalization. The
results are gauge invariant, and the Feynman gauge will be adopted in the
practical calculation.

As shown in Fig. 1, there are twelve distinct diagrams to calculate which can
be divided into three categories: quark-quark, quark-Wilson line, and Wilson
line-Wilson line.

We take Fig. 1(e) as the example to illustrate the calculation, in which the
one-loop corrections are

$\widetilde{O}_{e}=\left(\psi_{1}\left(z_{1}\right)\left(ig_{s}\int
d^{d}\eta_{1}\bar{\psi}_{1}\left(\eta_{1}\right)\not{A}\left(\eta_{1}\right)\psi_{1}\left(\eta_{1}\right)\right)\right){}^{T}\left(-ig_{s}\int_{0}^{1}dt_{1}z_{1}\cdot
A\left(t_{1}z_{1}\right)\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0).$
(56)

The color indexes and the parameter
$\left(\displaystyle\frac{\mu^{2}}{e^{\ln(4\pi)-\gamma_{E}}}\right)^{\epsilon}$
are not written out explicitly. The gluon and quark propagators in the
coordinate space are

$\displaystyle
G(x-y)=\frac{\Gamma(d/2-1)}{4\pi^{d/2}}\frac{-g_{\mu\nu}}{\left(-(x-y)^{2}+i\epsilon\right)^{d/2-1}},$
(57) $\displaystyle
Q(x-y)=\frac{\Gamma(d/2)}{2\pi^{d/2}}\frac{i(\not{x}-\not{y})}{\left(-(x-y)^{2}+i\epsilon\right)^{d/2}}.$
(58)

Following the tedious but standard routine: Schwinger parameterization, plus
field Fourier transformation, completing the square, parameter shifting, Dirac
algebra simplification, and parameter integration, one arrives at

$\displaystyle\widetilde{O}_{e}=g_{s}^{2}\frac{(-i)^{d/2-1}}{8\pi^{d/2}}\int
d^{d}k_{1}\int_{0}^{1}dt_{1}\int_{0}^{\infty}d\sigma_{1}\int_{0}^{\infty}d\sigma_{2}\sigma_{1}^{d/2-1}\sigma_{2}^{d/2-2}\left(\sigma_{1}+\sigma_{2}\right)^{-d/2}$
(59) $\displaystyle\times
e^{\frac{i\left(4\left(\sigma_{1}+\sigma_{2}t_{1}\right)z_{1}(k_{1}\cdot
n_{z})+k_{1}^{2}-4\sigma_{1}\sigma_{2}\left(t_{1}-1\right)^{2}(-z_{1}^{2}))\right)}{4\left(\sigma_{1}+\sigma_{2}\right)}}\psi_{1}^{T}\left(k_{1}\right)\left((-z_{1}^{2})-\frac{z_{1}(k_{1}\cdot
n_{z})+\left(\sigma_{1}+\sigma_{2}t_{1}\right)(-z_{1}^{2})}{\sigma_{1}+\sigma_{2}}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$

where $\sigma_{1}$ and $\sigma_{2}$ are Schwinger parameters, and $k_{1}$ is
from the Fourier transformation of $\psi_{1}^{T}(z_{1})$. Note that terms like
$k_{1}^{2}$ or $\not{k}_{1}\psi(z)$ have been neglected in the calculation due
to the equation of motion. By changing $(\sigma_{1},\sigma_{2})$ to
$(\sigma,\eta)$ with $\sigma_{1}=\displaystyle\frac{\sigma}{\eta}$ and
$\sigma_{2}=\displaystyle\frac{\sigma}{1-\eta}$, the above result can be
rearranged as

$\displaystyle\widetilde{O}_{e}$ $\displaystyle=$ $\displaystyle-
g_{s}^{2}\frac{(-i)^{d/2-1}}{8\pi^{d/2}}\int
d^{d}k_{1}\int_{0}^{1}dt_{1}\int_{0}^{1}d\eta\int_{0}^{\infty}d\sigma\sigma^{\frac{d}{2}-3}\left(\left(1-\eta\right)z_{1}(k_{1}\cdot
n_{z})+\sigma\left(t_{1}-1\right)(-z_{1}^{2})\right)$ (60)
$\displaystyle\times
e^{-iz_{1}n_{z}\cdot\left(k_{1}\left(-\left(\eta\left(t_{1}-1\right)\right)-1\right)+\sigma\left(t_{1}-1\right)^{2}z_{1}n_{z}\right)}\psi_{1}^{T}\left(k_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$

and then one can calculate the two parts.

For the term involving $(-z_{1}^{2})$ in Eq. (60), we further define

$\displaystyle\widetilde{O}_{e1}=$ $\displaystyle-
g_{s}^{2}\frac{1}{8\pi^{d/2}}\Gamma(d/2-1)\int
d^{d}k_{1}\int_{0}^{1}dt_{1}\int_{0}^{1}d\eta\left(1-t_{1}\right)^{3-d}\left(z_{1}^{2}\right)^{2-d/2}$
(61) $\displaystyle\times
e^{iz_{1}k_{1}\left(\eta\left(t_{1}-1\right)+1\right)}\psi_{1}^{T}\left(k_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$

and we can have the simplified form

$\displaystyle\widetilde{O}_{e1}=-g_{s}^{2}\frac{1}{8\pi^{d/2}}\Gamma(d/2-1)\left(z_{1}^{2}\right)^{2-d/2}$
$\displaystyle\int_{0}^{1}dt_{0}\int_{0}^{1}d\eta\left(t_{0}\right)^{3-d}\psi_{1}^{T}\left(\left(1-\eta
t_{0}\right)z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$ (62)

with $t_{0}=1-t_{1}$. The $t_{0}\to 0$ corresponds to a UV divergence since
that divergence is regularized by $d<4$ and one end of the Wilson line
approaches $z_{1}$ when $t_{0}\to 0$. One can separate this divergence from
the rest by using $\psi^{T}((1-\eta t_{0})z_{1})=\left(\psi^{T}((1-\eta
t_{0})z_{1})-\psi^{T}(z_{1})\right)+\psi^{T}(z_{1})$. Then it is
straightforward to obtain the results for these two parts

$\displaystyle\widetilde{O}_{e11}=$
$\displaystyle\displaystyle\frac{\alpha_{s}C_{F}}{4\pi}\left(\frac{1}{\epsilon_{\text{UV}}}+\log\left(\frac{1}{4}\mu^{2}z_{1}^{2}e^{2\gamma_{E}}\right)\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$
(63) $\displaystyle\widetilde{O}_{e12}=$
$\displaystyle\displaystyle\frac{\alpha_{s}C_{F}}{2\pi}\int_{0}^{1}d\eta\left(\frac{1-\eta}{\eta}\right)_{+}(\psi_{1}^{T}\left((1-\eta)z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0).$
(64)

For the $z_{1}(k_{1}\cdot n_{z})$ term in Eq. (60)

$\displaystyle\widetilde{O}_{e2}=$ $\displaystyle-
g_{s}^{2}\frac{(-i)^{d/2-1}}{8\pi^{d/2}}\int
d^{d}k_{1}\int_{0}^{1}dt_{1}\int_{0}^{1}d\eta\int_{0}^{\infty}d\sigma\sigma^{\frac{d}{2}-3}\left(1-\eta\right)z_{1}(k_{1}\cdot
n_{z})$ (65) $\displaystyle\times
e^{-iz_{1}n_{z}\cdot\left(k_{1}\left(-\left(\eta\left(t_{1}-1\right)\right)-1\right)+\sigma\left(t_{1}-1\right)^{2}z_{1}n_{z}\right)}\psi_{1}^{T}\left(k_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$

there is an IR divergence:

$\displaystyle\begin{aligned}
\widetilde{O}_{e2}=&-C_{F}\frac{g_{s}^{2}}{8\pi^{2}}\int_{0}^{1}d\eta\left(\left(\ln\left(\frac{1}{4}\mu^{2}z_{1}^{2}e^{2\gamma_{E}}\right)+\frac{1}{\epsilon_{\mathrm{IR}}}+2\right)\left(\frac{1-\eta}{\eta}\right)_{+}+\left(\frac{2\ln\eta}{\eta}\right)_{+}\right)\\\
&\times\psi_{1}^{T}\left(z_{1}(1-\eta)\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0).\end{aligned}$
(66)

Collecting all these pieces and removing the UV divergence in the
$\overline{\rm MS}$ scheme give the final result:

$\displaystyle\widetilde{O}_{e}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}\ln\left(\frac{1}{4}\mu_{\text{UV}}^{2}z_{1}^{2}e^{2\gamma_{E}}\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
(67)
$\displaystyle-\frac{\alpha_{s}C_{F}}{2\pi}\int_{0}^{1}d\eta\left(\frac{1-\eta}{\eta}\right)_{+}\left(\ln\left(\frac{1}{4}\mu_{\text{IR}}^{2}z_{1}^{2}e^{2\gamma_{E}}\right)+\frac{1}{\epsilon_{\text{IR}}}+1\right)\psi_{1}^{T}\left((1-\eta)z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
$\displaystyle-\frac{\alpha_{s}C_{F}}{\pi}\int_{0}^{1}d\eta\left(\frac{\ln\eta}{\eta}\right)_{+}\psi_{1}^{T}\left((1-\eta)z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$

where $\displaystyle\alpha_{s}=\frac{g^{2}_{s}}{4\pi}$. We have checked that
after making a Fourier transformation, the above results are consistent with
Ref. Deng:2023csv  in momentum space.

Then, in the same manner, results for the quark-Wilson-line diagrams are
derived as:

$\displaystyle\widetilde{O}_{d}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{8}\left(L_{12}^{\text{UV}}-L_{1}^{\text{UV}}\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
$\displaystyle-\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left((1-\eta)z_{1}+\eta
z_{2}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)\left\\{\left(L_{12}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}$
$\displaystyle+\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left((1-\eta)z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)\left\\{\left(L_{1}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\},$
(68)

$\displaystyle\widetilde{O}_{g}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{8}\left(L_{12}^{\text{UV}}-L_{2}^{\text{UV}}\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
$\displaystyle-\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(\eta
z_{1}+(1-\eta)z_{2}\right)\psi_{2}\left(z_{2}\right)\psi_{3}(0)\left\\{\left(L_{12}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}$
$\displaystyle+\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left((1-\eta)z_{2}\right)\psi_{3}(0)\left\\{\left(L_{2}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\},$
(69)

$\displaystyle\widetilde{O}_{e}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}L_{1}^{\text{UV}}\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
(70)
$\displaystyle-\frac{\alpha_{s}C_{F}}{2\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left((1-\eta)z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)\left\\{\left(L_{1}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\},$
$\displaystyle\widetilde{O}_{h}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{8\pi}L_{1}^{\text{UV}}\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
(71)
$\displaystyle-\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}\left(\eta
z_{1}\right)\left\\{\left(L_{1}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\},$
$\displaystyle\widetilde{O}_{f}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{4\pi}L_{2}^{\text{UV}}\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
(72)
$\displaystyle-\frac{\alpha_{s}C_{F}}{2\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left((1-\eta)z_{2}\right)\psi_{3}(0)\left\\{\left(L_{2}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\},$
$\displaystyle\widetilde{O}_{i}$
$\displaystyle=\frac{\alpha_{s}C_{F}}{8\pi}L_{2}^{\text{UV}}\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0)$
(73)
$\displaystyle-\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}\left(\eta
z_{2}\right)\left\\{\left(L_{2}^{\text{IR}}+1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\left(\frac{1-\eta}{\eta}\right)_{+}+2\left(\frac{\ln\eta}{\eta}\right)_{+}\right\\}.$

It should be mentioned that there are both IR and UV singularities in cases
(e, h, f, i), while there are no UV divergences for cases (d, g).

The last category is the Wilson line-Wilson line vertex, corresponding to Fig-
1 (j, k, l). After similar calculations, the one-loop results can be given as

$\displaystyle\widetilde{O}_{k}=\frac{\alpha_{s}C_{F}}{2\pi}\left(L_{1}^{\text{UV}}+2\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$
(74)
$\displaystyle\widetilde{O}_{l}=\frac{\alpha_{s}C_{F}}{2\pi}\left(L_{2}^{\text{UV}}+2\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0),$
(75)
$\displaystyle\widetilde{O}_{j}=-\frac{\alpha_{s}C_{F}}{4\pi}\left(L_{1}^{\text{UV}}+L_{2}^{\text{UV}}-L_{12}^{\text{UV}}+2\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}(0).$
(76)

Note that in these cases, only UV divergences arise when the two ends of the
gluons coincide with each other. The cases (e, f) and cases (d, g, h, i) have
different color coefficients Deng:2023csv . These color differences will also
appear in the following cases. More precisely, for cases (e, f, k, l), the
color algebra gives the same results as the meson case, which is $C_{F}$. For
other cases, the color parameter is $-\displaystyle\frac{C_{F}}{2}$.

The last pattern, quark-quark, is shown in Fig-1 (a, b, c). Following the same
routine, these results can be written down directly

$\displaystyle\widetilde{O}_{a}=$
$\displaystyle-\frac{\alpha_{s}C_{F}}{4\pi}\int_{0}^{1}d\eta_{1}\int_{0}^{1-\eta_{1}}d\eta_{2}\left(L_{12}^{\text{IR}}-3+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\psi_{1}^{T}\left(z_{1}\left(1-\eta_{1}\right)+z_{2}\eta_{1}\right)\Gamma\psi_{2}\left(z_{2}\left(1-\eta_{2}\right)+z_{1}\eta_{2}\right)\psi_{3}(0),$
(77) $\displaystyle\widetilde{O}_{b}=$
$\displaystyle-\frac{\alpha_{s}C_{F}}{8\pi}\int_{0}^{1}d\eta_{1}\int_{0}^{1-\eta_{1}}d\eta_{2}\left(L_{1}^{\text{IR}}-1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\psi_{1}^{T}\left((1-\eta_{1})z_{1}\right)\Gamma\psi_{2}\left(z_{2}\right)\psi_{3}\left(\eta_{2}z_{1}\right),$
(78) $\displaystyle\widetilde{O}_{c}=$
$\displaystyle-\frac{\alpha_{s}C_{F}}{8\pi}\int_{0}^{1}d\eta_{1}\int_{0}^{1-\eta_{1}}d\eta_{2}\left(L_{2}^{\text{IR}}-1+\frac{1}{\epsilon_{\mathrm{IR}}}\right)\psi_{1}^{T}\left(z_{1}\right)\Gamma\psi_{2}\left((1-\eta_{1})z_{2}\right)\psi_{3}\left(\eta_{2}z_{2}\right).$
(79)

The UV divergences have been subtracted in all these cases. It should be
mentioned that case (a) has an extra finite part compared to case (b) and case
(c).

Putting all these contributions together and sandwiching them between the
vacuum state $\left\langle 0|\right.$ and the lowest-order Fock state
$\left.|uds\right\rangle$, one can obtain the results in Eq. (2.2).

Additionally, the one-loop correction for the local matrix elements can also
be given. In this case, only the quark-quark category needs to be considered.
Since the $z^{2}$ is equal to $0$ in the local case, the local matrix elements
of quasi-DA and LCDA can be obtained by a similar calculation procedure. The
corresponding local matrix element can be given as

$\displaystyle{M}_{p}(0,0,0,P^{z},\mu)=\left(1-\frac{\alpha_{s}C_{F}}{4\pi}\frac{1}{\epsilon_{\rm
IR}}\right){M}_{0}\left(0,0,0,P^{z},\mu\right),$ (80)

in which the UV divergences have been subtracted.

##  Appendix B Explicit expressions for the hybrid counterterm

In the hybrid renormalization scheme, the counterterm in the matching kernel
can be split into different regions, such as $I_{\rm H}$, $I_{\rm HSI}$ and
$I_{\rm S}$ in Eq. (3.4). In this appendix, we provide the explicit
expressions for these terms. Those expressions can also be found in a
Mathematica notebook in the source package on arXiv.

We first give the necessary master formulas:

$\displaystyle\text{I1}\left(\left\\{L_{2},L_{1}\right\\},p\right)\equiv$
$\displaystyle\int_{L_{1}}^{L_{2}}\frac{dz}{2\pi}e^{ipz}\ln[z^{2}]$ (81)
$\displaystyle=$
$\displaystyle-\frac{i\left(2\left(\gamma_{E}+\log\left(-iL_{2}p\right)\right)+\left(-1+e^{iL_{2}p}\right)\log\left(L_{2}^{2}\right)+2\Gamma\left(0,-ipL_{2}\right)\right)}{2\pi
p}$
$\displaystyle+\frac{i\left(2\left(\gamma_{E}+\log\left(-iL_{1}p\right)\right)+\left(-1+e^{iL_{1}p}\right)\log\left(L_{1}^{2}\right)+2\Gamma\left(0,-ipL_{1}\right)\right)}{2\pi
p},$
$\text{I0}\left(\left\\{L_{2},L_{1}\right\\},p\right)\equiv\int_{L_{1}}^{L_{2}}\frac{dz}{2\pi}e^{ipz}=-\frac{i\left(-1+e^{iL_{2}p}\right)}{2\pi
p}+\frac{i\left(-1+e^{iL_{1}p}\right)}{2\pi p}.$ (82)

One can obtain the following results based on the above master integrals for
convenience ($L>0$)

$\displaystyle\text{I1t}\left(\left\\{L,-L\right\\},p\right)$
$\displaystyle\equiv
i[\text{I1}\left(\left\\{L,0\right\\},p\right)-\text{I1}\left(\left\\{0,-L\right\\},p\right)]$
(83)
$\displaystyle=\frac{2(-\text{Ci}(Lp)+\log(L)\cos(Lp)+\log(p)+\gamma_{E})}{\pi
p},$
$\text{I0}\left(\left\\{\infty,L\right\\},p\right)=\frac{\delta(p)}{2}+\frac{i}{2\pi
p}+\frac{i\left(-1+e^{iLp}\right)}{2\pi p},$ (84)
$\text{I0}\left(\left\\{-L,-\infty\right\\},p\right)=-\frac{i\left(-1+e^{-iLp}\right)}{2\pi
p}+\frac{\delta(p)}{2}-\frac{i}{2\pi p},$ (85)
$\text{I0}\left(\left\\{\infty,-\infty\right\\},p\right)=\delta(p).$ (86)

Based upon the above formulas, all the integrals in the Eq. (3.4) can be
written down ($z_{s}>0$).

$\displaystyle I_{\rm
H}[p_{1},p_{2}]\equiv\int\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{ip_{1}z_{1}+ip_{2}z_{2}}\left[\frac{7}{8}\ln\left(z_{1}^{2}\right)+\frac{7}{8}\ln\left(z_{2}^{2}\right)+\frac{3}{4}\ln\left((z_{1}-z_{2})^{2}\right)\right]$
(87)
$\displaystyle\quad\quad\quad\quad\times\left(\theta(2z_{s}-|z_{1}|)\theta(z_{s}-|z_{2}|)+\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-z_{s})\theta(2z_{s}-|z_{2}|)\right)$
$\displaystyle=\frac{7}{8}\left[\text{I0}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{2}\right)\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)+\left(\text{I0}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{1}\right)-\text{I0}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)\right)\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},p_{2}\right)\right.$
$\displaystyle\left.+\text{I0}\left(\left\\{z_{s},-z_{s}\right\\},p_{2}\right)\left(\text{I1}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{1}\right)-\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)\right)+\text{I0}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)\text{I1}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{2}\right)\right]$
$\displaystyle+\frac{3}{8\pi\left(p_{1}+p_{2}\right)}\left[\sin\left(\left(p_{1}+p_{2}\right)z_{s}\right)\left[\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)+\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},-p_{2}\right)\right.\right.$
$\displaystyle\left.\left.-\text{I1}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{1}\right)-\text{I1}\left(\left\\{2z_{s},-2z_{s}\right\\},-p_{2}\right)+\text{I1}\left(\left\\{3z_{s},-3z_{s}\right\\},p_{1}\right)+\text{I1}\left(\left\\{3z_{s},-3z_{s}\right\\},-p_{2}\right)\right]\right.$
$\displaystyle\left.+\sin\left(2\left(p_{1}+p_{2}\right)z_{s}\right)\left[-\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)-\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},-p_{2}\right)\right.\right.$
$\displaystyle\left.\left.+\text{I1}\left(\left\\{3z_{s},-3z_{s}\right\\},p_{1}\right)+\text{I1}\left(\left\\{3z_{s},-3z_{s}\right\\},-p_{2}\right)\right]\right.$
$\displaystyle\left.+\cos\left(2\left(p_{1}+p_{2}\right)z_{s}\right)\left[-\text{I1t}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)+\text{I1t}\left(\left\\{z_{s},-z_{s}\right\\},-p_{2}\right)\right.\right.$
$\displaystyle\left.\left.+\text{I1t}\left(\left\\{3z_{s},-3z_{s}\right\\},p_{1}\right)-\text{I1t}\left(\left\\{3z_{s},-3z_{s}\right\\},-p_{2}\right)\right]\right.$
$\displaystyle\left.+\cos\left(\left(p_{1}+p_{2}\right)z_{s}\right)\left[-\text{I1t}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)+\text{I1t}\left(\left\\{z_{s},-z_{s}\right\\},-p_{2}\right)-\text{I1t}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{1}\right)\right.\right.$
$\displaystyle\left.\left.+\text{I1t}\left(\left\\{2z_{s},-2z_{s}\right\\},-p_{2}\right)+\text{I1t}\left(\left\\{3z_{s},-3z_{s}\right\\},p_{1}\right)-\text{I1t}\left(\left\\{3z_{s},-3z_{s}\right\\},-p_{2}\right)\right]\right],$

$\displaystyle I_{\rm
HSI}[p_{1},p_{2}]\equiv\int\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{ip_{1}z_{1}+ip_{2}z_{2}}\theta(z_{s}-|z_{1}|)\theta(|z_{2}|-2z_{s})$
(88)
$\displaystyle\quad\quad\quad\quad\,\,\,\times\left[\frac{7}{8}\ln\left(z_{1}^{2}\right)+\frac{7}{8}\ln\left((2z_{s})^{2}\right)+\frac{3}{4}\ln\left((z_{1}-2z_{s}{\rm
sign}[z_{2}])^{2}\right)\right]$
$\displaystyle=\frac{1}{8}\left[6e^{2ip_{1}z_{s}}\text{I1}\left(\left\\{-z_{s},-3z_{s}\right\\},p_{1}\right)\text{I0}\left(\left\\{\infty,2z_{s}\right\\},p_{2}\right)+6e^{-2ip_{1}z_{s}}\text{I1}\left(\left\\{3z_{s},z_{s}\right\\},p_{1}\right)\text{I0}\left(\left\\{-2z_{s},-\infty\right\\},p_{2}\right)\right.$
$\displaystyle\left.+7\left(\text{I0}\left(\\{\infty,-\infty\\},p_{2}\right)-\text{I0}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{2}\right)\right)\left(\log\left(4z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)+\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},p_{1}\right)\right)\right],$

$\displaystyle I_{\rm
HSII}[p_{1},p_{2}]\equiv\int\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{ip_{1}z_{1}+ip_{2}z_{2}}\theta(|z_{1}|-2z_{s})\theta(z_{s}-|z_{2}|)$
(89)
$\displaystyle\quad\quad\quad\quad\quad\times\left[\frac{7}{8}\ln\left((2z_{s})^{2}\right)+\frac{7}{8}\ln\left(z_{2}^{2}\right)+\frac{3}{4}\ln\left(({\rm
sign}[z_{1}]2z_{s}-z_{2})^{2}\right)\right]$
$\displaystyle=\frac{1}{8}\left[6e^{2ip_{2}z_{s}}\text{I1}\left(\left\\{-z_{s},-3z_{s}\right\\},p_{2}\right)\text{I0}\left(\left\\{\infty,2z_{s}\right\\},p_{1}\right)+6e^{-2ip_{2}z_{s}}\text{I1}\left(\left\\{3z_{s},z_{s}\right\\},p_{2}\right)\text{I0}\left(\left\\{-2z_{s},-\infty\right\\},p_{1}\right)\right.$
$\displaystyle\left.+7\left(\text{I0}\left(\\{\infty,-\infty\\},p_{1}\right)-\text{I0}\left(\left\\{2z_{s},-2z_{s}\right\\},p_{1}\right)\right)\left(\log\left(4z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},-z_{s}\right\\},p_{2}\right)+\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},p_{2}\right)\right)\right],$

$\displaystyle I_{\rm
HSIII}[p_{1},p_{2}]\equiv\int\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{ip_{1}z_{1}+ip_{2}z_{2}}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}-z_{2}|)$
(90)
$\displaystyle\times\left[\frac{7}{8}\ln\left(\left(z_{s}+\left(z_{1}-z_{2}\right)\theta\left(z_{1}-z_{2}\right)\right){}^{2}\right)+\frac{7}{8}\ln\left(\left(z_{s}+\left(z_{2}-z_{1}\right)\theta\left(z_{2}-z_{1}\right)\right){}^{2}\right)+\frac{3}{4}\ln\left(\left(z_{1}-z_{2}\right){}^{2}\right)\right]$
$\displaystyle=\frac{1}{8}\text{I0}\left(\\{\infty,-\infty\\},p_{1}+p_{2}\right)\left[7\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},-z_{s}\right\\},\frac{1}{2}\left(p_{1}-p_{2}\right)\right)+6\,\text{I1}\left(\left\\{z_{s},-z_{s}\right\\},\frac{1}{2}\left(p_{1}-p_{2}\right)\right)\right.$
$\displaystyle\left.+7e^{\frac{1}{2}i\left(p_{1}-p_{2}\right)z_{s}}\text{I1}\left(\left\\{-z_{s},-2z_{s}\right\\},\frac{1}{2}\left(p_{1}-p_{2}\right)\right)+7e^{-\frac{1}{2}i\left(p_{1}-p_{2}\right)z_{s}}\text{I1}\left(\left\\{2z_{s},z_{s}\right\\},\frac{1}{2}\left(p_{1}-p_{2}\right)\right)\right]$
$\displaystyle+\frac{ie^{-i\left(p_{1}+p_{2}\right)z_{s}}}{16\pi\left(p_{1}+p_{2}\right)}\left[-7\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{0,-z_{s}\right\\},p_{1}\right)+7e^{2i\left(p_{1}+p_{2}\right)z_{s}}\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{0,-z_{s}\right\\},-p_{2}\right)\right.$
$\displaystyle\left.+7e^{2i\left(p_{1}+p_{2}\right)z_{s}}\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},0\right\\},p_{1}\right)-7\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},0\right\\},-p_{2}\right)-6\text{I1}\left(\left\\{0,-z_{s}\right\\},p_{1}\right)\right.$
$\displaystyle\left.+6e^{2i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{0,-z_{s}\right\\},-p_{2}\right)-7e^{ip_{1}z_{s}}\text{I1}\left(\left\\{-z_{s},-2z_{s}\right\\},p_{1}\right)+7e^{i\left(2p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{-z_{s},-2z_{s}\right\\},-p_{2}\right)\right.$
$\displaystyle\left.+6e^{2i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{z_{s},0\right\\},p_{1}\right)-6\text{I1}\left(\left\\{z_{s},0\right\\},-p_{2}\right)+7e^{i\left(p_{1}+2p_{2}\right)z_{s}}\text{I1}\left(\left\\{2z_{s},z_{s}\right\\},p_{1}\right)-7e^{ip_{2}z_{s}}\text{I1}\left(\left\\{2z_{s},z_{s}\right\\},-p_{2}\right)\right],$

$\displaystyle I_{\rm
HSIV}[p_{1},p_{2}]\equiv\int\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{ip_{1}z_{1}+ip_{2}z_{2}}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(z_{s}-|z_{1}+z_{2}|)$
(91)
$\displaystyle\times\left[\frac{7}{8}\ln\left(\left(z_{s}+\left(z_{1}+z_{2}\right)\theta\left(z_{1}+z_{2}\right)\right){}^{2}\right)+\frac{7}{8}\ln\left(\left(-z_{s}+\left(z_{2}+z_{1}\right)\theta\left(-z_{2}-z_{1}\right)\right){}^{2}\right)+\frac{3}{4}\ln\left(\left(2z_{s}+|z_{1}+z_{2}|\right){}^{2}\right)\right]$
$\displaystyle=\frac{1}{16}\text{I0}\left(\\{\infty,-\infty\\},\frac{1}{2}\left(p_{1}-p_{2}\right)\right)\left[7\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},-z_{s}\right\\},\frac{1}{2}\left(p_{1}+p_{2}\right)\right)\right.$
$\displaystyle\left.+6e^{i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{-2z_{s},-3z_{s}\right\\},\frac{1}{2}\left(p_{1}+p_{2}\right)\right)+7e^{\frac{1}{2}i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{-z_{s},-2z_{s}\right\\},\frac{1}{2}\left(p_{1}+p_{2}\right)\right)\right.$
$\displaystyle\left.+7e^{-\frac{1}{2}i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{2z_{s},z_{s}\right\\},\frac{1}{2}\left(p_{1}+p_{2}\right)\right)+6e^{-i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{3z_{s},2z_{s}\right\\},\frac{1}{2}\left(p_{1}+p_{2}\right)\right)\right]$
$\displaystyle+\frac{ie^{-i\left(p_{1}+p_{2}\right)z_{s}}}{16\pi\left(p_{1}-p_{2}\right)}\left[-7e^{2ip_{2}z_{s}}\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{0,-z_{s}\right\\},p_{1}\right)+7e^{2ip_{1}z_{s}}\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{0,-z_{s}\right\\},p_{2}\right)\right.$
$\displaystyle\left.+7e^{2ip_{1}z_{s}}\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},0\right\\},p_{1}\right)-7e^{2ip_{2}z_{s}}\log\left(z_{s}^{2}\right)\text{I0}\left(\left\\{z_{s},0\right\\},p_{2}\right)-6e^{2i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{-2z_{s},-3z_{s}\right\\},p_{1}\right)\right.$
$\displaystyle\left.+6e^{2i\left(p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{-2z_{s},-3z_{s}\right\\},p_{2}\right)-7e^{i\left(p_{1}+2p_{2}\right)z_{s}}\text{I1}\left(\left\\{-z_{s},-2z_{s}\right\\},p_{1}\right)+7e^{i\left(2p_{1}+p_{2}\right)z_{s}}\text{I1}\left(\left\\{-z_{s},-2z_{s}\right\\},p_{2}\right)\right.$
$\displaystyle\left.+7e^{ip_{1}z_{s}}\text{I1}\left(\left\\{2z_{s},z_{s}\right\\},p_{1}\right)-7e^{ip_{2}z_{s}}\text{I1}\left(\left\\{2z_{s},z_{s}\right\\},p_{2}\right)+6\text{I1}\left(\left\\{3z_{s},2z_{s}\right\\},p_{1}\right)-6\text{I1}\left(\left\\{3z_{s},2z_{s}\right\\},p_{2}\right)\right],$

$\displaystyle I_{\rm
S}[p_{1},p_{2}]\equiv\int\frac{dz_{1}}{2\pi}\frac{dz_{2}}{2\pi}e^{ip_{1}z_{1}+ip_{2}z_{2}}\theta(|z_{1}|-z_{s})\theta(|z_{2}|-z_{s})\theta(|z_{1}-z_{2}|-z_{s})\theta(|z_{1}+z_{2}|-z_{s})$
(92)
$\displaystyle\quad\quad\quad\quad\times\left[\frac{7}{8}\ln\left(z_{s}^{2}\right)+\frac{7}{8}\ln\left((2z_{s})^{2}\right)+\frac{3}{4}\ln\left(({\rm
sign}[z_{1}]z_{s}-{\rm sign}[z_{2}]2z_{s})^{2}\right)\right]$
$\displaystyle=\delta\left(p_{1}-p_{2}\right)\delta\left(p_{1}+p_{2}\right)\left(\frac{7}{4}\log\left(4z_{s}^{4}\right)+\frac{3}{4}\log\left(9z_{s}^{4}\right)\right)$
$\displaystyle-\frac{\delta\left(p_{1}+p_{2}\right)\left(6\log\left(z_{s}^{2}\right)+7\log\left(4z_{s}^{4}\right)\right)\sin\left(\frac{1}{2}\left(p_{1}-p_{2}\right)z_{s}\right)}{4\pi\left(p_{1}-p_{2}\right)}$
$\displaystyle-\frac{\delta\left(p_{1}-p_{2}\right)\left(6\log\left(9z_{s}^{2}\right)+7\log\left(4z_{s}^{4}\right)\right)\sin\left(\frac{1}{2}\left(p_{1}+p_{2}\right)z_{s}\right)}{4\pi\left(p_{1}+p_{2}\right)}$
$\displaystyle-\frac{\delta\left(p_{2}\right)\left(20\log\left(z_{s}\right)+\log(3456)\right)\sin\left(p_{1}z_{s}\right)}{4\pi
p_{1}}-\frac{\delta\left(p_{1}\right)\left(20\log\left(z_{s}\right)+\log(3456)\right)\sin\left(p_{2}z_{s}\right)}{4\pi
p_{2}}$
$\displaystyle-\frac{\left(20\log\left(z_{s}\right)+\log(128)\right)\left(p_{1}\cos\left(\left(p_{1}+2p_{2}\right)z_{s}\right)+p_{2}\cos\left(\left(2p_{1}+p_{2}\right)z_{s}\right)\right)}{8\pi^{2}p_{1}p_{2}\left(p_{1}+p_{2}\right)}$
$\displaystyle+\frac{\left(20\log\left(z_{s}\right)+\log(93312)\right)\left(p_{1}\cos\left(\left(p_{1}-2p_{2}\right)z_{s}\right)-p_{2}\cos\left(2p_{1}z_{s}-p_{2}z_{s}\right)\right)}{8\pi^{2}p_{1}\left(p_{1}-p_{2}\right)p_{2}}.$

##  Appendix C Double plus function

In this appendix, we demonstrate that the matching kernel at the one-loop
level can be properly expressed in the form of the double plus function in Eq.
(55). This is equivalent to proving the following lemma:
Lemma. Consider a generalized function $G(x_{1},x_{2},x_{10},x_{20})$
satisfying the following properties

$\displaystyle\lim_{|x_{1}-x_{10}|\rightarrow+\infty}|x_{1}-x_{10}|G(x_{1},x_{2},x_{10},x_{20})=_{a.e.}0,$
(93)
$\displaystyle\lim_{|x_{2}-x_{20}|\rightarrow+\infty}|x_{2}-x_{20}|G(x_{1},x_{2},x_{10},x_{20})=_{a.e.}0,$
$\displaystyle\lim_{|x_{3}-x_{30}|\rightarrow+\infty}|x_{3}-x_{30}|G(x_{1},x_{2},x_{10},x_{20})=_{a.e.}0,$
$\displaystyle\exists\,\epsilon>0,\,\left|\int_{x_{1}-\epsilon}^{x_{1}+\epsilon}dx_{10}\,|x_{1}-x_{10}|G(x_{1},x_{2},x_{10},x_{20})\right|<+\infty,$
(94)
$\displaystyle\exists\,\epsilon>0,\,\left|\int_{x_{2}-\epsilon}^{x_{2}+\epsilon}dx_{20}\,|x_{2}-x_{20}|G(x_{1},x_{2},x_{10},x_{20})\right|<+\infty,$
$\displaystyle\exists\,\epsilon>0,\,\left|\int_{x_{3}-\epsilon}^{x_{3}+\epsilon}dx_{30}\,|x_{3}-x_{30}|G(x_{1},x_{2},x_{10},x_{20})\right|<+\infty,$

and in a specific form

$\displaystyle G(x_{1},x_{2},x_{10},x_{20})=$ $\displaystyle
A(x_{1},x_{2},x_{10},x_{20})+\delta(x_{1}-x_{10})B(x_{2},x_{10},x_{20})+\delta(x_{2}-x_{20})C(x_{1},x_{10},x_{20})$
(95) $\displaystyle+$
$\displaystyle\delta(x_{3}-x_{30})D(x_{1},x_{2},x_{10},x_{20})+\delta(x_{1}-x_{10})\delta(x_{2}-x_{20})E(x_{10},x_{20}),$

where $x_{3}(x_{30})$ is shorthand for $1-x_{1}-x_{2}(1-x_{10}-x_{20})$ and
$A,B,C,D,E$ can be expanded into power series with logarithmic terms without
Dirac delta functions. $``a.e."$ is shorthand for “almost everywhere” which
means the results hold except for $x_{1}=x_{10}$, $x_{2}=x_{20}$ or
$x_{3}=x_{30}$ here. The lemma states that if the integral $\int
dx_{1}dx_{2}G(x_{1},x_{2},x_{10},x_{20})=0$, the generalized function
$G(x_{1},x_{2},x_{10},x_{20})$ can be written as the double plus function
defined in Eq. (18).

Proof. Based on the condition that $\int
dx_{1}dx_{2}G(x_{1},x_{2},x_{10},x_{20})=0$, we obtain

$\displaystyle E(x_{10},x_{20})=$ $\displaystyle-\int
A(x_{1},x_{2},x_{10},x_{20})dx_{1}dx_{2}-\int
B(x_{2},x_{10},x_{20})dx_{2}-\int C(x_{1},x_{10},x_{20})dx_{1}$ (96)
$\displaystyle-\int\delta(x_{3}-x_{30})D(x_{1},x_{2},x_{10},x_{20})dx_{1}dx_{2}.$

Plugging $E(x_{10},x_{20})$ into Eq. (95) gives:

$\displaystyle
G(x_{1},x_{2},x_{10},x_{20})=A(x_{1},x_{2},x_{10},x_{20})-\delta(x_{1}-x_{10})\delta(x_{2}-x_{20})\int
A(y_{1},y_{2},x_{10},x_{20})dy_{1}dy_{2}$ (97)
$\displaystyle+\delta(x_{1}-x_{10})B(x_{2},x_{10},x_{20})-\delta(x_{1}-x_{10})\delta(x_{2}-x_{20})\int\delta(y_{1}-x_{10})B(y_{2},x_{10},x_{20})dy_{1}dy_{2}$
$\displaystyle+\delta(x_{2}-x_{20})C(x_{1},x_{10},x_{20})-\delta(x_{1}-x_{10})\delta(x_{2}-x_{20})\int\delta(y_{2}-x_{20})C(y_{1},x_{10},x_{20})dy_{1}dy_{2}$
$\displaystyle+\delta(x_{3}-x_{30})D(x_{1},x_{2},x_{10},x_{20})-\delta(x_{1}-x_{10})\delta(x_{2}-x_{20})\int\delta(y_{3}-x_{30})D(y_{1},y_{2},x_{10},x_{20})dy_{1}dy_{2}$
$\displaystyle=\left[A(x_{1},x_{2},x_{10},x_{20})+\delta(x_{1}-x_{10})B(x_{2},x_{10},x_{20})+\delta(x_{2}-x_{20})C(x_{1},x_{10},x_{20})\right.$
$\displaystyle\left.+\delta(x_{3}-x_{30})D(x_{1},x_{2},x_{10},x_{20})\right]_{\oplus},$

which is the double plus function defined in Eq. (18). Because of Eqs. (93)
and (94), one obtains a finite result after convoluting the generalized
function $G(x_{1},x_{2},x_{10},x_{20})$ with a Schwartz function
$\Phi(x_{10},x_{20})$,

$\int dx_{10}dx_{20}G(x_{1},x_{2},x_{10},x_{20})\Phi(x_{10},x_{20})<\infty.$
(98)

## References

  * [1] H. H. Shih, S. C. Lee and H. n. Li, Phys. Rev. D 59, 094014 (1999) doi:10.1103/PhysRevD.59.094014 [arXiv:hep-ph/9810515 [hep-ph]].
  * [2] R. Aaij et al. [LHCb], Nature Phys. 11, 743-747 (2015) doi:10.1038/nphys3415 [arXiv:1504.01568 [hep-ex]].
  * [3] R. Aaij et al. [LHCb], JHEP 06, 115 (2015) [erratum: JHEP 09, 145 (2018)] doi:10.1007/JHEP06(2015)115 [arXiv:1503.07138 [hep-ex]].
  * [4] R. Aaij et al. [LHCb], Phys. Rev. D 105, no.5, L051104 (2022) doi:10.1103/PhysRevD.105.L051104 [arXiv:2111.10194 [hep-ex]].
  * [5] V. L. Chernyak and I. R. Zhitnitsky, Nucl. Phys. B 246, 52-74 (1984) doi:10.1016/0550-3213(84)90114-7
  * [6] I. D. King and C. T. Sachrajda, Nucl. Phys. B 279, 785-803 (1987) doi:10.1016/0550-3213(87)90019-8
  * [7] V. M. Braun, S. E. Derkachov, G. P. Korchemsky and A. N. Manashov, Nucl. Phys. B 553, 355-426 (1999) doi:10.1016/S0550-3213(99)00265-5 [arXiv:hep-ph/9902375 [hep-ph]].
  * [8] I. V. Anikin, V. M. Braun and N. Offen, Phys. Rev. D 88, 114021 (2013) doi:10.1103/PhysRevD.88.114021 [arXiv:1310.1375 [hep-ph]].
  * [9] M. Gockeler, R. Horsley, T. Kaltenbrunner, Y. Nakamura, D. Pleiter, P. E. L. Rakow, A. Schafer, G. Schierholz, H. Stuben and N. Warkentin, et al. Phys. Rev. Lett. 101, 112002 (2008) doi:10.1103/PhysRevLett.101.112002 [arXiv:0804.1877 [hep-lat]].
  * [10] V. M. Braun et al. [QCDSF], Phys. Rev. D 79, 034504 (2009) doi:10.1103/PhysRevD.79.034504 [arXiv:0811.2712 [hep-lat]].
  * [11] G. S. Bali, V. M. Braun, M. Göckeler, M. Gruber, F. Hutzler, A. Schäfer, R. W. Schiel, J. Simeth, W. Söldner and A. Sternbeck, et al. JHEP 02, 070 (2016) doi:10.1007/JHEP02(2016)070 [arXiv:1512.02050 [hep-lat]].
  * [12] G. S. Bali et al. [RQCD], Eur. Phys. J. A 55, no.7, 116 (2019) doi:10.1140/epja/i2019-12803-6 [arXiv:1903.12590 [hep-lat]].
  * [13] Z. F. Deng, C. Han, W. Wang, J. Zeng and J. L. Zhang, JHEP 07, 191 (2023) doi:10.1007/JHEP07(2023)191 [arXiv:2304.09004 [hep-ph]].
  * [14] X. Ji, Phys. Rev. Lett. 110, 262002 (2013) doi:10.1103/PhysRevLett.110.262002 [arXiv:1305.1539 [hep-ph]].
  * [15] X. Ji, Sci. China Phys. Mech. Astron. 57, 1407-1412 (2014) doi:10.1007/s11433-014-5492-3 [arXiv:1404.6680 [hep-ph]].
  * [16] K. Cichy and M. Constantinou, Adv. High Energy Phys. 2019, 3036904 (2019) doi:10.1155/2019/3036904 [arXiv:1811.07248 [hep-lat]].
  * [17] Y. Zhao, Int. J. Mod. Phys. A 33, no.36, 1830033 (2019) doi:10.1142/S0217751X18300338 [arXiv:1812.07192 [hep-ph]].
  * [18] X. Ji, Y. S. Liu, Y. Liu, J. H. Zhang and Y. Zhao, Rev. Mod. Phys. 93, no.3, 035005 (2021) doi:10.1103/RevModPhys.93.035005 [arXiv:2004.03543 [hep-ph]].
  * [19] G. Martinelli, C. Pittori, C. T. Sachrajda, M. Testa and A. Vladikas, Nucl. Phys. B 445, 81-108 (1995) doi:10.1016/0550-3213(95)00126-D [arXiv:hep-lat/9411010 [hep-lat]].
  * [20] K. Zhang et al. [$\chi$QCD], Phys. Rev. D 104, no.7, 074501 (2021) doi:10.1103/PhysRevD.104.074501 [arXiv:2012.05448 [hep-lat]].
  * [21] Y. K. Huo et al. [Lattice Parton Collaboration (LPC)], Nucl. Phys. B 969, 115443 (2021) doi:10.1016/j.nuclphysb.2021.115443 [arXiv:2103.02965 [hep-lat]].
  * [22] X. Ji, Y. Liu, A. Schäfer, W. Wang, Y. B. Yang, J. H. Zhang and Y. Zhao, Nucl. Phys. B 964, 115311 (2021) doi:10.1016/j.nuclphysb.2021.115311 [arXiv:2008.03886 [hep-ph]].
  * [23] J. Hua et al. [Lattice Parton], Phys. Rev. Lett. 127, no.6, 062002 (2021) doi:10.1103/PhysRevLett.127.062002 [arXiv:2011.09788 [hep-lat]].
  * [24] X. Gao, A. D. Hanlon, S. Mukherjee, P. Petreczky, P. Scior, S. Syritsyn and Y. Zhao, Phys. Rev. Lett. 128, no.14, 142003 (2022) doi:10.1103/PhysRevLett.128.142003 [arXiv:2112.02208 [hep-lat]].
  * [25] J. Hua et al. [Lattice Parton], Phys. Rev. Lett. 129, no.13, 132001 (2022) doi:10.1103/PhysRevLett.129.132001 [arXiv:2201.09173 [hep-lat]].
  * [26] C. Y. Chou and J. W. Chen, Phys. Rev. D 106, no.1, 014507 (2022) doi:10.1103/PhysRevD.106.014507 [arXiv:2204.08343 [hep-lat]].
  * [27] J. Hua et al. [Lattice Parton (LPC)], PoS LATTICE2021, 322 (2022) doi:10.22323/1.396.0322
  * [28] X. Gao, A. D. Hanlon, N. Karthik, S. Mukherjee, P. Petreczky, P. Scior, S. Shi, S. Syritsyn, Y. Zhao and K. Zhou, Phys. Rev. D 106, no.11, 114510 (2022) doi:10.1103/PhysRevD.106.114510 [arXiv:2208.02297 [hep-lat]].
  * [29] F. Yao et al. [Lattice Parton], [arXiv:2208.08008 [hep-lat]].
  * [30] Y. Su, J. Holligan, X. Ji, F. Yao, J. H. Zhang and R. Zhang, Nucl. Phys. B 991, 116201 (2023) doi:10.1016/j.nuclphysb.2023.116201 [arXiv:2209.01236 [hep-ph]].
  * [31] X. Ji, [arXiv:2209.09332 [hep-lat]].
  * [32] X. Gao, A. D. Hanlon, S. Mukherjee, P. Petreczky, P. Scior, S. Syritsyn and Y. Zhao, PoS LATTICE2022, 104 (2023) doi:10.22323/1.430.0104
  * [33] X. Gao, A. D. Hanlon, J. Holligan, N. Karthik, S. Mukherjee, P. Petreczky, S. Syritsyn and Y. Zhao, Phys. Rev. D 107, no.7, 074509 (2023) doi:10.1103/PhysRevD.107.074509 [arXiv:2212.12569 [hep-lat]].
  * [34] Y. Ji, F. Yao and J. H. Zhang, [arXiv:2212.14415 [hep-ph]].
  * [35] R. Zhang, doi:10.13016/dspace/gggm-uvss
  * [36] J. Holligan, X. Ji, H. W. Lin, Y. Su and R. Zhang, Nucl. Phys. B 993, 116282 (2023) doi:10.1016/j.nuclphysb.2023.116282 [arXiv:2301.10372 [hep-lat]].
  * [37] R. Zhang, J. Holligan, X. Ji and Y. Su, Phys. Lett. B 844, 138081 (2023) doi:10.1016/j.physletb.2023.138081 [arXiv:2305.05212 [hep-lat]].
  * [38] X. Gao, W. Y. Liu and Y. Zhao, [arXiv:2306.14960 [hep-ph]].
  * [39] T. Izubuchi, X. Ji, L. Jin, I. W. Stewart and Y. Zhao, Phys. Rev. D 98, no.5, 056004 (2018) doi:10.1103/PhysRevD.98.056004 [arXiv:1801.03917 [hep-ph]].
  * [40] K. Orginos, A. Radyushkin, J. Karpie and S. Zafeiropoulos, Phys. Rev. D 96, no.9, 094503 (2017) doi:10.1103/PhysRevD.96.094503 [arXiv:1706.05373 [hep-ph]].
  * [41] A. V. Radyushkin, Phys. Rev. D 96, no.3, 034025 (2017) doi:10.1103/PhysRevD.96.034025 [arXiv:1705.01488 [hep-ph]].
  * [42] A. V. Radyushkin, Phys. Lett. B 781, 433-442 (2018) doi:10.1016/j.physletb.2018.04.023 [arXiv:1710.08813 [hep-ph]].
  * [43] J. W. Chen, X. Ji and J. H. Zhang, Nucl. Phys. B 915, 1-9 (2017) doi:10.1016/j.nuclphysb.2016.12.004 [arXiv:1609.08102 [hep-ph]].
  * [44] X. Ji, J. H. Zhang and Y. Zhao, Phys. Rev. Lett. 120, no.11, 112001 (2018) doi:10.1103/PhysRevLett.120.112001 [arXiv:1706.08962 [hep-ph]].
  * [45] T. Ishikawa, Y. Q. Ma, J. W. Qiu and S. Yoshida, Phys. Rev. D 96, no.9, 094019 (2017) doi:10.1103/PhysRevD.96.094019 [arXiv:1707.03107 [hep-ph]].
  * [46] J. Green, K. Jansen and F. Steffens, Phys. Rev. Lett. 121, no.2, 022004 (2018) doi:10.1103/PhysRevLett.121.022004 [arXiv:1707.07152 [hep-lat]].
  * [47] X. D. Ji, [arXiv:hep-ph/9507322 [hep-ph]].
  * [48] M. Beneke, Phys. Rept. 317, 1-142 (1999) doi:10.1016/S0370-1573(98)00130-6 [arXiv:hep-ph/9807443 [hep-ph]].
  * [49] C. Bauer, G. S. Bali and A. Pineda, Phys. Rev. Lett. 108, 242002 (2012) doi:10.1103/PhysRevLett.108.242002 [arXiv:1111.3946 [hep-ph]].
  * [50] G. S. Bali, C. Bauer, A. Pineda and C. Torrero, Phys. Rev. D 87, 094517 (2013) doi:10.1103/PhysRevD.87.094517 [arXiv:1303.3279 [hep-lat]].
  * [51] F. Karbstein, M. Wagner and M. Weber, Phys. Rev. D 98, no.11, 114506 (2018) doi:10.1103/PhysRevD.98.114506 [arXiv:1804.10909 [hep-ph]].