Pseudo and quasi quark PDF in the BFKL approximation

I examine the high-energy behavior of the Ioffe-time distribution for the quark bi-local space-like separated operator using the high-energy operator product expansion. These findings have significant implications for lattice calculations, which require extrapolation for large Ioffe-time values. I perform an explicit Fourier transform for both the pseudo-PDF and quasi-PDF, and investigate their behavior within the first two leading twist contributions. I show that the quark pseudo-PDF captures the BFKL resummation (resummation of all twists) and exhibits a rising behavior for small $x_B$ values, while the quasi-PDF presents a different behavior. I demonstrate that an appropriate small-$x_B$ behavior cannot be achieved solely through DGLAP dynamics, emphasizing the importance of all-twist resummation. This study provides valuable insights into quark non-local operators' high-energy behavior and the limitations of lattice calculations in this context.


Introduction
Over the past decade, the investigation of Euclidean-separated, gauge-invariant, bi-local operators via lattice gauge formalism has garnered significant interest due to its ability to provide direct access to parton distribution functions (PDFs) from first principles (for a review see [1][2][3][4]).The foundation for this surge in activity lies in the observation that space-like separated operators can be examined through lattice QCD formalism, and that in the infinite momentum frame, they reduce to the conventional light-cone operators through which PDFs are defined [5].Deviations from the infinite momentum frame emerge as inverse powers of the large parameter of the boost, suggesting that such corrections can be systematically suppressed.
The distributions introduced in [5], known as quasi-PDFs, were later complemented by alternative PDFs called pseudo-PDFs in [6].The Bjorken-x (x B ) dependence of these two distributions has been extensively studied in recent years [7][8][9][10][11]; however, lattice formalism is unlikely to provide access to their behavior at small x B values.
We are going to demonstrate (preliminary result have been presented in [12]) that despite being defined through the same space-like separated bi-local operators [13][14][15][16], quark quasi-and pseudo-PDFs exhibit markedly different behavior at small x B values.
The importance of understanding the small-x B behavior of quark pseudo-and quasi-PDFs cannot be overstated, even when the corresponding behavior for gluon pseudo-and quasi-PDFs has been established [17].This is primarily due to the fact that most lattice calculations focus on quark distributions rather than gluon ones.Comprehensive knowledge of small-x B behavior for quark distributions is essential for providing a more complete picture of the underlying dynamics in QCD, ultimately leading to more accurate predictions and a deeper understanding of the fundamental interactions within the strong force.Moreover, insights into the quark pseudo-and quasi-PDFs' behavior at small-x B values will be indispensable for the interpretation of experimental results from future QCD colliders probing this regime, such as the Electron-Ion Collider (EIC) in the USA [18], the EIC in China [19], and the Large Hadron electron Collider (LHeC) in Europe [20].
The small-x B behavior of deep inelastic scattering (DIS) structure functions can be calculated using the high-energy operator product expansion (OPE) [21], wherein the Tproduct of two electromagnetic currents is expanded in terms of coefficient functions, also known as impact factors, convoluted with the matrix elements of infinite Wilson lines.The evolution equation of the matrix elements of infinite Wilson lines with respect to the rapidity factorization parameter is the BK equation [21][22][23][24][25], which addresses smallx B leading-log resummation through its linear term and the unitarity property of the theory through its non-linear term.Due to its non-linear nature, the BK equation is a generalization of the BFKL equation.
The high-energy OPE is a complementary procedure to the OPE in non-local (finite gauge-link) operators adopted in the Bjorken limit, where the factorization scale is the transverse momentum of the fields and the evolution equation of the finite gauge-link with respect to the factorization parameter is the DGLAP evolution equation.
In the DGLAP regime, the structure functions of DIS are governed by the anomalous dimensions of twist-2 operators, while in the BFKL regime, the structure functions receive contributions from an infinite series of twist expansion.By moving beyond the BFKL limit to the overlapping region, one can obtain the anomalous dimension of twist-2 (and potentially higher twist) operators in the small-x B limit from the small-x B structure function.With the analytic continuation of the twist-2 local operator to non-integer spin, the local operators become non-local light-ray operators, making the prediction of the small-x B limit of the twist-2 anomalous dimension from the BFKL resummation more transparent [26].In this work, we demonstrate that this formalism enables us to obtain the leading (LT) and next-to-leading twist (NLT) corrections from the full BFKL resummed results for the pseudo-and quasi-PDFs.
In the Appendix we will study correlation function of two light-ray operators [26][27][28][29][30][31][32], in the BFKL limit.Investigating this correlation function in this limit is of great importance, particularly as recent works have primarily focused on the gluon case.Studying the quark correlation function provides a more comprehensive understanding of the underlying physics in the BFKL regime and is crucial for elucidating the behavior of four-point functions in the general case of conformal field theories.The organization of the paper is as follows.In Section 2, we provide a concise introduction to the high-energy OPE formalism and demonstrate its application to the Ioffe-time distribution.In Section 3, we examine the Ioffe-time distribution in the context of leading and next-to-leading twist approximations.In order to reinforce our conclusions pertaining to the high-energy behavior of the Ioffe-time amplitude reached using the Golec-Biernat Wusthoff model, in section 4, our focus will shift to the photon impact factor model, which distinctly, does not embody saturation mechanisms.Sections 5 and 6 are dedicated to the investigation of pseudo and quasi PDFs, respectively, within the BFKL limit and in the leading and next-to-leading twist approximations.Finally, in Section 7, we present our conclusions and discuss the implications of our findings.

