First Glimpse into the Kaon Gluon Parton Distribution Using Lattice QCD

In this work, we present the first results on the gluon parton distribution of the kaon from lattice quantum chromodynamics. We carry out the lattice calculation at pion mass around 310~MeV and two lattice spacings, 0.15 and 0.12~fm, using $2+1+1$-flavor HISQ ensembles generated by MILC Collaboration. The kaon correlators are calculated using clover fermions and momentum-smearing sources with maximum boost momentum around 2~GeV and high statistics (up to 324,000 measurements). We study the dependence of the resulting reduced Ioffe-time pseudo-distributions at multiple boost momenta and lattice spacings. We then extract the kaon gluon distribution function in the $\overline{\text{MS}}$ scheme at $\mu = 2$~GeV, neglecting the mixing between the gluon and singlet-quark sectors. Our results at the smaller lattice spacing are consistent with phenomenological determinations.


I. INTRODUCTION
The study of pseudoscalar meson structure could shed light into the longstanding quantum-chromodynamics (QCD) mystery of the origin of mass. In the Standard Model, the Higgs boson provides most of the masses of all noncomposite fermions. However, the masses of the valence quarks of mesons contribute only a small fraction of the total meson mass. The remainder is thought to arise from the dynamics of the strong interaction. Studying kaon and pion structure can greatly aid our understanding of the emergence of hadronic mass. One commonly studied type of structure is the parton distribution function (PDF), which gives the probability to find quarks and gluons inside a hadron as a function of their momentum fraction x of the mother hadron. Since the gluon field is responsible for almost all of the pion's mass, studying the gluon PDFs of mesons would be very valuable in understanding the structure of mesons. There have been a number of attempts to extract pion gluon PDFs [1][2][3][4], but not much progress has been made on the kaon PDFs. The future experimental studies, for example, U.S.-based Electron-Ion Collider (EIC) [5,6], the Electron-Ion Collider in China (EicC) [7], the Drell-Yan and J/ψ-production experiments COMPASS++/AMBER [8] in Europe will aim at improving our knowledge of both the pion and kaon gluon and quark PDFs. We refer readers to recent reviews [9][10][11] for more details on the current and future prospects for studies of pion and kaon structure.
Lattice QCD provides a method for determining Standard-Model inputs to improve our knowledge of nonperturbative meson gluon structure. However, due to the notorious noise-to-signal issue associated with gluon observables, there have been only a few efforts to determine the first moment of pion gluon PDF [12,13]  with 450 MeV as the lightest pion mass calculated. No kaon moments have yet been calculated. There is little chance that we will be able to use lattice calculations of increasingly higher moments to reconstruct the x-dependence of the meson PDFs. Recently, direct lattice calculations of the gluon PDF of hadrons have been proposed: pseudo-PDF [14] and quasi-PDF [15,16] approaches, following the proposal of Large-Momentum Effective Theory (LaMET) [17][18][19]. Ongoing efforts using other methods, such as hadronic tensor [20], the operator product expansion [21], lattice cross sections [22,23] and the pseudo-PDF approach [24,25], have also had some success in extracting x-dependent hadron structure. Most such PDF calculations, however, are done using the popular quasi-PDF and pseudo-PDF techniques and mostly limited to isovector quark distributions in the nucleon and the valence-quark distribution in the pion and kaons [16, 25-36, 36-50, 50-63, 63-75].
The first exploratory study applying the quasi-PDF approach to gluon PDFs [39] used overlap valence fermions on gauge ensembles with 2+1 flavors of domainwall fermion and spacetime volume 24 3 × 64, a = 0.1105(3) fm, and M sea π = 330 MeV. The gluon operators were calculated for all spacetime lattice sites and high statistics: 207,872 measurements were taken of the two-point functions with valence quarks at the light sea and strange masses (corresponding to pion masses 340 and 678 MeV, respectively). Both pion and nucleon coordinate-space gluon quasi-PDF matrix-element ratios were compared with global fits after Fourier transformation; the lattice results were consistent with global fits within large uncertainties. Further progress was made using the pseudo-PDF method. Fan et al. studied the nucleon gluon PDF using clover fermions on N f = 2 + 1 + 1 MILC lattices, with a single lattice spacing 0.12 fm and two valence pion masses: 310 and 690 MeV [69]. This was the first x-dependent gluon PDF calculated using lattice QCD. Khan et al. used distillation nucleon interpolating operators and extracted the nucleon gluon PDF on a 2+1-flavor 358-MeV pion mass ensemble at lattice spacing 0.094 fm [76]. Both works neglected the mixing of the gluon operator with the quark-singlet sector.

