High-energy dipole scattering amplitude from evolution of low-energy proton light-cone wave functions

The forward scattering amplitude of a small dipole at high energies is given in the mean field approximation by the Balitsky-Kovchegov (BK) evolution equation. It requires an initial condition $N(r; x_0)$ describing the scattering of a dipole with size $r$ off the target that is probed at momentum fraction $x_0$. Rather than using ad hoc parameterizations tuned to high-energy data at $x\ll x_0$, here we attempt to construct an initial scattering amplitude that is consistent with low-energy, large-$x$ properties of the proton. We start from a non-perturbative three quark light-cone model wave function from the literature. We add ${\cal O}(g)$ corrections due to the emission of a gluon, and ${\cal O}(g^2)$ virtual corrections due to the exchange of a gluon, computed in light-cone perturbation theory with exact kinematics. We provide numerical data as well as analytic parameterizations of the resulting $N(r; x_0)$ for $x_0=0.01 - 0.05$. Solving the BK equation in the leading logarithmic (LL) approximation towards lower $x$, we obtain a fair description of the charm cross section in deeply inelastic scattering measured at HERA by fitting one parameter, the coupling constant $\alpha_s\simeq 0.2$. However, without the option to tune the initial amplitude at $x_0$, the fit of the high precision data results in $\chi^2/N_\text{dof} = 2.3$ at $N_\text{dof} =38$, providing clear statistical evidence for the need of systematic improvement e.g. of the photon wave function, evolution equation, and initial condition.


I. INTRODUCTION
In Deep Inelastic Scattering (DIS) a pointlike virtual photon probes the rich QCD dynamics taking place inside the proton or a nucleus. At high energies, where the small Bjorken-x part of the target wave function is probed, one observes very large gluon densities [1]. When the gluon densities become of the same order as inverse coupling, non-linear QCD dynamics start to dominate and multiple scattering effects are important [2]. In the high-energy limit, the scattering process is most conveniently described in the dipole picture in a frame where the virtual photon has a large momentum [3], and its partonic Fock states, such as |qq at leading order (LO), have a long lifetime as they scatter from the color field of the target.
Describing the QCD dynamics in this high-density domain is natural in the Color Glass Condensate [4] framework. Here the center-of-mass energy or Bjorken-x dependence of various observables (and as such the target structure) is described in the large-N c limit by the perturbative Balitsky-Kovchegov (BK) renormalization group equation [5,6]. It describes how the dipole-target scattering amplitude, which contains information about the target structure, changes with increasing energy. The dipole amplitude (a correlator of two Wilson lines) is actually a convenient degree of freedom at high energies: all cross sections computed at high energy in the CGC framework are expressed in terms of the dipole amplitude or higher-point correlators which can be written, in a Gaussian approximation, in terms of the two-point dipole amplitude [7].
The initial condition for the dipole-proton scattering amplitude depends on non-perturbative properties of the proton. A typical approach in the field has been to assume an intuitive functional form at an initial x 0 1 and fit various unknown parameters to the HERA total cross section data; see, e.g., Refs. [8][9][10] where a very good description of small-x HERA data is obtained at leading order, resumming powers of α s ln 1/x via BK evolution with running coupling corrections [11]. Recent developments to full NLO accuracy have also allowed for a simultaneous description of total and heavy quark production data [12,13]. The drawback of this approach is that one is sensitive to the assumed functional form of the initial dipole amplitude and that the model parameters need to be re-fitted if the evolution is initialized at different x 0 . Furthermore, there is no relation to the low energy ( or "large-x") proton structure.
In this work, we take a complementary approach aim-arXiv:2303.16339v2 [hep-ph] 26 Jun 2023 ing to compute the initial dipole-proton scattering amplitude at moderate x 0 . As we will discuss in more detail next, the necessary non-perturbative input consists in a proton valence quark wave function that is constrained by low-energy data. The x 0 -dependent initial condition is then obtained by computing the dipole-target scattering amplitude including one perturbative gluon emission in the target, with the gluon longitudinal momentum fraction regulated by x 0 [14,15]. The advantages of this approach are that we do not assume an ad-hoc functional form of the scattering amplitude and that the initial condition can be computed and the BK evolution initialized at any (moderate) x 0 without a need to perform new fits. Also, this approach largely eliminates the freedom of tuning initial conditions in order to optimally match the evolution equation to the small-x data. This may reveal quantitative evidence for the need for improvements beyond leading-log, or even running coupling BK evolution. Finally, we would like to point out that light-cone Hamiltonian calculations of wave functions have been employed previously to set initial conditions for QCD scale evolution to high virtuality Q 2 , in order to describe DIS in the "dipole approach" using correlators of eikonal Wilson lines as degrees of freedom [16][17][18][19]. Our approach is similar in spirit although here the goal is to determine initial conditions for evolution to small x.