Ioffe-time distribution in the BFKL limit
The quark distribution is defined through the light-cone matrix element [27]: with z 2 = 0, P the hadronic momentum, and with the gauge link defined as (2. 2) The distribution Q(x B ) is defined through the quark distribution, D q , and anti-quark distribution, D q as In order to calculate the low-x B behavior of the Deep Inelastic Scattering (DIS) structure function, one can employ the high-energy Operator Product Expansion (OPE) framework.In this approach, the T-product of two electromagnetic currents, TJ µ (x)J ν (y), is expanded in terms of coefficient functions (impact factors), convoluted with the matrix elements of infinite Wilson lines, which are the relevant operators in the high-energy (Regge) limit.The evolution equations governing the matrix elements of infinite Wilson lines with respect to the factorization parameter, or rapidity, are the BK-JIMWLK evolution equations [21][22][23][24][25].This equation accounts for both leading-log resummation and preservation of the unitarity property of the theory.The BK equation is considered a generalization of the BFKL equation, which, although it handles leading-log resummation, results in a violation of unitarity.
The high-energy OPE is analogous to the OPE in non-local (finite gauge-link) operators employed in the Bjorken limit.In this case, the factorization scale is in transverse momentum, and the evolution equation of the finite gauge-link with respect to the factorization parameter is the DGLAP evolution equation.
To examine the high-energy behavior of the Ioffe-time distribution and, consequently, the low-x B behavior of pseudo-and quasi-PDFs, we consider Lorentz decomposition of the structure function defined through the bi-local operator [9]: (2.4) Here, ̺ = z • P , and P represents the hadronic momentum.In the high-energy limit, the 0 component and the 3 component are indistinguishable.Let us introduce the light-cone vectors n and n ′ such that n 2 = n ′2 = 0, n • n ′ = 1, and In this way, a generic vector can be written as where the transverse coordinate is x µ ⊥ = (0, x 1 , x 2 , 0).We adopt the high-energy limit in which the largest light-cone component of the target's momentum is P − .Under this large boost, we have ψ(x In the high-energy limit Consequently, the component of the structure function (2.4) that survives after the boost is M(̺, z 2 ) which contains the leading twist as well as higher twist.On the other hand, the component N (̺, z 2 ) contains only higher twists which are suppressed in the high-energy limit.
At high-energy, the bi-local operator is written in terms of integration along the longitudinal direction as where, as discussed above, at high-energy we can write / z 2z•P → γ − 2P − .We can now apply the high-energy OPE to (2.6).This method was also used in [17] for the study of the high-energy behavior of the gluon Ioffe-time distribution.
To begin calculating the impact factor diagram, we consider the operator in equation (2.6) in the external gluon field, which, at high-energy, shrinks to a shock-wave (represented by the red vertical band in the figure).In the high-energy (Regge) limit, the main degree of freedom are gluons, so it is natural to assume that the field of the target state is predominantly made of gluons .For this purpose, we need the quark propagator in the shock-wave in coordinate space [21] ψ(x) ψ(y) (2.7) where Functionally integrating over the quantum field, which are the quark fields at point x and point y, the Ioffe-time non-local operator, eq.(2.6), becomes The subscript A on the angle bracket means that the operator is being evaluated in the background of the gluon field.In eq.(2.8) we have reduced (expanded) the operator on the left-hand-side (LHS) of the equal sign into a convolution between a coefficient (the impact factor) and new operators given by the trace in the fundamental representation of two infinite Wilson lines, one at the transverse point z 1⊥ and another one at the transverse point z 2⊥ .Let us define X 1⊥ = x ⊥ − z 1⊥ and Y 1⊥ = y ⊥ − z 1⊥ , and from eq. (2.8) we arrive at where we also used tr{γ Note that we indicated with z 2 the point at which the shock wave cuts the quark propagator.The point z 1 , instead, is the point at which the shock-wave cut the straight dotted blue line which represents the gauge-link, and it is a point that runs from point x to point y.
The evolution of the trace of two infinite Wilson lines with respect to the rapidity is the BK equation.However, for our purpose, it will be enough the linearization of the BK eq.i.e. the BFKL equation which we can write as where U (z ⊥ ) is the forward matrix element which depends only on its transverse size and it is obtained from In eq.(2.10) we made use of the two dimensional scalar product defined as (x, y) ⊥ = x 1 y 1 + x 2 y 2 .
Evolution equation (2.10) is with respect to the coordinate space evolution parameter defined as which is reminiscent of the composite Wilson lines operators introduced to restore the Möbius conformal invariance lost at NLO level (for details see Refs.[33][34][35][36]).The peculiarity of the evolution parameter a is that it is in coordinate space and it suits well our purpose.
It is convenient to define the parameters u = |y + | ∆ + , ū = x + ∆ + , and use them to rewrite the solution of the BFKL equation, eq.(2.13) as U a (z 12 ), and ∆ = (x − y).Let us convolute the solution of the evolution equation of forward matrix element, eq.(2.14), with eq.(2.9) and we arrive at The GBW notation within the angle bracket signifies that the matrix element will be evaluated using the GBW model [37].For an alternative model with saturation we can use the McLerran-Venugopalan (MV) [38] model.In the next section, instead, we will consider a model without saturation mechanism.We need the projection of the impact factor on the power like eigenfunctions which is Result (2.16) is obtained using , where z 1⊥ is a point which runs from x ⊥ to y ⊥ , so it can be parametrized as z 1⊥ = x u = ux ⊥ + ūy ⊥ .
In (2.16) where we used The α s correction omitted in equation (2.18) would contribute if we were to also compute the next-to-leading order (NLO) impact factor and employ the solution of the NLO BFKL equation, i.e., perform a comprehensive NLO calculation.
As mentioned above, we will use the GBW model to evaluate the forward dipole matrix element where Q s is the saturation scale and σ 0 = 29.12mb is the dimension of the dipole cross section whose numerical value was obtained from fitting Hera data [37].In high-energy QCD, the saturation scale Q s is particularly important when studying the behavior of hadronic matter at small Bjorken-x.The saturation scale is related to the gluon density inside hadrons and offers a quantitative measure of the transition between the linear regime of parton evolution and the non-linear regime dominated by parton saturation.In our case, we can view the saturation scale as a non-perturbative parameter.Using model (2.19) in eq.(2.17) and integrating over ω ⊥ we arrive at Now, we can finally rewrite result (2.20) using the space like vector z µ , z 2 < 0, and obtain the high-energy behavior of the operator defined in eq.(2.5)As previously mentioned, at high energy, we do not distinguish the 0 component from the 3 component, allowing us to define the Ioffe-time parameter as The evaluation of the integral over the parameter ν in result (2.21) can be performed numerically or in the saddle point approximation.In eq.(2.21), the factor γ3 sin 2 (πγ) is a slowly varying function of ν.Thus, in the saddle point approximation, eq.(2.21) is The outcome presented in eq.(2.22) embodies the characteristic logarithm resummation, featuring the exponentiation of the well-known Pomeron intercept within the saddle point approximation.In this particular instance, the resummed logarithms correspond to the large values of the Ioffe-time parameter ρ.
In Fig. 3, we have plotted both the real and imaginary parts of the quark Ioffe-time distribution using the saddle point approximation, as given by equation (2.22) (represented by the blue dashed curve), and the numerical evaluation of equation (2.21) (depicted by the orange curve).The plots were obtained using . However, for the Ioffe-time amplitude, we will use where ̺ 0 is the starting point of the evolution, which, as can be observed from Fig. 3, is ̺ 0 = 8, so Q s = 0.33 GeV.The value of the strong coupling we used is ᾱs = αsNc π = 0.2.
an integral along the imaginary axes by changing the integration variable from parameter ν, to γ = 1 2 + iν and take its Mellin transform, thus obtaining Performing the Mellin transform we have The integration over γ can be calculated by taking the residue closing the contour to the right because 0 < Where, we recall, ω = j − 1.If we take the inverse Mellin transform of (3.3), we would reproduce again the result of the Ioffe-time amplitude in the BFKL approximation that we already obtained in the previous section.
If we relax the BFKL condition αs ω ≃ 1, and consider the limit α s ≪ ω ≪ 1, we can demonstrate that the result (3.3) can be computed by resumming an infinite series of residues, with each contribution having a higher power of ∆.Since we consider 0 < ∆ 2 ⊥ < 1, the series consists of contributions that become less significant as the power of ∆ 2 ⊥ increases.This is essentially a twist expansion in coordinate space.Our objective is then to calculate the first two twist contributions, which are the ones most likely to be computed within the lattice formalism.
To proceed let us observe that In the limit γ → 1 we have ℵ(γ) → ᾱs 1−γ and 1 where we used the small α s expansion for The next-to-leading residue, i.e. the next-to-leading contribution in terms of Indeed, as done for the leading residue, in the limit of ᾱs ≪ ω ≪ 1, it is easy to show that we just need to make the substitution ℵ(γ) → ᾱs 2−γ and 1 . So from eq. ( 3.2) we take the residue at γ = 2 − ᾱs ω and obtain where expanding in ᾱs we used Besides the series of the moving (dynamical) poles of which we calculated the first two leading ones, there is also a non moving pole at γ = 1.The contribution of this pole does not contribute as explain in the Appendix and as also demonstrated in the gluon case [17].
Summing the two leading residues, results (3.5) and (3.7), we finally arrive at The two leading residues we just calculated organize themselves as an expansion in x 2 ⊥ : the sub-leading term is suppressed by an extra power of x 2 ⊥ which is typical of the coordinate space twist expansion.
We now need to perform the inverse Mellin transform.We will consider only the case 0 < < 1 which is consistent with twist expansion.Indeed, in performing the inverse Mellin transform, one could also consider the case > 1, but this case is inconsistent with twist expansion because higher powers of x 2 ⊥ cannot be disregarded.The inverse Mellin transform of (3.9) is where I 1 (u) is the modified Bessel function with The result given in eq.(3.10) represents the Ioffe-time amplitude, including next-toleading twist corrections.In Fig. 4, we show the real and imaginary components of the leading (eq.(2.22), green dashed curve) and next-to-leading (eq.(3.9), solid red curve) twist corrections, in comparison with the BFKL resummation result.
For the real part, the leading and next-to-leading twist corrections exhibit a value of approximately −1 at Ioffe-time ̺ = 8 and display a modest decrease for larger values of ̺.Conversely, the BFKL resummed curves possess a value of roughly −2 for ̺ = 8 and demonstrate a swift decline for higher values of ̺.Regarding the imaginary part, although an overlap between the leading and next-to-leading twist corrections and the numerical outcome of the BFKL resummation (blue dashed curve) is observed at ̺ = 12, the leading and next-to-leading twist corrections exhibit a gentle increase, while the BFKL resummed result rises more rapidly.
It is noteworthy that the behavior observed for large values of ̺ in the leading and next-to-leading twist, as illustrated in Fig. 4, is consistent with the lattice calculation results shown in Fig. 1 of Ref. [9] (see also [39]).We can observe a similar consistency also with the results presented in Fig. 2 and 3 of Ref. [40] even though in [40] the pseudo-Ioffe-time distribution has been calculated for a maximum value the Ioffe-time parameter ̺ = 8, thus making the comparison (especially for the imaginary part) for large Ioffe-time values a bit difficult.
This consistency suggests that lattice calculations may not adequately capture higher twist contributions.Indeed, the behavior of the pseudo Ioffe-time distribution with BFKL resummation (all twist resummed) represented by the orange solid curve and dashed blue curve in Fig. 4 suggests a stepper decrease for the real part and a steeper increase for the imaginary part.

Ioffe-time amplitude with the photon impact factor model
In this section, our main goal lies in demonstrating that the twist expansion, as derived in the prior section, is not simply an incidental result of the chosen model.To maintain the focus on the salient findings, we confine the intricate details of the computation to the Appendix, while here, we concentrate on the presentation of the primary outcomes.
The Ioffe-time amplitude, when considered in relation to the photon impact factor, is as follows (please see Fig. 5).
Applying the high-energy OPE to (4.1) as done in the case of the GBW model in the previous section, we arrive at (see section A for the details of the calculation) where we remind that ∆ 2 ⊥ = (x − y) 2 ⊥ , and the saddle-point approximation and arrive at In the case of photon impact factor the role of the Ioffe-time parameter is ̺ = Lq − and the large logarithms resummed are of the type ln ̺ for large ̺.
To obtain the leading and next-to-leading twist correction from (4.2) we perform first a Mellin transform, calculate the first two residues and then Mellin transform back and obtain (see section A for the details of the calculation) where I 2 (h) is the modified Bessel function with Result (4.4) shows that the twist corrections organize themselves as an expansion in Q 2 ∆ 2 ⊥ like the one we obtained using the GBW model, thus proving that this feature is not peculiar to the model GBW model only.Moreover, from Fig. 6, where we plot only the leading residue (the next-to-leading residue is just a very small shift), we observe that the behavior of the Ioffe-time amplitude at leading twist (and next-to-leading) with the photon impact factor model is not only similar to the one obtained with the GBW model which is presented in Fig. 4 but also quite close to the one obtained through Lattice calculation in ref. [9].

Quark pseudo-PDF
In this section, we will perform the Fourier transform of the Ioffe-time distribution in the BFKL limit to derive the quark pseudo PDF.Following this, we will calculate the same transformation for both leading and next-to-leading twist corrections.Finally, we will compare the pseudo-PDF in the BFKL approximation with the leading and next-toleading twists by plotting their corresponding values within the x B range of [0.01, 0.1].Here we plot the real and imaginary parts of the numerical evaluation of eq.(5.1) (orange curve) with its saddle point approximation, eq. ( 5.3) (blue curve).

Pseudo PDF in the BFKL approximation
We will perform the Fourier transform of eq.(2.21) following the pseudo PDF definition.
The integration over the parameter ν can be evaluated in saddle point approximation.Using again ℵ(γ = 1/2) = ᾱs 4 ln 2 and we have where we can use again Γ[1 + ᾱs 4 ln 2] sin 2π ᾱs ln 2 ≃ 2π ᾱs ln 2 + O(ᾱ s ).We observe that the behavior of the pseudo-PDF in the BFKL limit is primarily dictated by the well-known exponentiation of the Pomeron intercept, which effectively resums the logarithms of 1/x B .This point is described in Fig. 7, in which we present the results from eq. ( 5.3) along with the numerical evaluation of (5.1).The comparison demonstrates the efficacy of our saddle point approximation as described in eq. ( 5.3).

Pseudo PDF in the leading and next-to-leading approximation
Let us proceed by calculating the pseudo-PDF in the leading and next-to-leading approximation, that is we perform the Fourier transform of the leading and next-to-leading twist correction calculated for the Ioffe-time distribution.
It is convenient to perform the inverse Mellin transform after the Fourier transform.So, our starting point is eq.(3.9) (here L > 0) Since we are in the approximation ω ≪ 1 we can use Γ(1 + ω) ∼ 1 in (5.4) and we obtain where I 1 is the modified Bessel function and we defined . (5.6) We again performed the inverse Mellin transform only in the region 0 < < 1 where the higher twist effects can be considered small corrections.in Fig. 9, we compare result (5.4) with eq.(5.5),where the approximation Γ(1 + ω) ≃ 1 has been employed.Our analysis demonstrates the validity of the approximation Γ(1 + ω) ≃ 1.
In Fig. 8, we observe that the first two twist contributions, which have only the real part, as described by eq.(5.5), exhibit strong agreement with the all-twists resummed BFKL results (5.3) and (5.1) for larger values of x B .However, as we approach smaller x B values, the two twist corrections deviate from the full BFKL result, suggesting that, lattice calculations, now available only at large values of x B , may describe well only the region where leading twist contributions dominate, which is the domain of DGLAP dynamics.On the other hand, the small-x B region, where all twist corrections contribute equally and are described by the BFKL dynamics, is not attainable by lattice calculations.

Quark quasi-PDF in the BFKL approximation
In this section we are going to obtain the quark quasi-PDF by performing the Fourier transform of the Ioffe-time distribution given in eq.(2.21).To this end, following Ref.[17], we introduce the real parameter ς, with −z 2 = ς 2 > 0, and the four-vector Under the large boost we can then identify (LP − ) 2 = (z µ P µ ) 2 = ς 2 P 2 ξ .So, we can perform the following substitution Therefore, starting from (2.5) and (2.21), and using / z 2z•p = 1 2P ξ / ξ, the quasi-distribution, i.e. the pseudo Ioffe-time distribution rewritten á la quasi-PDF is In Fig. 10 we plot the quasi-Distribution, eq.(6.2), for P ξ = 10 GeV (red dashed curve), P ξ = 4 GeV (Blue dashed curve), and P ξ = 2 GeV (orange dashed curve).Result presented in Fig. 10 can be compared with the ones obtained, although with smaller values of P ξ , in lattice calculations in refs.[41,42].
Evaluating the integration over the ν parameter in the saddle point approximation the quasi-distribution (6.2) becomes 3) The quark quasi-PDF is obtained performing the Fourier transform with respect to the ς parameter.We notice that, contrary to the pseudo-PDF case, if we perform the Fourier transform of eq.(6.2) first, and then the integration with respect to the ν parameter, we end up with a divergent integral.So, in this case we will proceed as follow.We first integrate the ν parameter in the saddle point approximation, and then perform the Fourier 0.00 In the left and the right panel we plot, respectively, the real and the imaginary part of the numerical evaluation of eq. ( 6.4) with P ξ = 10 GeV (red curve), P ξ = 4 GeV (blue curve), P ξ = 2 GeV (orange curve).For the numerical evaluation we used ᾱs = 0.2, M N = 1 GeV.and Q s = 0.33 GeV.
transform.So, the numerical evaluation of the Fourier transform of eq. ( 6.3) is In Fig. 11, we exhibit the numerical evaluation of eq.(6.4) with P ξ = 10 GeV (red curve), P ξ = 4 GeV (blue curve) and P ξ = 2 GeV 1 (orange curve), demonstrating the distinctive behavior of the quark quasi-PDF.It is crucial to emphasize that this behavior significantly deviates from the quark pseudo-PDF, as represented in Fig. 8.This divergence can be ascribed to the lack of the typical logarithm resummation that is a hallmark of the BFKL formalism.We should also enphasize that the saddle point approximation is a valid approximation for large values of P ξ .