arXiv:2112.03124v2 [hep-lat] 21 Dec 2021
Since then, Fan has continued to pursue gluon PDF calculations of the nucleon and pion using N f = 2 + 1 + 1 MILC lattices with multiple lattice spacings and pion masses. Reference [77] reported the first lattice-QCD pion x-dependent PDFs at two lattice spacings: 0.15 and 0.12 fm and M π ≈ 220, 310 and 690 MeV); the work used the global-fit PDF quark contribution to estimate the effect of the quark mixings in the gluon PDF determination. However, kaon PDF calculations remain limited to the valence-quark distribution [46].
In this work, we present the first lattice-QCD calculation of the x-dependent kaon gluon parton distribution using the pseudo-PDF method. We study the kaon PDFs using two lattice spacings, 0.12 and 0.15 fm, and one pion mass of 310 MeV. The rest of the paper is organized as follows. In Sec. II, we discuss the operator used in this calculation, the numerical setup of lattice simulation, and how the ground-state kaon gluon matrix elements are obtained on the lattice. In Sec. III, we discuss our strategy for extracting the kaon gluon PDF from our lattice calculations. The systematic uncertainties introduced by different steps are studied, and the lattice-spacing and pion-mass dependence are investigated.

II. LATTICE MATRIX ELEMENTS
The lattice gluon matrix elements of the kaon's groundstate |0(P z ) are calculated for several boost momenta P z and Wilson-line displacement lengths z using the unpolarized gluon operator first introduced in Ref. [14], where O(F µν , F αβ ; z) = F µ ν (z)U (z, 0)F α β (0), the Wilson link length is denoted by z, and the field strength F µν is given by where a is the lattice spacing, g 0 is the strong coupling constant, with the plaquette P µ,ν = U µ (x)U ν (x + aμ)U † µ (x+aν)U † ν (x) and P [µ,ν] = P µ,ν −P ν,µ . The operator i =z,t O(F ti , F zi ; z) corresponds to the same matching kernel as in Ref. [14]. However, this operator is not employed, since it vanishes at P z = 0 due to kinematic reasons, which increases the difficulty of obtaining the distributions. Having calculated the aforementioned matrix elements, we investigate their dependence on Ioffe time, ν = zP z .
On the lattice, we employ clover valence fermions on three ensembles with N f = 2 + 1 + 1 highly improved staggered quarks (HISQ) [78], generated by the MILC Collaboration [79,80] at two lattice spacings (a ≈ 0.12 and 0.15 fm) and pion mass of 310 MeV. The clover quark masses were adjusted to reproduce the lightest strange and light sea pseudoscalar meson utilized by the PNDME collaboration [81][82][83][84]. We apply five steps of HYP smearing [85] for the gluon loops to reduce the statistical noise, as elucidated in Ref. [39]. In order to reach higher meson boost momenta, we use Gaussian momentum smearing [86] for the quark fields. Table I shows the lattice spacing a, valence pion mass M val π and kaon mass M val K , lattice size L 3 × T , number of configurations N cfg , and number of total two-point correlator measurements N 2pt meas for the three ensembles.  For a meson Φ, the two-point correlator is given by where P z is the meson momentum in the z-direction, χ Φ =q 1 γ 5 q 2 is the pseudoscalar-meson interpolation operator, and t is the Euclidean time. |A Φ,i | 2 and E Φ,i denote the amplitude and energy for the ground state (i = 0) and the first excited state (i = 1), respectively. We generate the three-point gluon correlators by multiplying the kaon two-point correlators with the gluon operators. The matrix elements of the gluon operators are then extracted by fitting the three-point correlators to the first terms of the energy-eigenstate expansion, where t sep is the source-sink time separation, and t is the gluon-operator insertion time. The amplitudes A Φ,0 , A Φ,1 , and energies E Φ,0 , E Φ,1 are extracted from twostate fits to the two-point correlators. 0|O|0 , 0|O|1 ( 1|O|0 ), and 1|O|1 denote the ground-state matrix element, the ground-excited matrix element, and the excited-state matrix element, respectively. The ground-state matrix element was obtained from a two-state fit to the three-point correlators or from a two-state simultaneous (two-sim) fit across multiple time separations. In order to corroborate that the fitted matrix elements were extracted correctly, the ratio of the three-point to the two-point correlator was computed. The first column of Figs. 1 and 2 show the kaon R ratio for the ensembles a12m310 and a15m310 for selected momenta P z and Wilson-line length z. If there were no excited states, Eq. 5 would yield the groundstate matrix element. It can be seen that R ratio increases as the source-sink separation increases. In some cases, R ratio starts to converge at large separation, which shows that the contribution from excited states diminishes. The gray bands depict the ground-state matrix elements obtained from a "two-sim" fit to the three-point correlators at t sep ∈ [5,9]. In the second and third columns, one can see the variation of the fitted ground-state matrix elements for different choices of t min sep and t max sep .
In the examples shown in Fig. 2, our ground-state fitted matrix elements are stable for the a15m310 ensemble regardless of the choices we made for t min,max sep . However, for the example from the a12m310 ensemble, we found that the ground-state matrix elements are sensitive to the selection of t sep . Our choices of t sep ∈ [5,9] for the analysis do provide a reliable ground-state extraction.

