Proton Structure Functions at Next-to-Leading Order in the Dipole Picture with Massive Quarks

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Proton Structure Functions at Next-to-Leading Order in the Dipole Picture with Massive Quarks © Authors 2023. Published by the American Physical Society Published version Hänninen, Henri; Mäntysaari, Heikki; Paatelainen, Risto; Penttala, Jani Hänninen, H., Mäntysaari, H., Paatelainen, R., & Penttala, J. (2023). Proton Structure Functions at Next-to-Leading Order in the Dipole Picture with Massive Quarks. Physical Review Letters, 130(19), Article 192301. https://doi.org/10.1103/PhysRevLett.130.192301 2023 Proton Structure Functions at Next-to-Leading Order in the Dipole Picture with Massive Quarks Henri Hänninen,1,2,3,* Heikki Mäntysaari ,2,3,†Risto Paatelainen ,3,4,‡and Jani Penttala 2,3,§ 1Department of Mathematics and Statistics, University of Jyväskylä, P.O. Box 35, 40014 University of Jyväskylä, Finland 2Department of Physics, University of Jyväskylä, P.O. Box 35, 40014 University of Jyväskylä, Finland 3Helsinki Institute of Physics, University of Helsinki, P.O. Box 64, 00014 University of Helsinki, Finland 4Department of Physics, University of Helsinki, P.O. Box 64, 00014 University of Helsinki, Finland (Received 11 November 2022; revised 3 February 2023; accepted 18 April 2023; published 8 May 2023) We predict heavy quark production cross sections in deep inelastic scattering at high energy by applying the color glass condensate effective theory. We demonstrate that, when the calculation is performed consistently at next-to-leading order accuracy with massive quarks, it becomes possible, for the first time in the dipole picture with perturbatively calculated center-of-mass energy evolution, to simultaneously describe both the light and heavy quark production data at small xBj. Furthermore, we show how the heavy quark cross section data provides additional strong constraints on the extracted nonperturbative initial condition for the small-xBj evolution equations. DOI: 10.1103/PhysRevLett.130.192301 Introduction.—Probing the properties of the nonlinearly behaving gluonic matter in protons and nuclei at high energies is a major science goal of the future Electron-Ion Collider (EIC) [1–3]. Measuring the total and heavy quark production cross sections in deep inelastic scattering (DIS) off nuclei is especially intriguing, as nonlinear saturation effects are enhanced in heavy nuclei [4]. The EIC will be able to perform very precise total cross section measurements over a relatively wide kinematical domain characterized by the gluon longitudinal momentum fraction xBj and the photon virtuality Q2. Nonlinear gluon saturation effects are expected to have a modest effect on structure functions in the EIC kinematics (see, e.g., [5–7]). To unambiguously determine the existence of nonlinear QCD dynamics at collider energies and to quantify its role on the small-xBj structure of protons and nuclei, it is likely necessary to perform a global analysis of the future proton and nuclear DIS data at small xBj.In particular, it will be important to include both the inclusive and heavy quark production data that have different sensitivities on saturation effects in order to extract in detail the properties of the QCD matter at extremely large parton densities. Charm production is an especially powerful process as the charm mass is large enough to suppress nonperturbative effects, but simultaneously light enough to allow one to access QCD dynamics in the nonlinear regime. To describe QCD dynamics at high energies, where parton densities are very large and emergent nonlinear phenomena dominate, it is convenient to use the color glass condensate (CGC) [8,9] effective field theory framework. The DIS process is then naturally described in the dipole picture [10,11], where the photon splits into a quarkantiquark pair long before the interaction with the target. The interaction of the quark dipole with the target is then taken to be eikonal, i.e., the transverse coordinates of the partons do not change when they traverse through the target color field. In this picture, leading-order (LO) calculations including a resummation of the high-energy logarithms αsln 1=x to all orders (where αsis the strong coupling) within the CGC framework have been successful in describing the precise proton structure function data from HERA [12–14]. This suggests that the HERA data is compatible with the hypothesis that gluon saturation is manifest at HERA energies. In addition, calculations based on collinear factorization have also found the resummation of the high-energy logarithms to be important in order to describe the details of the HERA data [15]. The structure function data is used to constrain the nonperturbative initial condition for the small-xBj evolution equations. Therefore, a good description of the total cross section data is crucial when applying the CGC framework to describe any other scattering process (e.g., protonnucleus collisions at the LHC [13,16–20]). Compatibility with the available cross section data is also required when developing a realistic description for the early stages of heavy-ion collisions [21], needed to extract the fundamental properties of the quark-gluon plasma. In this Letter, we present predictions for heavy quark production cross sections in DIS using the nonperturbative Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW LETTERS 130, 192301 (2023) 0031-9007=23=130(19)=192301(7) 192301-1 Published by the American Physical Society initial condition for the perturbative Balitsky-Kovchegov (BK) small-xBj evolution equation [22,23], determined from the fits to total DIS cross section data in [24]. The predicted heavy quark cross sections are shown to be in excellent agreement with the HERA data [25]. This is the first time in the CGC framework that a simultaneous description of total and heavy quark production data is achieved in calculations where the energy dependence is obtained by solving the small-xBj evolution equation. A crucial ingredient, here, is the next-to-leading order (NLO) accuracy in αsrecently achieved for the massive impact factors from first-principle light-cone perturbation theory calculations [26–28]. We also demonstrate how the heavy quark production data can provide additional constraints for the extracted nonperturbative initial condition of the BK evolution. The results presented here are from the first-ever numerical calculation of the heavy quark structure functions in the dipole picture at NLO. The successful description of the HERA data demonstrates that future global analyses are feasible and can be applied to probe in detail gluon saturation at the LHC and future EIC, where nuclear targets with larger saturation scales are available. Structure functions at high energy.—Using the optical theorem, the total virtual photon ðγÞ—proton (p) cross section can be obtained from the forward elastic γþp→ γþpscattering amplitude. In the dipole picture, the γþpscattering is described in terms of eikonal interactions between the partonic Fock states of the photon and the target color field, and perturbatively calculable impact factors describing the photon fluctuations to the given partonic states. Eikonal interactions with the target are encoded in the Wilson lines, which are the scattering matrix elements for bare partons propagating through the target color field. At NLO, the contributing photon Fock states are the quark-antiquark jq¯ qiand quark-antiquark-gluon jq¯ qgi states. Therefore, at NLO, the total virtual photon cross section can be schematically decomposed into two parts. The first contribution [illustrated in Fig. 1(a)] corresponds to the case where the q¯ qg system crosses the shockwave σγ q¯ qg ¼Kq¯ qg ⊗N012:ð1Þ The second contribution [illustrated in Fig. 1(b)], which includes the lowest-order part (interaction with an unevolved target) and the one-gluon-loop QCD corrections to it, reads σγ q¯ q¼Kq¯q⊗N01:ð2Þ Here, Kq¯ qand Kq¯ qg refer to the perturbatively computed NLO impact factors obtained with massive quarks in [26– 28] and in the massless quark limit in [29–31]. In addition, the notation ⊗refers to an integral over the parton transverse coordinates xiand longitudinal momentum fractions in the mixed space. Additionally, N01 and N012 are correlators of two or three Wilson lines, where the subscripts 0,1,2 refer to the transverse coordinates of the quark, antiquark, and the gluon. In terms of the Wilson lines, VðxÞ in the fundamental representation, these correlators read S01 ¼1 NchTrfVðx0ÞV†ðx1Þgi;ð3Þ S012 ¼Nc 2CFS02S21 − 1 N2 c S01:ð4Þ Here, hirefers to the average over the target color charge configurations, Ncis the number of colors, CF¼ðN2 c−1Þ= ð2NcÞ,Sij ¼1−Nij, and Sijk ¼1−Nijk. In addition, we have used the mean-field limit (which is a precise approximation [32]) to factorize the expectation value of the product to a product of expectation values. The Wilson lines and their correlators satisfy small-xBj evolution equations describing their dependency on the center-of-mass energy (see Ref. [33] for a detailed discussion of the evolution variable). The dipole amplitude N01 satisfies the BK equation [22,23] and via Eq. (4) N012 also depends on the center-of-mass energy. The evolution rapidity depends on the lower limit of the emitted gluon longitudinal momentum fraction [24,34]. The integration over the emitted gluon phase space in Eq. (1) contributes a large logarithm of energy that modifies the scattering amplitude of the original dipole N01. These logarithms are resummed into the BK equation [34]. The BK equation and a numerical solution to it are known at NLO [35–37]. We use the initial condition fitted to the HERA data in [24] including only massless quarks, where the full (numerically heavy) NLO BK equation has been approximated by evolution equations that use different schemes to resum the most important higher-order corrections. The same evolution equations, ResumBK [38,39], KCBK [40], and TBK [33], referring to different resummation schemes, are used in this Letter as in [24]. The structure functions are written in terms of the total virtual photon-target cross sections as F2¼ðQ2=4π2αemÞ× ðσγ Tþσγ LÞ, and FL¼ðQ2=4π2αemÞσγ L. Here, the subscripts Tand Lrefer to the transverse and longitudinal (a) (b) FIG. 1. Example diagrams contributing to the elastic γp amplitude at NLO. The blue band represents the dipole– shockwave interaction. PHYSICAL REVIEW LETTERS 130, 192301 (2023) 192301-2 virtual photon polarization, respectively, and σγ T;L correspond to a sum of q¯ qand q¯ qg contributions. The experimental data is reported in terms of the reduced cross section σrðy; x; Q2Þ¼F2ðx; Q2Þ−y2 1þð1−yÞ2FLðx; Q2Þ;ð5Þ where y¼Q2=ðsxÞis the inelasticity and ffiffiffi s pis the leptonnucleon center-of-mass energy. Results.—We calculate the proton reduced cross section σrand the charm and bottom contributions to it (σr;c and σr;b). We use the NLO dipole-proton scattering amplitudes determined in [24], available at [41]. In particular, we use the “light quark”fits of [24] where only the massless u,d, and squarks are included and the nonperturbative initial condition is fitted to the light quark contribution of the reduced cross section data measured at HERA [42]. This contribution is determined in [24] by subtracting interpolated charm and bottom quark contributions from the total cross section data. We do not include the fits to the inclusive HERA data as they use the massless quark cross sections to fit the inclusive data containing a substantial heavy quark contribution. In [24], multiple different fits are reported, corresponding to different choices for the initial evolution rapidity Y0;BK and different schemes for the coordinate space running coupling and resummations of particular higherorder corrections. In total, there are 12 fits reported for massless quarks. All different fits result in an approximately equally good description of the light quark contribution to the HERA structure function data. We calculate predictions for the charm production cross section in the region xBj <0.01;2.5GeV2≤Q2<50 GeV2 using all the different fits from [24], and compare the result to the HERA data from [25] in order to find which fits (if any) are allowed by the heavy quark production data. The charm mass (in the pole mass scheme used in the calculation of [27]) is allowed to vary within 1.1GeV <m c< 1.6GeV. We consider a fit to be compatible with the HERA charm production data if one obtains χ2 c=N ≲2.5 with the optimal charm mass. We find that predictions calculated by using three of the 12 fits are in excellent agreement with the charm production data. This is illustrated in Fig. 2, where a comparison to the HERA reduced cross section data in a few selected Q2bins is shown. The H1 and ZEUS collaborations have also measured inclusive bquark production [25], but due to the larger uncertainties and more limited kinematical coverage, we do not use this dataset to determine which NLO fits from [24] are allowed. However, we note that, using each of the three fits discussed above, an excellent description of the bquark production data is obtained. In each case, we find χ2 b=N ≲ 1.6when the bquark mass is also fitted to this data. The excellent agreement with the predicted heavy quark production cross sections and the HERA measurements shows that, at NLO, it is possible to simultaneously describe all small-xBj proton structure function data. The results also demonstrate that the inclusion of the heavy quark production data to the extraction of the nonperturbative initial condition for the high-energy evolution equation provides additional tight constraints. Similar conclusions have also been made in calculations of exclusive heavy quarkonium production [43,44]. The advantage of the charm reduced cross section studied in this Letter is that one does not need to introduce an additional model uncertainty related to the nonperturbative vector meson structure. The fits that are found to be compatible with the charm quark production data are summarized in Table Ialong with the determined optimal heavy quark masses. The fact that the heavy quark data provides additional strong constraints for the determination of the initial condition for the BK evolution is expected. The heavy quark cross section is sensitive to much smaller dipoles than the inclusive one which can not discriminate fits that differ only at small dipole sizes. We note that the heavy quark production data FIG. 2. Charm reduced cross section predictions calculated using the different NLO fits from [24] for the dipole amplitude that result in a good description of the charm data. The results are compared to the combined HERA data from [25]. TABLE I. Fitted initial conditions for the small-xBj evolution at NLO from [24] that are compatible with the heavy quark production data from HERA. The corresponding charm and bottom masses are also shown. The terminology used to specify the resummation scheme and the running coupling prescription follows that of [24], and the abbreviation PD refers to parent dipole and BSD to Balitsky þsmallest dipole [45] running coupling. # Resummation scheme αsY0;BK mc (GeV) χ2 c=N mb (GeV) χ2 b=N χ2 tot=N 1 ResumBK PD 0 1.42 1.86 4.83 1.37 1.25 2 KCBK PD 0 1.49 2.55 4.96 1.58 1.23 3 TBK BSD 0 1.29 1.02 5.04 1.12 1.83 PHYSICAL REVIEW LETTERS 130, 192301 (2023) 192301-3 only allows fits where the BK evolution is started at initial rapidity Y0;BK ¼0. In the second class of fits considered in [24], the dipole is frozen in the low-energy region 0<Y<Y 0;BK ¼lnð1=0.01Þwhere Yis the evolution rapidity. This is not completely consistent as the q¯ qg production cross section (1) in the soft gluon limit results in a (leading order) BK evolution for the dipole. Additionally, we note that (in the case of ResumBK and KCBK evolutions formulated in terms of the projectile rapidity) the parent dipole prescription for the running coupling is preferred. We interpret that these physical constraints from the heavy-quark production data are because charm and bottom production probe dipole amplitudes in the perturbative region, and contribution from large dipoles dominating in light quark production with N∼1, is suppressed, see, e.g., [6]. Based on the observations above, we argue that the fits summarized in Table Iare the ones that should be used in all NLO CGC calculations. The potential deviation between the predictions is then a measure of the model uncertainty after the nonperturbative input is constrained by all HERA structure function data. To more clearly illustrate the compatibility of the NLO CGC calculation with the most recent precise HERA data from [46], we show, in Fig. 3, the total reduced cross sections computed using the dipole amplitude fits allowed by the charm data. We emphasize that this is the first time in the CGC framework that a simultaneous description of both the total and heavy quark production cross section is obtained when a perturbative small-xBj evolution equation is used to describe the center-of-mass energy dependence. Previous LO analyses have found it impossible to perform such a global fit to the HERA data without introducing, for example, additional parameters that render the proton probed by a charm quark dipole different from the proton probed by a light quark dipole [12]. A similar approximative NLO evolution equation as in this Letter was used in [14] but coupled to the LO impact factor. In that case, it was also found impossible to simultaneously describe the inclusive and heavy quark production data. When the computation is promoted to full NLO accuracy, the mass dependence is modified for two reasons. First, after including higher-order corrections to the BK equation (in projectile rapidity), the dipole amplitude no longer evolves toward an asymptotic shape with an anomalous dimension γ<1(at small dipole sizes r, the amplitude behaves as N∼r2γ)[47]. Instead, the anomalous dimension (which is γ≳1in the fits reported in [24]) remains approximatively constant suppressing the dipole amplitudes at small dipoles [24,37]. Hence, the heavy quark production cross section is suppressed relative to light quark production. Second, adding the NLO corrections to the massive impact factor enhances the heavy quark production. With TBK evolution, we have opposite systematics: a small γis developed and the impact factor suppresses heavy quark production. The net effect of these two competing NLO corrections is such that the mass dependence of the cross section matches that of the HERA data when the three fits identified in this Letter are used. Finally, we illustrate the remaining theory uncertainty when performing NLO CGC calculations. We calculate predictions for the proton longitudinal structure function FL, and for the charm and bottom quark contributions to it, in the EIC kinematics. We take xBj ¼2×10−3, and show, in Fig. 4, the structure functions as a function of Q2 calculated using the three fits determined above. For the bottom structure function, the different fits result in almost identical predictions for the EIC, whereas for charm production, the predictions begin to differ at Q2≳20 GeV2. On the other hand, in the total longitudinal cross section a significant difference up to 20% is seen at all Q2. Therefore, an inclusion of the future FLdata in the global analysis will FIG. 3. Total reduced cross section calculated using the dipole amplitude fits allowed by the heavy quark production data. Note that, as the σrdepends on inelasticity y, the theory curves connecting the calculated points are not smooth. FIG. 4. Total (solid lines), charm (dashed lines), and bottom (dotted lines) longitudinal structure functions as a function of photon virtuality in the EIC kinematics calculated using the three dipole amplitude fits compatible with the heavy quark data. PHYSICAL REVIEW LETTERS 130, 192301 (2023) 192301-4 provide further constraints for the initial condition of the small-xBj evolution. The currently available FLdata from HERA [48] is not able to distinguish between the different fits. Discussion.—We have calculated heavy quark production cross sections in DIS at NLO in the CGC framework. The Bjorken-xBj dependence is obtained by solving the BK evolution equation with an initial condition extracted by fitting the total DIS cross section data in [24]. We identify a small subset of the fits reported in [24] that result in predictions for the charm and bottom structure functions which are in excellent agreement with the HERA data [25]. These three fits, constrained by both the total and heavy quark cross section data summarized in Table I, should be used in all future phenomenological applications at NLO accuracy. This is the first time in the CGC framework with perturbative energy evolution when a simultaneous description of all small-xBj proton structure function data is obtained. A good agreement with the HERA measurements is a crucial test for the gluon saturation physics incorporated in the CGC framework and enables rigorous studies of nonlinear QCD dynamics in DIS and other scattering processes. In particular, we demonstrate that global analyses including all small-xBj structure function data are feasible at NLO and that the heavy quark production data can provide additional constraints in such analyses. As an application, we have calculated predictions for the proton longitudinal structure function FL, which will be measured accurately at the future Electron-Ion Collider. We reported predictions separately for the inclusive and heavy quark production cross section and showed that the remaining model uncertainty is moderate. Including the FLdata to the global analysis will further constrain the nonperturbative initial condition for the small-xBj evolution equations. To fully explore the model uncertainties, one should perform a global analysis to the HERA inclusive and heavy quark production data, taking into account the correlated experimental uncertainties, and extract the nonperturbative model parameters with their uncertainties directly from such an analysis. Additional constraints and more detailed probes of nonlinear dynamics can be obtained by including other observables such as diffractive structure functions and exclusive cross sections. Such studies are becoming feasible thanks to the extensive progress toward NLO accuracy in the CGC framework, see, e.g., Refs. [43,44,49–58].In the future, we plan to perform a full Bayesian analysis to determine the likelihood distribution for all the model parameters, which will also enable one to fully take into account the propagation of model uncertainties. We thank T. Lappi for discussions, and V. Apaja for computational support. This work was supported by the Academy of Finland, the Centre of Excellence in Quark Matter and Projects No. 338263, No. 346567, No. 321840, No. 347499, and No. 353772. This work was also supported under the European Union’s Horizon 2020 research and innovation programme by the European Research Council (ERC, Grant Agreement No. ERC-2018-ADG- 835105 YoctoLHC) and by the STRONG-2020 project (Grant Agreement No. 824093). J. P. is supported by the Finnish Cultural Foundation. The content of this article does not reflect the official opinion of the European Union and responsibility for the information and views expressed therein lies entirely with the authors. Computing resources from CSC–IT Center for Science in Espoo, Finland and from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533) were used in this work. *[email protected] †[email protected] ‡[email protected] §[email protected] [1] R. Abdul Khalek et al., Science requirements and detector concepts for the electron-ion collider: EIC Yellow Report, Nucl. Phys. A1026, 122447 (2022). [2] E. C. Aschenauer, S. Fazio, J. H. Lee, H. Mäntysaari, B. S. Page, B. Schenke, T. Ullrich, R. Venugopalan, and P. Zurita, The electron-ion collider: Assessing the energy dependence of key measurements, Rep. Prog. Phys. 82, 024301 (2019). [3] A. Accardi et al., Electron Ion Collider: The next QCD frontier: Understanding the glue that binds us all, Eur. Phys. J. A 52, 268 (2016). [4] H. Kowalski, T. Lappi, and R. Venugopalan, Nuclear Enhancement of Universal Dynamics of High Parton Densities, Phys. Rev. Lett. 100, 022303 (2008). [5] C. Marquet, M. R. Moldes, and P. Zurita, Unveiling saturation effects from nuclear structure function measurements at the EIC, Phys. Lett. B 772, 607 (2017). [6] H. Mäntysaari and P. Zurita, In depth analysis of the combined HERA data in the dipole models with and without saturation, Phys. Rev. D 98, 036002 (2018). [7] N. Armesto, T. Lappi, H. Mäntysaari, H. Paukkunen, and M. Tevio, Signatures of gluon saturation from structurefunction measurements, Phys. Rev. D 105, 114017 (2022). [8] F. Gelis, E. Iancu, J. Jalilian-Marian, and R. Venugopalan, The color glass condensate, Annu. Rev. Nucl. Part. Sci. 60, 463 (2010). [9] H. Weigert, Evolution at small xBj: The color glass condensate, Prog. Part. Nucl. Phys. 55, 461 (2005). [10] N. N. Nikolaev and B. G. Zakharov, Color transparency and scaling properties of nuclear shadowing in deep inelastic scattering, Z. Phys. C 49, 607 (1991). [11] A. H. Mueller, Unitarity and the BFKL pomeron, Nucl. Phys. B437, 107 (1995). [12] J. L. Albacete, N. Armesto, J. G. Milhano, P. Quiroga-Arias, and C. A. Salgado, AAMQS: A non-linear QCD analysis of new HERA data at small-xincluding heavy quarks, Eur. Phys. J. C 71, 1705 (2011). PHYSICAL REVIEW LETTERS 130, 192301 (2023) 192301-5 [13] T. Lappi and H. Mäntysaari, Single inclusive particle production at high energy from HERA data to protonnucleus collisions, Phys. Rev. D 88, 114020 (2013). [14] B. Duclou´e, E. Iancu, G. Soyez, and D. N. Triantafyllopoulos, HERA data and collinearly-improved BK dynamics, Phys. Lett. B 803, 135305 (2020). [15] R. D. Ball, V. Bertone, M. Bonvini, S. Marzani, J. Rojo, and L. Rottoli, Parton distributions with small-xresummation: Evidence for BFKL dynamics in HERA data, Eur. Phys. J. C 78, 321 (2018). [16] P. Tribedy and R. Venugopalan, QCD saturation at the LHC: Comparisons of models to p þp and A þA data and predictions for p þPb collisions, Phys. Lett. B 710, 125 (2012);718, 1154(E) (2013). [17] A. M. Stasto, B.-W. Xiao, and D. Zaslavsky, Towards the Test of Saturation Physics Beyond Leading Logarithm, Phys. Rev. Lett. 112, 012302 (2014). [18] J. L. Albacete, P. Guerrero Rodríguez, and Y. Nara, Ultraforward particle production from color glass condensate and Lund fragmentation, Phys. Rev. D 94, 054004 (2016). [19] H. Mäntysaari and H. Paukkunen, Saturation and forward jets in proton-lead collisions at the LHC, Phys. Rev. D 100, 114029 (2019). [20] Y. Shi, L. Wang, S.-Y. Wei, and B.-W. Xiao, Pursuing the Precision Study for Color Glass Condensate in Forward Hadron Productions, Phys. Rev. Lett. 128, 202302 (2022). [21] B. Schenke, P. Tribedy, and R. Venugopalan, Fluctuating Glasma Initial Conditions and Flow in Heavy Ion Collisions, Phys. Rev. Lett. 108, 252301 (2012). [22] Y. V. Kovchegov, Small xF 2structure function of a nucleus including multiple Pomeron exchanges, Phys. Rev. D 60, 034008 (1999). [23] I. Balitsky, Operator expansion for high-energy scattering, Nucl. Phys. B463, 99 (1996). [24] G. Beuf, H. Hänninen, T. Lappi, and H. Mäntysaari, Color glass condensate at next-to-leading order meets HERA data, Phys. Rev. D 102, 074028 (2020). [25] H. Abramowicz et al. (H1, ZEUS Collaboration), Combination and QCD analysis of charm and beauty production cross-section measurements in deep inelastic ep scattering at HERA, Eur. Phys. J. C 78, 473 (2018). [26] G. Beuf, T. Lappi, and R. Paatelainen, Massive quarks in NLO dipole factorization for DIS: Transverse photon, Phys. Rev. D 106, 034013 (2022). [27] G. Beuf, T. Lappi, and R. Paatelainen, Massive Quarks at One Loop in the Dipole Picture of Deep Inelastic Scattering, Phys. Rev. Lett. 129, 072001 (2022). [28] G. Beuf, T. Lappi, and R. Paatelainen, Massive quarks in NLO dipole factorization for DIS: Longitudinal photon, Phys. Rev. D 104, 056032 (2021). [29] H. Hänninen, T. Lappi, and R. Paatelainen, One-loop corrections to light cone wave functions: the dipole picture DIS cross section, Ann. Phys. (Amsterdam) 393, 358 (2018). [30] G. Beuf, Dipole factorization for DIS at NLO: Loop correction to the γ T;L →q¯ qlight-front wave functions, Phys. Rev. D 94, 054016 (2016). [31] G. Beuf, Dipole factorization for DIS at NLO: Combining the q¯ qand q¯ qg contributions, Phys. Rev. D 96, 074033 (2017). [32] Y. V. Kovchegov, J. Kuokkanen, K. Rummukainen, and H. Weigert, Subleading-Nccorrections in non-linear small-x evolution, Nucl. Phys. A823, 47 (2009). [33] B. Duclou´e, E. Iancu, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Non-linear evolution in QCD at highenergy beyond leading order, J. High Energy Phys. 04 (2019) 081. [34] B. Duclou´e, H. Hänninen, T. Lappi, and Y. Zhu, Deep inelastic scattering in the dipole picture at next-to-leading order, Phys. Rev. D 96, 094017 (2017). [35] I. Balitsky and G. A. Chirilli, Next-to-leading order evolution of color dipoles, Phys. Rev. D 77, 014019 (2008). [36] T. Lappi and H. Mäntysaari, Direct numerical solution of the coordinate space Balitsky-Kovchegov equation at next to leading order, Phys. Rev. D 91, 074016 (2015). [37] T. Lappi and H. Mäntysaari, Next-to-leading order Balitsky- Kovchegov equation with resummation, Phys. Rev. D 93, 094004 (2016). [38] E. Iancu, J. D. Madrigal, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Resumming double logarithms in the QCD evolution of color dipoles, Phys. Lett. B 744, 293 (2015). [39] E. Iancu, J. D. Madrigal, A. H. Mueller, G. Soyez, and D. N. Triantafyllopoulos, Collinearly-improved BK evolution meets the HERA data, Phys. Lett. B 750, 643 (2015). [40] G. Beuf, Improving the kinematics for low-xQCD evolution equations in coordinate space, Phys. Rev. D 89, 074039 (2014). [41] G. Beuf, H. Hänninen, T. Lappi, and H. Mäntysaari, Color glass condensate at next-to-leading order meets HERA data, 2020, 10.5281/zenodo.4229269. [42] F. D. Aaron et al. (H1, ZEUS Collaboration), Combined measurement and QCD analysis of the inclusive ep scattering cross sections at HERA, J. High Energy Phys. 01 (2010) 109. [43] H. Mäntysaari and J. Penttala, Complete calculation of exclusive heavy vector meson production at next-to-leading order in the dipole picture, J. High Energy Phys. 08 (2022) 247. [44] H. Mäntysaari and J. Penttala, Exclusive heavy vector meson production at next-to-leading order in the dipole picture, Phys. Lett. B 823, 136723 (2021). [45] I. Balitsky, Quark contribution to the small-xevolution of color dipole, Phys. Rev. D 75, 014001 (2007). [46] H. Abramowicz et al. (H1, ZEUS Collaboration), Combination of measurements of inclusive deep inelastic ep scattering cross sections and QCD analysis of HERA data, Eur. Phys. J. C 75, 580 (2015). [47] J. L. Albacete and Y. V. Kovchegov, Solving high energy evolution equation including running coupling corrections, Phys. Rev. D 75, 125021 (2007). [48] V. Andreev et al. (H1 Collaboration), Measurement of inclusive ep cross sections at high Q2at ffiffiffi s p¼225 and 252 GeV and of the longitudinal proton structure function FLat HERA, Eur. Phys. J. C 74, 2814 (2014). [49] R. Boussarie, A. V. Grabovsky, D. Y. Ivanov, L. Szymanowski, and S. Wallon, Next-to-Leading Order Computation of Exclusive Diffractive Light Vector Meson PHYSICAL REVIEW LETTERS 130, 192301 (2023) 192301-6 Production in a Saturation Framework, Phys. Rev. Lett. 119, 072002 (2017). [50] P. Caucal, F. Salazar, and R. Venugopalan, Dijet impact factor in DIS at next-to-leading order in the color glass condensate, J. High Energy Phys. 11 (2021) 222. [51] H. Mäntysaari and J. Penttala, Exclusive production of light vector mesons at next-to-leading order in the dipole picture, Phys. Rev. D 105, 114038 (2022). [52] G. Beuf, H. Hänninen, T. Lappi, Y. Mulian, and H. Mäntysaari, Diffractive deep inelastic scattering at NLO in the dipole picture: The q¯ qg contribution, Phys. Rev. D 106, 094014 (2022). [53] F. Bergabo and J. Jalilian-Marian, One-loop corrections to dihadron production in DIS at small x,Phys. Rev. D 106, 054035 (2022). [54] F. Bergabo and J. Jalilian-Marian, Single inclusive hadron production in DIS at small x: Next to leading order corrections, J. High Energy Phys. 01 (2023) 095. [55] P. Caucal, F. Salazar, B. Schenke, and R. Venugopalan, Back-to-back inclusive dijets in DIS at small x: Sudakov suppression and gluon saturation at NLO, J. High Energy Phys. 11 (2022) 169. [56] E. Iancu and Y. Mulian, Forward dijets in proton-nucleus collisions at next-to-leading order: The real corrections, J. High Energy Phys. 03 (2021) 005. [57] K. Roy and R. Venugopalan, NLO impact factor for inclusive photon þdijet production in eþADIS at small x,Phys. Rev. D 101, 034028 (2020). [58] P. Taels, T. Altinoluk, G. Beuf, and C. Marquet, Dijet photoproduction at low xat next-to-leading order and its back-to-back limit, J. High Energy Phys. 10 (2022) 184. PHYSICAL REVIEW LETTERS 130, 192301 (2023) 192301-7