Quasi PDF in the leading and next-to-leading approximation
To obtain the quasi-PDF in the leading and next-to-leading approximation, we follow a similar approach as employed for the pseudo-PDF.Thus, our starting point is eq.(3.9), 1 The allowed values of P ξ are the ones such that ᾱs ln The value of P ξ = 2 GeV, although barely satisfies this condition, it is used because it is a value attainable in lattice calculations.One should also consider that plots in Fig. 11 are obtained using saddle point approximation which works better at large values of P ξ .Here we present the plots for the LT (magenta curve) and NLT (green dashed curve), eq.(6.5), and compare them with the BFKL result (blue dashed curve), eq. ( 6.2).In the left panel we plot the real parts, while in the right panel we plot the imaginary parts.The plots are obtained using P ξ = 4 GeV, ᾱs = 0.2, M N = 1 GeV.and Q s = 0.33 GeV.
which we rewrite using the notation specific to the quasi-PDF as we did in the previous subsection.The result (3.9) must undergo an inverse Mellin transformation and a Fourier transformation in the quasi-PDF definition, i.e., maintaining the orientation of the z µ vector fixed.First, let us obtain the quasi-distribution in the LT and NLT as inverse Mellin transform of eq.(3.9) which is (recall that ω = j − 1): The Inverse Mellin transform (6.5) has been performed only in the region 0 < Q 2 s ς 2 < 1 which is consistent with higher twist expansion.
In Fig. 12 we plot the real and imaginary parts of the LT and NLT of eq.(6.5), and compare them with the BFKL result eq.(6.2).We notice that the first two leading twist contributions (the green and magenta curves which are one on top of the other), are consistent with the BFKL result (resummation of large logarithms of P ξ ) within the region ς ∈ [0.01, 1].
The quark quasi-PDF in the LT and NLT is obtained performing the Fourier transform of (6.5) with respect to the ς parameter.So we have We should recall that now we are working in the α s ≪ ω ≪ 1 approximation, which means that we are slightly stepping away from the BFKL limit to enter in the DGLAP one.To perform the Inverse Mellin we have to distinguish two ranges of the values x B .For x B < Qs M N √ 2 we have In the left panel we plot the real part of the LT (dark green dashed curve) and NLT (magenta dashed curve) of the quasi-PDF distribution eq.( 6.8) and BFKL resummation (blue curve) eq. ( 6.4).In the right panel we plot the imaginary parts.The plots are obtained using P ξ = 4 GeV, ᾱs = 0.2, M N = 1 GeV, and Q s = 0.33 GeV.
In Fig. 13, we present the leading and next-to-leading twist contributions, as described by eqs.(6.8) and (6.10), for the quasi quark PDF as well as the BFKL resummed result eq. ( 6.8).The integral in eqs.(6.8) and (6.10) is performed by closing the contour to the left.To this end, we have to impose ln The plots for the LT and NLT corrections in Fig. 13 are obtained using P ξ = 4 GeV.A different choice, such as P ξ = 2 GeV, would lead to a divergent integral.This means that the formalism we adopted to study the small-x B behavior of the twist corrections imposes a limit on the allowed values of P ξ .The behavior of the quasi-PDF at LT and NLT is particularly noteworthy due to two primary factors.Firstly, it exhibits a stark deviation from the pseudo-PDF case, and secondly, the next-to-leading twist correction displays a divergent trend that opposes the direction of the leading one at small-x B .This unusual behavior is attributable to the peculiar nature of the higher twist corrections, which are characterized by the 1/(x 2 B P 2 ξ ) factor.Indeed, such higher twist corrections, instead of being suppressed as was hoped, are enhanced at small-x B .We also notice that the full BFKL result is consistent with the LT correction.