III. KAON GLUON PDFS FROM THE PSEUDO-PDF METHOD
The Ioffe-time pseudo-distribution (pITD) [24,25] is given by The reduced pITD (RpITD) [15,25,87] was constructed to remove the ultraviolet divergences in the pITD by taking the ratio of the pITD to its corresponding zdependent matrix element with P z = 0, then normalizing the ratio by the matrix element at z 2 = 0, The renormalization of O(z) and kinematic factors are cancelled in the RpITDs. By construction, the RpITD double ratios employed here are normalized to one at z = 0. Figure 3 shows the lattice-spacing dependence of the kaon gluon RpITD. The RpITDs are plotted as a function of Wilson link length z at approximately the same momentum P z for ensembles a12m310 and a15m310, which have lattice spacings around 0.12 and 0.15 fm. We observe that the RpITD from the a12m310 ensemble has higher values for each Wilson-link length and has a slower rate in decrease of its values. Most points are within one sigma of their corresponding points on the other ensemble, with slightly larger central values for most z. As the Wilson-link length increases, we see the difference increases.
Using the matrix elements obtained from Eq. 4 and Eq. 7, we compute the RpITD for both ensembles at the various Wilson link lengths z and meson momenta P z . Figure 4 shows these results compared to the RpITD previously calculated for the pion on the a12m310 ensemble in Ref. [77]. We see that the kaon RpITDs for both ensembles, shown in the left and middle plots, are congruent. Additionally, the RpITD for the kaon and pion, shown in the rightmost plot, from a12m310 are very similar, in accordance with the results from Ref. [10], where the ratio of the kaon PDF to the pion PDF was calculated and found to be around 1.
We can then extract the gluon PDF distribution through the pseudo-PDF matching condition [14] that connects the RpITD M to the lightcone gluon PDF g(x, µ 2 ) through where µ is the renormalization scale in the MS scheme and x g = 1 0 dx xg(x, µ 2 ) is the gluon momentum fraction of the kaon. R gg is the gluon-in-gluon matching kernel where α s is the strong coupling at scale µ, N c = 3 is the number of colors, and γ E = 0.5772 is the Euler-Mascheroni constant. For the term R gg (y, z 2 , µ 2 ), z was chosen to be 2e −γ E −1/2 /µ so that the logarithmic term vanishes, which suppresses the residuals containing higher order logarithmic terms, following the previ-a12m310 two-sim . The gray band shown in all plots corresponds to the extracted ground-state matrix element from the two-sim fit using tsep ∈ [5,9]. From left to right, the columns represent: the ratio of the three-point to two-point correlators with the reconstructed fit bands from the two-sim fit using tsep ∈ [5,9], shown as functions of t − tsep/2, the two-sim fit results using tsep ∈ [t min sep , 9] as functions of t min sep , and the two-sim fit results using tsep ∈ [5, t max sep ] as functions of t max sep .
a15m310 two-sim two-sim fit two-sim . The gray band shown in all plots corresponds to the extracted ground-state matrix element from the two-sim fit using tsep ∈ [5,9]. From left to right, the columns represent: the ratio of the three-point to two-point correlators with the reconstructed fit bands from the two-sim fit using tsep ∈ [5,9], shown as functions of t − tsep/2, the two-sim fit results using tsep ∈ [t min sep , 9] as functions of t min sep , and the two-sim fit results using tsep ∈ [5, t max sep ] as functions of t max sep . ous publication regarding the one-loop evolution of the pseudo-PDF [64]. Equation 8 and the terms R B (y), R L (y), R C (y) are defined in Eqs. 7.21-23 and in the paragraph below Eq. 7.23 in Ref. [14].
Note that the lattice-calculated RpITDs are also connected to the singlet quark-PDF q s of the kaon via the quark-gluon matching kernel R gq with an additional x g R gq (xν, z 2 µ 2 ) term added to Eq. 8. In the past pion gluon PDF study [77], it was found that the quark PDF contribution is small; therefore in this work, we will neglect the kaon quark PDF. One can obtain the gluon PDF g(x, µ 2 ) by fitting the RpITD through the matching condition in Eq. 8; a similar procedure has also been used by HadStruc Collaboration [63,67,76].
To obtain the gluon PDF g(x, µ 2 ) on the right-hand side of Eq. 8, we adopt the phenomenologically motivated form for x ∈ [0, 1] and zero elsewhere. The beta function B(A + 1, C + 1) = 1 0 dx x A (1 − x) C is used to normalize the area to unity. Such a form is also used in global fits to obtain the pion gluon PDF by JAM [1,2]. We fit the lattice RpITDs M lat (ν, z 2 , a, M π ) obtained in Eq. 7 to the parametrization form M fit (ν, µ, z 2 , a, M π ) in Eq. 8 by minimizing the χ 2 function, The reconstructed fit bands of the kaon RpITDs for the a15m310 and a12m310 ensembles, and pion RpITDs for a12m310 at each z 2 , compared with the lattice calculation points, are shown in Fig. 4 from left to right. We see almost no z 2 -dependence (labeled in different colors) in the reconstructed bands in both a12m310 ensembles, but slightly more dependence in the a15m310 case. The RpITD data points fluctuate more in the a15m310 kaon data and as a result the fit quality is not as good as in the a12m310 case. We found the a12m310 fit to both the pion and kaon gluon PDF to have very stable quality with χ 2 /dof around 1 with consistent output of f g (x, µ), regardless of the choice of the maximum value of the Wilson-line displacement z. However, for the a15m310 ensemble, χ 2 /dof can go as large as 9(3). By reducing the fitted z to smaller number z = 1, χ 2 /dof can be reduced to 5(3). We suspect that higher-twist effects are enhanced at this coarse lattice spacing such that the fit fails to accurately describe the lattice data. Possible future work including NNLO matching may help to improve the fit on this ensemble.
Our results on the kaon and pion gluon PDFs xg(x, µ)/ x g as a function of x from fitting the lattice data are shown in Fig. 5. In the inset of the left-hand plot, we show the kaon xg(x, µ)/ x g at µ = 2 GeV in the MS scheme with M π ≈ 310 MeV from two lattice spacings a ≈ {0.12, 0.15} fm. Note that we use z = 1 for the coarser lattice spacing due to the fit quality becoming significantly worse as z increases, whereas we are able to get a good fit for the finer-lattice RpITD data for z up to 5. With these choices of data input, our best determination of the kaon xg(x, µ)/ x g are consistent within the statistical errors. We also compare our finer-lattice result with the same quantity obtained from the Dyson-Schwinger equation (DSE) approach [88], depicted as the purple band on the left-hand side. Reference [88] provides recent DSE calculations of the pion gluon PDF and the ratio of kaon to pion gluon PDFs as function of x. Although the kaon gluon PDF is not directly provided, it can be obtained from the product of the pion gluon PDF and the kaon-to-pion ratio. Accordingly, we can make a direct comparison between our lattice calculations and the kaon gluon PDF from DSE calculations. On the other hand, the DSE'20 error shown is likely overestimated due to lack of correlations in the published parameters for their kaon and pion PDFs. Our kaon gluon PDF results from the smallest lattice spacing, a ≈ 0.12 fm, is consistent with the DSE results [88] within one-sigma error for x > 0.15. On the right-hand side of Fig. 5, we compare the obtained kaon result at smaller lattice spacing with the pion results obtained from the same ensembles; we note the kaon gluon PDF is slightly smaller than the one obtained for the pion, similar to its corresponding quark valence PDF and DSE [88] results. We predict the kaon x 2 g / x g and x 3 g / x g to be 0.0779(94) and 0.0187 (42), while in good agreement with corresponding results from DSE [88]: 0.075 and 0.015. On the pion x 2 g / x g , our 310-MeV results gives .092 (15) and 0.0250 (75), while results from DSE [88], JAM [1,2] and xFitter [3] are 0.076, 0.103, 0.158 and 0.015, 0.024, 0.048, respectively. Future study including finer lattice spacing and lighter pion mass will be important to refine this calculation and provide better predictions on this poorly known meson quantity.