II. DIPOLE-PROTON SCATTERING AT MODERATE x
We first provide an overview of our approach to the light-cone structure of the proton. We employ a truncated Fock space description which starts with a three quark state. The corresponding Fock space amplitude (wave function) Ψ qqq corresponds to a non-perturbative solution of the QCD light-front Hamiltonian. To date, exact solutions for the light-cone wave functions are not available. In the future, lattice gauge theory may provide numerical solutions for moderate parton momentum fractions x i and transverse momenta k i via a large momentum expansion of equal-time Euclidean correlation functions in instant quantization [20,21]; see ref. [22] for a recent lattice computation of the wave function of the leading qq state of the pion. Also, the MAP collaboration [23] has recently extracted the wave functions of the first four Fock states of the pion from fits to its parton distribution functions and electromagnetic form factor.
Here, we rely on solutions of effective light-cone Hamiltonians for guidance on the low energy and low virtuality Q 2 structure of the proton. Specifically, we shall employ the HO wave function of Refs. [24,25]. In these references, the authors fixed the parameters of the three quark wave function to the proton "radius", or Dirac form factor at Q 2 → 0, to the anomalous magnetic moments of the proton and neutron, and to the axial vector coupling g A . The wave function also matches rea-sonably well the empirical knowledge of the longitudinal and transverse momentum distribution of single quarks in the valence quark regime. Finally, the wave function of Refs. [24,25] also provide predictions for quark momentum correlations.
At next-to-leading order (NLO) in the Fock expansion we add the three quarks and one gluon state with amplitude Ψ qqqg , as well as the virtual corrections to Ψ qqq due to the exchange of a gluon by two quarks in the proton. These corrections are obtained via light-cone perturbation theory calculations [14,15]. The presence or exchange of the gluon extends the range of parton lightcone momentum fractions to lower x, and pushes their transverse momenta into the perturbative regime. It also affects their momentum correlations.
The central element of our analysis is the (imaginary part of the) eikonal scattering amplitude N of a small dipole of transverse size r. The real part of N corresponds to two-gluon exchange, Here K is the momentum transfer which is Fourier conjugate to the impact parameter b. As explained below, we will eventually average N (r, b) over a suitable range of impact parameters. We emphasize that the expression above accounts only for a single, perturbative two-gluon exchange (see its derivation in Ref. [26]), it does not resum the Glauber-Mueller multiple scattering series. This restricts its applicability to the regime of weak scattering, N (r, b) 1. Furthermore, N (r, b) actually acquires an imaginary part due to the perturbative exchange of three gluons; its magnitude has been shown to be much smaller than its real part [27,28], and in practice it is of interest only for processes involving C-conjugation odd exchanges [29]. For the present purposes, it can be neglected.
The coupling of the two static gluons to the proton is described in terms of the color charge density correlator The color charge density operator corresponds to the light-cone plus component of the color current on the x + = 0 light front, integrated over x − , ρ a (q) ≡ J +a (q), when the proton carries positive P z . Eqs. (1,2) correspond to the leading twist contribution to the matrix element of the dipole operator in the proton. Dozens of diagrams contribute to this correlator at NLO, their explicit expressions are listed in Ref. [14]. We point out that G 2 (q 1 , q 2 ) satisfies a Ward identity due to the color neutrality of the proton; it vanishes when either q 2 1 or q 2 2 → 0 so that N (r, b) in Eq. (1) is free of IR divergences. However, G 2 does exhibit a collinear singularity which is regularized by assigning a mass to the quarks in the light-cone energy denominators for the q → qg and qg → q vertices; see Ref. [14] for details. All the results presented here were obtained with m coll = 0.2 GeV. This is consistent with the quark mass and transverse momentum scales which appear in the non-perturbative three quark wave function of Refs. [24,25]. The color charge correlator also exhibits a soft singularity when the lightcone momentum fraction x g of the gluon goes to zero. This is regularized with a cutoff x on x g , and the resummation of yet softer gluons will be performed through the BK equation. Note that at x = 0.1 the NLO contribution to G 2 truly is a reasonably small O(g 2 ) perturbative correction [15]. However, by x = 0.01 its magnitude grows to essentially O(1), a leading-log correction. Hence, at such x resummation is required and it is justified to use the computed dipole as an initial condition for the leading order BK evolution.
We recall, also, that at the given order ultraviolet divergences cancel [14], so that G 2 is independent of the renormalization scale, and that the coupling does not run. Lastly, let us mention that the angular dependence of the correlator G 2 , as well as the dependence of its Fourier transform on impact parameter, has been analyzed numerically in detail in Ref. [15].