Conclusions
By employing the high-energy operator product expansion, we have derived the highenergy behavior of the Ioffe-time distribution for the quark non-local operator.This finding bears considerable importance for lattice calculations, which are not optimally suited for computing large Ioffe-time behavior.In order to obtain the pseudo distribution from lattice calculations, an extrapolation for extensive Ioffe-time values is required prior to conducting the Fourier transform.
In Fig. 4, we presented the Ioffe-time behavior, demonstrating that the Ioffe-time amplitude at leading and next-to-leading twist corrections exhibits a large Ioffe-time parameter behavior that aligns with the results obtained through lattice calculations in Ref. [9].Nonetheless, our analysis also indicates that the leading and next-to-leading twist corrections are insufficient to accurately characterize the behavior determined from the complete BFKL result.Indeed, the higher twist corrections become more important at larger Ioffe-time parameters.While lattice calculations are, in principle, capable of capturing all twist effects, their restriction to lower values of Ioffe-time parameters prohibit them from adequately capturing the domain where the higher twist effects bear the same weight as the leading ones.Consequently, lattice calculations limited to low Ioffe-time parameters likely remain valid only for the first two twist effects which align well with the dynamics as described by DGLAP.To strengthen our findings, we evaluated a model without saturation effects in Section 4 (see Fig. 6).The assessment confirmed that both the pattern of the twist expansion and the behavior of the first two leading twists at high Ioffe-time parameters are not exclusive to the GBW model.
Subsequently, we performed an explicit Fourier transform for both the pseudo-PDF and quasi-PDF from the Ioffe-time amplitude.The pseudo-PDF, demonstrated in Fig. 8, follows the expected result in the saddle point approximation, capturing the BFKL resummation (resummation of all twists) and exhibiting an expected rising behavior for small x B values.Conversely, the quasi-PDF resulted in a different and unusual behavior, reaffirming that they are unsuited for small-x B studies as was previously observed in the gluon case [17].
We also examined the behavior of the pseudo-PDF and quasi-PDF within the first two leading twist contributions, which are potentially attainable from lattice calculations.Although the pseudo-PDF's next-to-leading twist contribution exhibits similar behavior to that with BFKL resummation, it does not approach the BFKL resummation result significantly faster than the leading contribution.This observation confirms that an appropriate small-x B behavior cannot be achieved solely through DGLAP dynamics; rather, an all-twist (BFKL) resummation is required.Indeed, although lattice calculations may in principle capture all twist contributions, the inability to study them in the small-x B region, where the higher twist are as important as the leading one, make lattice calculations relevant only for the leading twist effect which are described by DGLAP dynamics.Furthermore, while the equivalence between the pseudo-PDF and quasi-PDF formalisms can be demonstrated at moderate x B values, the situation differs for small-x B .Our findings show that the higher twist contributions in the quasi-PDF case are not suppressed but instead enhanced, scaling as O(1/(x 2 B P 2 ξ )), as shown in equations (6.8) and (6.10).Consequently, this equivalence between the two formalisms is not expendable to include small x B values.We also studied the quasi-Distribution (see Fig. 12), and we established that at LT and NLT it is in good agreement with the full BFKL resummed result.A similar conclusion can be drawn for the quasi-PDF displayed in Fig. 13.Indeed, although the NLT correction diverges from the LT one for reasons previously discussed, the behavior of the LT, instead, is in good agreement with the full BFKL resummation.impact factor, bottom part of Fig. 5 we use result of section B and arrive at where in this case γ ′ = 1 2 + iν ′ , and Ũ (ν) is defined in eq.(C.32).Notice that in the photon IF case, the initial point of the evolution, a 0 of the evolution parameter (2.11), of the BFKL equation (2.10) has changed.Indeed we have where ∆E = Q 2 √ 2q + is the inverse of the life-time of the quark-anti-quark pair which fluctuates from the virtual photon.The projection of the dipole operator V(ω ⊥ ) onto to the LO eigenfunctions is So, eq.(A.3) becomes We need the dipole-dipole amplitude at LO which in the Mellin space is So, using (A.5) and the LO dipole-dipole scattering (A.7)we arrive at At this point we may proceed similarly to the case GBW model.We take the Mellin transform, calculate the first two leading residue and then Mellin transform back.Starting with the Mellin transform of (A.8) we have In the limit α s ≪ ω ≪ 1 use ℵ(γ) → ᾱs 1−γ and γ → 1 − ᾱs ω , the leading residue is with ω = j − 1 and where for γ → 1 − ᾱs ω we used The next-to-leading residue is obtain in the same way, we use ℵ(γ) → ᾱs 2−γ and calculate the residue closing the contour to the right.Summing the first two residue and performing the inverse Mellin transform for the relevant case ln with h defined in (4.5) and ̺ = Lq − .