IV. SUMMARY
In this work, we made the first lattice study of kaon gluon parton distribution, using the pseudo-PDF approach. We used clover fermions as valence action and 310-MeV 2+1+1 HISQ configurations generated by the MILC collaboration at two lattice spacings, 0.15 and 0.12 fm. We used momentum smearing and highstatistics measurements, up to O(324,000), to reach the kaon boost momentum around 2 GeV. We carefully studied the excited-state contributions to the matrix elements using a two-state fitting strategy and made sure that our ground-state matrix elements were stably obtained. We then calculated the reduced pseudo-ITD using the obtained fitted ground-state matrix elements and extracted the gluon parton distribution. We found that the kaon gluon PDF at the finer lattice spacing is consistent with the DSE result [88] within statistical uncertainties, except in the small-x region, which our zP z is too small to constrain. When comparing with the pion PDF result, we found the kaon PDF to be slightly smaller in central value for most of the x > 0.2 region. We found that the kaon gluon PDFs show potential discretization error at the coarse lattice spacing of 0.15 fm; future study using an additional lattice spacing of 0.09 fm would give us a better estimate of the systematics uncertainty in the 0.12 fm results. We suspect the quark-gluon mixing is smaller than our statistical error based on the prior pion gluon calculation. Other systematics from sources such as higher-twist contributions, finite-volume effects, unphysical pion-mass effects are not included in this pioneer study. Future studies should aim at improving these systematics and provide a better determination of the kaon gluon PDFs for the upcoming experimental efforts.
in part by Michigan State University through computational resources provided by the Institute for Cyber-Enabled Research (iCER). The work of AS, ZF and HL are partially supported by the US National Science Foundation under grant PHY 1653405 "CAREER: Constraining Parton Distribution Functions for New-Physics Searches".