III. SMALL-x EVOLUTION OF THE PROTON LIGHT-CONE WAVE FUNCTION
In order to obtain an initial condition for bindependent BK evolution 1 we average the dipole-target scattering amplitude obtained from Eq. (1) over the impact parameter b, Throughout this work, we denote the magnitudes of the transverse vectors as b = |b| and r = |r|. The resulting amplitude is dominated by perturbative contributions when the dipole size r is small. In this region there is a small cos(2φ) dependence on the angle φ between r and b [15] which vanishes when we integrate over b. Here S T is the proton transverse area. Inclusive cross sections considered in this work are not sensitive to the actual shape of the target but only to the total transverse size. The proton geometry is most directly probed in exclusive vector meson production process where the total momentum transfer K which is Fourier conjugate to the impact parameter is measurable. Parametrizing the J/ψ production cross section in HERA kinematics as e −B D K 2 one obtains B D = 4 GeV −2 [32]. Assuming a Gaussian impact parameter profile for the proton, this corresponds to a two-dimensional root-meansquare radius b Gaussian = √ 2B ≈ 0.56 fm and a proton area S T = 2πB. On the other hand, if we assume a step function (hard sphere) profile for the proton, the same diffractive slope is obtained when the proton radius is b Hard sphere = 2 √ B ≈ 0.79 fm, which corresponds to S T = 4πB.
Although exclusive vector meson data favors the Gaussian density profile over the hard sphere one (see e.g. [33]), the current data does not constrain the proton shape precisely. We also note that if the b-dependent dipole amplitude from Eq. (1) is directly used to compute exclusive J/ψ production cross section, the resulting spectra differs from the Gaussian profile case only in the region where there are no experimental constraints [34]. In this work the results shown below by default correspond to the Gaussian density profile (with b max = b Gaussian ) unless otherwise stated, but we also study the dependence on the b max cut by using a step function profile with b max = b Hard sphere = √ 2b Gaussian . The proton transverse area S T has also been extracted by fitting a parameterized initial condition for the BK evolution equation to the HERA structure function data. Leading order analyses [9,10] typically obtain S T ∼ 16 mb. In recent fits at NLO accuracy [12,13] proton areas S T ∼ 10 . . . 20 mb were obtained depending on the details of the analysis setup. We test this uncertainty in the proton small-x transverse profile by showing some results for both the Gaussian and hard sphere profiles with transverse areas 9.8 mb and 19.6 mb, respectively.
Before performing the impact parameter average we first study the impact parameter profile from the NLO light-cone wave function seen by a perturbative probe: The normalization condition bmax d 2 b T (b) = 1 is used to fix the constant C. We will refer to T (b) as the transverse "density" profile to match standard terminology from the literature. As the dipole amplitude is a rapidly increasing function of the dipole size r, this integral is dominated by dipoles of size r ∼ r max , as long as r max is in the perturbative domain. The extracted density profiles up to b max = 0.8 fm for different r max are shown in Fig. 1. For reference, a Gaussian profile, as used e.g. in the popular IPsat parametrization [35] for the dipole amplitude with the slope B = 4 GeV −2 , is also shown. We observe a similar transverse profile except for very central b 0.2 fm where the computed profile is more steeply falling. This region can only be probed at high momentum transfer |t| > ∼ 1 GeV 2 [34], which is not covered in the currently available coherent vector meson production data. The high-b tails of T (b) resulting from the LCPT one gluon emission corrections are exponential rather than Gaussian. However, in all we conclude that for the present  purposes the Gaussian profile used to match S T to b max is a reasonable approximation.
The b-averaged dipole amplitudes (using a Gaussian profile) are shown in Fig. 2 (linear scale) and Fig. 3 (logarithmic scale) at x 0 = 0.05, x 0 = 0.025 and x 0 = 0.01. Here we also show the dependence on the diffractive slope B: the bands correspond to varying B by ±10% which changes both S T and b max . The results depend weakly on this cut especially in the perturbative small-r domain. The dipole amplitude increases with r, approximately proportional to r 2 , as expected. For r > ∼ 0.4 fm the color neutrality of the proton, and the fact that the dipole scatters from a target of finite transverse extent, begin to slow the growth of N (r); a model that does not ac- count for the finite size of the proton in impact parameter space would attribute this to power corrections. Finally, when the size of the dipole becomes comparable to that of the target the amplitude is found to decrease again (not shown) as the end points of the dipole essentially "miss" the target. However, we emphasize that this behavior occurs at large r ∼ few fm where in any case the perturbative calculation of the scattering amplitude is not valid. Figures 2 and 3 confirm that down to x = 0.01 scattering of small dipoles with r significantly less than 1 fm remains quite weak, at least for α s = 0.2 which we determine below from a fit to the charm cross section in DIS. Therefore, it appears reasonable to start small-x evolution with this initial condition at x in the range 0.01 . . . 0.05.
To obtain analytic parameterizations of the dipole amplitude we fit our numerical data for the b-averaged scattering amplitude to the following expression which is inspired by the McLerran-Venugopalan (MV) model [36]: where Λ = 0.241 GeV is a fixed infrared scale. Such a parameterization has been used previously e.g. in Refs. [9,10,12] to fit the initial condition for BK evolution to the HERA data. While our fit is restricted to r < 0.5 fm, the parameterization forces N (r) → 1 in the large-r region. Of course, the behavior at large r can not be trusted, and other extrapolations would be possible. It is important, however, that the large-r extrapolation is such that the Fourier transform of 1 − N (r) at high k (and, consequently, the forward particle production cross section, for example) will be insensitive to the assumed form. We will also demonstrate below that perturbative observables, in our case the charm production cross section, are only sensitive to the perturbative regime of  small dipoles where our calculation should apply, and not to the extrapolation to large r. The free parameters in Eq. (5), Q 2 s,0 , γ and e c , are fit to the calculated dipole amplitude in the region 0.01 fm < r < 0.5 fm (we actually fit the logarithm of the dipole in order to give equal weight to small and intermediate r). The upper limit restricts to the perturbative domain, and the lower limit is imposed in order to give some weight to the region of intermediate r as well. The resulting dipole amplitudes are shown in Figs. 2 and 3 as dotted lines. The fit parameters are listed in Table I. In the fit we require that e c > e −1 in order to enforce positivity of the logarithm in Eq. (5), and all fit results give e c = e −1 within numerical accuracy, i.e. they require as small an infrared cutoff as allowed. The MV-model inspired parameterization is found to describe the dipoleproton scattering amplitude quite well, for all dipole sizes in the perturbative r 0.5 fm region. Here, of course, the linearized version of Eq. (5) is sufficient, as it should be: recall that Eq. (1) does not resum multiple scattering.
The momentum scale Q 2 s,0 remains non-perturbative down to x = 0.01; see below for an extraction of a "saturation scale" at lower x. However, it increases approximately as Q 2 s,0 ∼ 1/x 0.47 . The "anomalous dimension" of the dipole amplitude is γ = 1 within numerical accuracy. It appears reasonable to us that the initial condition for the evolution equation admits a power series expansion in r 2 , starting at its first power. On the other hand, leading order fits to HERA total cross section data [9,10] require γ ∼ 1.1 − 1.2 in order to obtain as slow a Q 2 dependence of the cross section as required by the HERA data [1,37]; recent fits at next-to-leading order accuracy performed in Ref. [12,13] also prefer γ 1 when the heavy quark production data is included. A problem with γ > 1 is that it renders the (dipole) unintegrated gluon distribution function [7,38,39] and the forward particle production cross section negative [10,40] in some range of transverse momentum k T . The dipole amplitude obtained here does not display this issue.
Next we solve the leading order BK equation with fixed coupling, using the numerical data for N (r, x 0 ) as an initial condition at x 0 = 0.01. Note that at this order in α s the coupling constant does not run in the LCPT calculation of the initial condition, and consequently we also limit ourselves to the fixed coupling case here. Evolution over 6 units of rapidity is shown in Fig. 4. For comparison, we also solve the BK equation using the modified MV-model initial condition with parameters as shown in Table I. This parameterized initial condition has a completely different behavior in the infrared region with N (r) → 1 at large r whereas the numerical data gives a decreasing N (r) when r exceeds a few fm, as already mentioned above. However, as can be seen in Fig. 4 the resulting BK-evolved dipole amplitudes are basically identical in the perturbative r 0.5 fm domain. In fact, due to the approach to the fixed point of the BK equation [41][42][43], at high rapidity the difference between the scattering amplitudes evolved with the two initial conditions diminishes. This demonstrates that the BK-evolved amplitude at small r is not affected by the uncontrolled large-r extrapolation of the initial condition.
One may define a saturation radius r s , and a corresponding saturation momentum Q s = √ 2/r s from the condition that N (r s ) = 1 − exp − 1 2 0.4. For this to be a perturbative scale requires about 6 units of rapidity evolution, as can also be seen from Fig. 4. This corresponds to x 2.5 · 10 −5 , where r s 0.3 fm, and Q s 1 GeV. These values are not very far from the first "saturation model" fit to HERA DIS data by Golec-Biernat and Wüsthoff [44] from 25 years ago. Many more recent fits mentioned above have since confirmed that reaching the strong scattering regime with a small dipole and a proton target requires deep evolution to rather small x. Also, some studies [17] of diffractive small-x scattering of a qq−g state from the proton have indicated that the regime of "color transparency" sets in when the typical transverse distance of the gluon from the qq is around 0.2 fm.
Since scattering at x = 0.01 is fairly weak we have also evolved our initial condition with the linear BFKL . After a few units of rapidity evolution, the linear equation begins to violate unitarity, N (r) ≤ 1, at large r. However, this regime of large dipoles is not under control in any case. More importantly though, at y = 2 − 4 the absence of the nonlinear correction begins to affect the solution significantly even at r substantially less than 1 fm. With BFKL we also noticed a greater difference between evolving the actual numerical data vs. the analytic modified MV-model parametrization (not shown), which differ in their large-r extrapolation. Therefore, for accurate results it appears to be rather important to evolve with the non-linear BK equation even if one restricts to r < 1 fm.
Let us finally study how the x dependence obtained from the direct, fixed order NLO LCPT calculation differs to the one obtained by solving the BK equation. We note that in the LCPT calculation x is an explicit cutoff for the longitudinal momentum of the emitted gluon, and this gluon emission is calculated in exact kinematics. On the other hand, in BK evolution multiple soft gluon emissions are resummed. This comparison is done by calculating the dipole amplitude at x = 0.01 directly from the LCPT using Eq. (1), and comparing that to the dipole amplitude obtained by solving the BK equation with the initial condition computed at x 0 = 0.05. The resulting dipole amplitudes are shown in Fig. 6. The most significant difference between the fixed order O(g 2 ) LCPT amplitude and the BK evolved dipole is that the evolution equation decreases the anomalous dimension γ towards the asymptotic value γ ∼ 0.6. On the other hand, the emission of one gluon in the direct LCPT calculation does not modify the extracted anomalous dimension, as can also be seen from Table. I. This is, of course, the expected behavior. As already mentioned above, DIS phenomenology does not appear to support γ < 1 at x = 0.01 or greater, so it seems reasonable to treat at least the emission of the first gluon with x g > 0.01 ex- plicitly in fixed order light-cone perturbation theory with exact kinematics 2 .