B Photon Impact factor
Using result of ref. [36], the T product of two electromagnetic currents can be expanded in terms of the impact factor and matrix elements of the trace of two Wilson lines.
2 s where we defined −q 2 = Q 2 > 0. Further, we define The symbol tilde indicates that the Wilson lines run along x − direction as In eq.(B.1) we defined the two tensor structures P µν 1 and P µν 2 as

C Correlation function with light-ray operators
Correlation functions in Conformal Field Theory play an essential role: they encode information about the spectrum of the operators and their scaling dimensions.In twodimensional CFTs, the primary fields and their descendants can be organized into representations of the Virasoro algebra, and the correlation functions of these primary fields can be determined by conformal symmetry and the operator product expansion (OPE) [43].
In higher-dimensional CFTs, correlation functions can be analyzed using the conformal bootstrap method, which is based on the OPE and the crossing symmetry of the correlation functions.This approach has led significant progress in determining the scaling dimensions and OPE coefficients of higher-dimensional CFTs [44].
In CFTs, local operators, the observables of the theory, are classified according to their scaling dimensions and their transformation properties under the conformal group, and can be organized into representations of the conformal group called conformal multiplets.
Correlation functions, which are the expectation values of products of local operators at distinct spacetime points, are important for the understanding of the dynamics and symmetries of the theory.In CFTs the conformal invariance imposes strong constraints on the form of the correlation functions, leading to a set of algebraic relations known as the conformal bootstrap equations [45].
Conformal invariance completely fixes the form of the two-and three-point functions up to constants, known as the operator product expansion (OPE) coefficients.These coefficients play an essential role in the structure of the theory and can be calculated using various techniques, such as the conformal bootstrap method and integrability [44].
N =4 super Yang-Mills theory is a four-dimensional quantum field theory with maximal supersymmetry, described by a gauge field, four Weyl fermions, and six scalar fields, all transforming in the adjoint representation of the gauge group [46].The N =4 SYM Lagrangian is invariant under supersymmetry transformations, which relate the bosonic and fermionic degrees of freedom.This invariance leads to a rich structure of non-renormalization theorems and exact results that make N =4 SYM particularly tractable [47].
Similar to CFTs, correlation functions in N=4 SYM are central to the understanding the dynamics of the theory.However, the presence of supersymmetry imposes additional constraints on the correlation functions, simplifying their calculation.
In this section we are going to consider the corelation functions of the supermultiplet of twist-two conformal operators, which is a fundamental concept in N =4 super-Yang-Mills theory.A supermultiplet is a collection of states that transform under the same irreducible representation of the superconformal algebra.Twist-two operators are a class of operators with specific scaling dimensions and transformation properties under conformal transformations.They are important because they allow us to study the renormalization properties of the theory by revealing the structure of the theory under conformal transformations.
The twist-two supermultiplets can be classified according to their collinear twist.The term "collinear twist" refers to the difference between the scaling dimension and the collinear spin of the operator, which determines how the operator behaves under conformal transformations.Operators sharing the same collinear twist form a multiplet with the same properties.For example, they have the same scaling behavior under renormalization and similar transformation properties under conformal transformations.
In the context of the superconformal algebra, the irreducible representations are the fundamental building blocks that describe how operators transform under the action of the superconformal algebra.Each irreducible representation is characterized by a set of quantum numbers, such as scaling dimension, spin, and R-charge, which determine the transformation properties of the operators in that representation.The classification of twist-two multiplets by collinear twist is essential for identifying the irreducible representations of the superconformal algebra.When studying the conformal structure of N =4 super-Yang-Mills theory, it is important to understand how the operators transform under the action of the superconformal algebra.By classifying the twist-two multiplets according to their collinear twist, it becomes possible to identify the corresponding irreducible representations and to analyze their transformation properties in a systematic way [48].