IV. TOTAL CROSS SECTION AT SMALL x
Next, we consider the DIS structure functions at small Bjorken-x. The overall normalization of the dipole amplitude depends on the strong coupling constant α s = g 2 /(4π), see Eq. (1). The same coupling constant also affects the Bjorken-x dependence of the dipole amplitude via the BK evolution. In this work, our strategy is to fix the value of α s by calculating the total charm production cross section, and comparing it to the HERA reduced cross section data from Ref. [48]. We set the initial condition for the BK evolution at x 0 = 0.01, and compare it to the HERA data in the region x < 0.01, Q 2 < 100 GeV 2 (note that the smallest Q 2 bin in the data is Q 2 = 2.5 GeV 2 ). In this region, there are N = 39 data points. The experimental data is reported as reduced cross section Here y = Q 2 /(sx) is the inelasticity variable, not to be confused with the evolution rapidity. The proton structure functions F 2 and F L are expressed in terms of the total cross section for the virtual photon-proton cross section σ γ * p : In the dipole picture, the total cross section for the virtual photon-proton scattering can be expressed as a convolution of the photon wave function and the dipole amplitude as [2] σ γ * A T,L = 2 (9) Here f is the quark flavor, Q 2 the photon virtuality and Ψ γ * →qq is the leading order wave function for the qq Fock state of the virtual photon. In this equation we replace N (r, b,x) by the impact parameter averaged dipole amplitude N (r,x), as described above, and d 2 b → S T . We fix the mass of the c quark to 1.4 GeV. The dipole amplitude in Eq. (9) is evaluated atx = x(1 + 4m 2 q /Q 2 ), where m q is the quark mass which enforces a smooth approach to the photoproduction limit [9,44].
In order to confirm that the charm production cross section is not sensitive to non-perturbatively large dipoles we show in fig. 7 the fraction of the total cross section at x = 0.0056, Q 2 = 10 GeV 2 as a function of the upper limit r max for the r integral in Eq. (9). It is evident that the charm cross section is saturated by small dipoles whereas the inclusive cross section (calculated using m q = 0.14 GeV for the light quarks) at Q 2 = 10 GeV 2 is sensitive to larger dipoles beyond sizes where we may trust our calculation. When using a modified MV-model parameterization as an initial condition for the evolution with different extrapolation in the infrared region, one needs to integrate up to even larger r in order to recover the full result for F 2 . The charm production cross section is not affected by the different large-r extrapolation (not shown). Qualitatively similar results have been obtained with the commonly used IPsat parameterization for the dipole-proton amplitude where, typically, even larger dipole sizes contribute as compared to the setup with factorized impact parameter dependence applied here [35,49]. For these reasons we shall focus on charm production. In the future, our approach could be applied to other hard, perturbative processes such as, for example, single-inclusive particle production at high enough transverse momentum. Considering only the strong coupling constant α s as free parameter we obtain a reasonably good description of the charm production data. The value of χ 2 /(N − 1) as a function of α s is shown in Fig. 8 using two different density profiles (Gaussian and hard sphere) for the proton. These two setups have different upper limits for the impact parameter b and correspondingly different transverse areas for the proton. The extracted optimal values for the strong coupling constant are α s = 0.200 for the Gaussian proton and α s = 0.181 for the hard sphere profile. These values are used throughout this paper. We note that fits of the total (rather than charm) cross section with tuned initial conditions [9,10] and running coupling corrections to the BK equation have achieved lower χ 2 /N ≈ 1, without being able to simultaneously describe the charm data [9]. However, with our calculated initial condition there is room for the expected improvements of the photon wave function, evolution equation, and, of course, of the initial condition.
In this analysis, we have fixed the collinear regulator to m coll = 0.2 GeV in the LCPT calculation of the initial condition. As the charm cross section is dominated by small dipoles, our results are not highly sensitive to the actual value of this regulator: changing m coll by a factor of 2 changes χ 2 /(N − 1) to HERA data by only 2 . . . 5% when using the optimal α s . We also keep the charm mass fixed to 1.4 GeV. The optimal value for α s and the achieved χ 2 /(N −1) naturally will depend on this choice. We choose to work with fixed quark mass and collinear regulator and do not attempt to fit these simultaneously with α s , as the purpose of this work is to demonstrate the feasibility of computing the initial condition for the BK equation, and we emphasize that numerically potentially important higher order effects such as running coupling are still missing from the setup.
A comparison to the HERA charm production data in different Bjorken-x bins is shown in Fig. 9 as a function of the photon virtuality. We have checked that these results remain the same if we use the analytic parameterization (5), with parameters from table I, as initial condition; this confirms the insensitivity of the charm cross section to the large-r extrapolation of the scattering amplitude.
At x = 0.008 there is only very little (≤ 0.2 units of rapidity when x 0 = 0.01) evolution, so the dipole amplitude is almost completely determined by our initial condition. On the other hand, we also show results at lower x which is dominated by the BK evolution. In addition to our standard setup where the initial condition for the BK evolution is set at x 0 = 0.01, we also show results using an initial condition computed at larger x 0 = 0.05. Note the weak dependence of this observable, at least, on where the "hand-off" from the x-dependent initial condition to the evolution equation occurs. In contrast, ad-hoc initial condition parametrizations have to be re-tuned when x 0 is changed. Fig. 9 shows a fair agreement of the reduced cross section obtained from our light-cone wave function with the HERA charm data. Close to the initial condition we obtain a slightly slower Q 2 dependence than seen in the data. As a result of the evolution this changes into faster virtuality dependence at very small x. This is because the BK evolution at fixed coupling develops a small anomalous dimension γ ≈ 0.6 for the dipole amplitude and a smaller anomalous dimension results in faster Q 2 dependence.
Lastly, in Fig. 10 we study how sensitive the charm production cross section is on the chosen proton density profile, and as such on the maximum impact parameter b max used in Eq. (3). The cross section is calculated at x = 0.008 which is close to the initial condition for the BK evolution, again set at x 0 = 0.01. In both cases, we use the optimal value for the strong coupling constant extracted above. The cross section increases only slightly when the hard sphere profile with larger b max is used, but the dependence on the virtuality is not affected. This weak dependence on the selected proton profile confirms that our results are not sensitive to non-perturbatively large impact parameters.

V. SUMMARY AND DISCUSSION
The present work represents a first attempt at relating the short-distance structure of the proton at high energies to its low-energy properties, covering several orders of magnitude in energy. We start from an effective three quark light-cone wave function which models the nonperturbative longitudinal and transverse momentum distributions of quarks at x > ∼ 0.1 as well as some of their correlations. The next step involves the computation, using exact kinematics, in light-cone perturbation theory of the O(g) correction to the light-cone wave function due to the emission of a gluon, and the O(g 2 ) virtual corrections due to the exchange of a gluon by two quarks. This provides a leading twist contribution to the dipole scattering amplitude, N (r) ∼ r 2γ with γ = 1 at small r. Optimistically, the LCPT correction extends the validity of the resulting light-cone wave function into the regime of perturbative transverse momenta, and parton momentum fractions x > ∼ 0.01. The corresponding dipole scattering amplitude is then evolved to yet higher energies (lower x) by solving the BK equation, which resums emissions of additional soft gluons, and generates power corrections and an anomalous dimension.
The convolution of the LO photon light-cone wave function with the impact parameter averaged BK dipole scattering amplitude at leading logarithmic accuracy provides a fair description of the reduced DIS charm cross section measured at HERA, for α s 0.2. This value of the strong coupling was obtained from a fit to the charm cross section at Q 2 < 100 GeV 2 and x < 0.01. None of the parameters of the low-energy three-quark model wave function were re-tuned to the high-energy data. Despite the fair description of the highly accurate data the resulting χ 2 /N dof ≈ 2.27 with N dof = 38 implies a very low statistical significance, i.e. a very low probability that the data represents statistical fluctuations about the model predictions: the integral over the χ 2 distribution from χ 2 = 2.27 × 38 to infinity, commonly denoted as the "pvalue", is 1.3 × 10 −5 . However, the very low statistical significance of the fit should not be confused with a need for large corrections, Fig. 9 shows that this is clearly not the case. This is entirely expected since there are multiple known sources of corrections such as, for example, of the photon wave function [50][51][52], of the evolution equation [53][54][55][56][57], and, of course, of the initial condition for the evolution equation (our proton light-cone wave function) which, e.g., may be improved with running coupling corrections. The data requires fairly moderate but systematic improvements of the model predictions across the relevant ranges of x and Q 2 .
We have also provided analytic parameterizations of  the impact parameter averaged dipole scattering amplitude for x = 0.01 . . . 0.05 which accurately fit the numerical data in the regime of perturbative dipoles, r 0.5 fm. These parameterizations can be used in practice to estimate the corrections predicted by more accurate evolution equations. In the supplementary material we also provide the tabulated numerical data for N (r, x). Their large-r extrapolation differs from that of the analytic parameterizations which allows for tests of the (in-)sensitivity to the uncontrolled non-perturbative regime of large dipoles. The quest for more accurate theoretical predictions at high energy for the upcoming EIC at BNL [58] and the proposed LHeC/FCC-he at CERN [59] requires initial conditions for the evolution equations which do not absorb theoretical improvements into a re-tune of their parameters.