C.1 Twist-two non-local operator with non-integer spin j
In this section it is convenient to introduce the light-cone vectors p µ 1 = s/2 n µ and p µ 2 = s/2 n ′µ such that p 1 •p 2 = s/2.We will also use the notation x p 1 = p µ 1 x µ = s/2 x − , and x p 2 = p µ 2 x µ = s/2 x + .To study the correlation function of operators with spin j and j ′ in the high-energy limit, i.e. in the limit in which the contributions αs j−1 ∼ 1 are dominant and needs to be resummed with BFKL, we need to consider the analytic continuation of local operators to non-integers values of spin j.In supersymmetric N =4 gauge theory one can construct the super-multiplet of non-local operators as an analytic continuation of the local ones to non-integer values of spin j [31,48] (see also [17]) The corresponding multiplicatively renormalizable light-ray operators with non-integer j for forward matrix elements are As demonstrated in ref. [28], the analytic continuation of anomalous dimensions of local operator to non integer j gives the anomalous dimension of light-ray operators.So, the anomalous dimensions are In CFT the two-point correlators of light-ray operators is entirely fixed by conformal symmetry up to some unknown structure constant like for the correlators of local operators.So, we may write where in (C.11) ∆ is the dimension of the operator and γ an is the anomalous dimension.
Our aim is to calculate the gluino correlation function (see Fig. 14) in the BFKL limit.However, it is known that correlation functions of non-local operator on the light-cone are divergent in the BFKL limit i.e. in the limit of j → 1 with αs j−1 ∼ 1.For this reason, to regulate the UV divergences one may adopt the point-splitting regulator, thus defying the set of point-splitting super-multiplet of non-local operator with non-integer j as Let us rewrite the operators Λ j and F j in terms of operators S 1 , S 2 , S 3 using (C.7), (C.8), and (C.9) Using (C.20) the correlator with F j operators can be written in terms of the S's operators as , and ∆ ⊥ = x ⊥ − y ⊥ and the same for ∆ ′ ⊥ .Comparing (C.21) with (C.11) we can appreciate how the point-splitting act as an UV regulator.Indeed, the UV regulator µ in (C.21) is now given by |∆ ⊥ ∆ ′ ⊥ |: if the transverse distances go to zero we get UV divergence.
We now observe that, taking the limit j → 1, which is the high-energy limit as we discussed above, the correlation function with F j operators is given only in terms of the S 1 operator so we have Now, le us consider the correlation function for gluino operator (C.19) Therefore, the correlation function with gluinos we have to consider is and study it at high-energy (Regge limit) x + , x ′+ → ∞ and y − , y ′− → −∞ keeping all other components fixed.In this limit the correlation function (C.24) factorizes as [34] (see Fig. 14) Before proceeding with the calculation, we rewrite the operator Λ j p 1 as Now, each factor of (C.25) we can apply the HE-OP which is diagrammatically shown for further details on the dipole-dipole scattering calculation) where, as done above, we defined X ⊥ = x ⊥ +y ⊥ 2 and the same for X ′ ⊥ and we used with resummation parameter the rapidities [34] and the same result with replacement ∆ ⊥ → ∆ ′ ⊥ for the integration over x ′+ and y ′+ for the impact factor with operators along the light-cone vector n ′ .In eq.(C.37) we used B(γ) = Γ 2 (γ) Γ(2γ) .We need the correlation function for the of the gluino operators with specific spins j and j ′ , so using eq.(C.2) we have with j = a + iς, j ′ = a + iς ′ , a ∈ ℜ, and a > 1 + ℵ(γ).The next step is to perform the integration over γ.To this end, we recall that we are in the limit of small ∆ 2 ⊥ and ∆ ′2 ⊥ , so we can calculate the integral by taking the residues closing the contour to the right.Therefore, we need to consider the residue at ℵ(γ) − ω = 0 and at γ = 1.
From eq. (C.36), we notice that besides the moving (dynamical) poles, like the one we just calculated, there ia also residue at γ = 1.To calculate its contribution, we start indeed from eq. (C.36), and expand near γ = 1 obtaining Notice that we now have a simple pole and an essential pole both at γ = 1.Calculating the residue, keeping in mind that we are in the α s ≪ ω limit, we have where we used N 2 − 1 ≃ N 2 c at large N c .
The occurrence of such spurious poles would compromise the anticipated result from conformal symmetry, as presented in eq.(C.11).In my previous work [17] and in ref. [31], it was illustrated that this contribution is offset by the diagrams not encompassed in the HE-OPE formalism, thus validating the expected conformal result (C.11).In a similar manner, here, these spurious poles will also be counteracted by contributions from other diagrams not accounted for within the HE-OPE formalism.At present, the mechanism enabling the cancellation of such spurious poles has yet to be discovered.

Figure 1 .
Figure 1.Diagram for the LO impact factor.We indicate in blue the quantum fields and in red the classical background field which, in the leading accuracy, is made by gluons.

Figure 2 .
Figure 2. In this diagrammatic illustration of the high-energy OPE, the target state is depicted by a green blob in the left panel, while the target is evaluated using the GBW model in the right panel.The dashed black line in the left panel signifies the factorization in rapidity, separating the fields into quantum (upper part of the diagram) and classical (lower part of the diagram).In the right panel, the red band represents the classical field of the target, consisting of gluon fields that undergo Lorentz contraction and time dilation, ultimately forming a shock wave.The resummation of the ladder diagrams, which represent the logarithms of the Ioffe-time parameter, is incorporated through the exponentiation of the Pomeron intercept.

Figure 3 .
Figure 3.In the left and right panel we plot the real and imaginary part, respectively, of the Ioffe-time amplitude; we compare the numerical evaluation of eq.(2.21) (orange curve) with its saddle point approximation, eq.(2.22) (blue dashed curve).To obtain the plots we used |z| = 0.5, and M N = 1 GeV.

Figure 4 .
Figure 4.In the left and right panel we plot the real and imaginary part, respectively of the Ioffe-time amplitude; we compare the numerical evaluation of eq.(2.21) (orange curve) with its saddle point approximation, eq.(2.22) with the LT, eq.(3.5) (green dashed curve), and the NLT (3.9) (red solid curve).

Figure 5 .
Figure 5.Here we diagrammatically illustrate the high-energy OPE with the photon impact factor model.The two black dashed lines represent the factorization in rapidity.The diagram at the bottom of the right panel is the diagram for the photon impact factor.

Figure 6 .
Figure 6.In the left and right panel we plot the real and imaginary part, respectively of the Ioffe-time amplitude at the leading twist eq.(4.4); the plot is obtained using the photon-impact factor model with Q = 0.33 GeV which is the same value as Q s at Ioffe-time ̺ = 8.

Figure 7 .
Figure7.Here we plot the real and imaginary parts of the numerical evaluation of eq.(5.1) (orange curve) with its saddle point approximation, eq.(5.3) (blue curve).

ᾱs2 ln 2 Γ 1 ,Figure 8 .
Figure 8.The plot presents the quark pseudo PDF by comparing the numerical evaluation of eq.(5.1) (illustrated by the orange curve) with its saddle point approximation derived from eq. (5.3) (portrayed by the blue curve).Furthermore, we display the LT (marked by the green dashed curve) and the NLT (signified by the solid red curve) obtained from eq. (5.4).

Figure 9 .
Figure 9.The plot illustrates the quark pseudo PDF for both LT and NLT, considering the cases with and without the approximation Γ(1 + ω) ≃ 1 for ω ≪ 1.The blue and orange curves represent the results obtained without employing the approximation Γ(1 + ω) ≃ 1.

. 1 )Figure 10 .
Figure 10.The plot illustrates the quark quasi-Distribution, eq.(6.2) for two different values of P ξ .The red dashed curve is obtained with P ξ = 10 GeV, the blue dashed one with P ξ = 4 GeV, and the orange dashed one with P ξ = 2 GeV.For the numerical evaluation we also used ᾱs = 0.2, Q s = 0.33 GeV, and M N = 1 GeV.

Figure 11 .
Figure 11.In the left and the right panel we plot, respectively, the real and the imaginary part of the numerical evaluation of eq.(6.4) with P ξ = 10 GeV (red curve), P ξ = 4 GeV (blue curve), P ξ = 2 GeV (orange curve).For the numerical evaluation we used ᾱs = 0.2, M N = 1 GeV.and Q s = 0.33 GeV.

Figure 12 .
Figure12.Here we present the plots for the LT (magenta curve) and NLT (green dashed curve), eq.(6.5), and compare them with the BFKL result (blue dashed curve), eq.(6.2).In the left panel we plot the real parts, while in the right panel we plot the imaginary parts.The plots are obtained using P ξ = 4 GeV, ᾱs = 0.2, M N = 1 GeV.and Q s = 0.33 GeV.

Figure 13 .
Figure13.In the left panel we plot the real part of the LT (dark green dashed curve) and NLT (magenta dashed curve) of the quasi-PDF distribution eq.(6.8) and BFKL resummation (blue curve) eq.(6.4).In the right panel we plot the imaginary parts.The plots are obtained using P ξ = 4 GeV, ᾱs = 0.2, M N = 1 GeV, and Q s = 0.33 GeV.

Figure 14 .
Figure 14.Here we present a diagrammatic representation of the four-point fermion correlation function.The left diagram illustrates the four-point function utilizing the BFKL approximation, whereas the right diagram demonstrates the four-point function with the application of the highenergy OPE.