Decoupling of the leading contribution in the discrete BFKL Analysis of High-Precision HERA Data

We analyse, in NLO, the physical properties of the discrete eigenvalue solution for the BFKL equation. We show that a set of eigenfunctions with positive eigenvalues, \omega, together with a small contribution from a continuum of eigenfunctions with negative \omega, provide an excellent description of high-precision HERA F_2 data in the region, x<0.001, Q^2>6 GeV^2. The phases of the eigenfunctions can be obtained from a simple parametrisation of the pomeron spectrum, which has a natural motivation within BFKL. The data analysis shows that the first eigenfunction decouples completely or almost completely from the proton. This suggests that there exist an additional ground state, which is naturally saturated and may have the properties of the soft pomeron.


Introduction
The aim of this paper is to apply, for the first time, the complex BFKL Green Function approach developed in our two previous papers [1,2] to the analysis of HERA data,. The new approach, although seemingly equivalent to the discrete BFKL solution developed in our previous papers [3,4], exhibits some differences. The most important of which is that the normalisation of the eigenfunctions is now determined analytically instead of being determined only numerically, as was the case in ref. [3,4]. As shown in [2], this seemingly minor technical difference has an important consequence: the convergence of the eigenfunctions is now much more rapid than previously. Instead of using O(100) eigenfunctions, as in ref. [3,4], we need to use only O(10) eigenfunctions to properly represent the Green Function. This constrains substantially the new BFKL solution, exhibits more clearly its physical properties, and leads to new results.
To obtain a good description of the HERA F 2 data it is necessary to define a nonperturbative boundary condition defined in terms of phases of eigenfunctions at low gluon transverse momenta k, close to k ∼ Λ QCD . Since in [3,4] we were using a large number of eigenfunctions, O(100), it was easy to find a simple, ad hoc, parametrisation for these phases. However, this parametrisation had no physical interpretation.
The first task of this paper is to find a simple parametrisation for the phases of much fewer eigenfunctions, O (10). In the search for such a condition we are guided by the principle of simplicity and some analogy to the Balmer series. In the QCD version of Regge theory developed in our papers, the BFKL equation is considered to be analogous to the Schrödinger equation for the wavefunction of the pomeron. The BFKL kernel corresponds to the Hamiltonian and the eigenvalues ω to the energy eigenvalues. In this paper, we find that we can specify the boundary condition in terms of a relation between the eigenvalues ω n of the BFKL operator and the principal quantum number n. This relation then determines the boundary condition in terms of the phases η n of the eigenfunctions, close to the non-perturbative region, k ∼ Λ QCD . In addition, the relation between ω and n is very simple and, for large n, has a good physical motivation within the context of the BFKL formalism.
We show in this paper that this new approach leads to unexpected results and gives a new insight into the role of gluon density. We recall that the BFKL Green Function is directly related to the gluon density (see below). The properties of this gluon density are very interesting for the LHC and cosmic ray physics. They are also interesting in themselves, because in contrast to the DGLAP evolution [5], the BFKL equation describes a system of quasi-bound self-interacting gluons. Such a system is sensitive to confinement effects and also has some sensitivity to Super-Symmetry effects (in the gluon sector), as was first observed in ref. [3,4] and is also valid in the present approach.
The paper is organised as follows: In Section 2 we recall the main properties of the BFKL Green Function and of their eigenfunctions, determined in our last papers [1,2]. We also indicate here the differences between the approach of ref. [3,4] and our present approach. In Section 3 we introduce the NLO corrections to BFKL and evaluate the properties of eigenvalues and eigenfunctions at NLO. In Section 4 we apply this formalism to HERA data and describe the search for a proper boundary condition and the new results. Finally, in Section 5 we summarise the results and conclude.

BFKL Green Function
The Green Function approach considered here is highly appropriate since it does not require any cutoff on the BFKL dynamics and provides a direct relation to the measurements at low-x. Thus, the deep inelastic structure function F 2 (x, Q 2 ) can be directly calculated as a convolution of the Green function with impact factors that encode the coupling of the Green function to the external particles that participate in that process.
where, Y = ln(1/x), t = ln(k 2 /Λ 2 QCD ), t = ln(k 2 /Λ 2 QCD ); k, k being the transverse momenta of the gluons entering the BFKL amplitude. Φ DIS (Q 2 , t) describes the (perturbativly calculable) coupling of the gluon with transverse momentum k to a photon of virtuality Q 2 and Φ P (t ) describes the coupling of a gluon of transverse momentum k to the target proton, see In [1] we determined the BFKL Green Function G ω (t, t ) (in Mellin space) from the equation whereΩ denotes the BFKL operator, which was given in terms of the LO characteristic function, χ(α s (t), ν), bŷ with α s ≡ C A α s /π. By placing α s (t) on either side of the differential operator we assured the hermiticity of the whole operator.
We have shown in [1,2] that the Green Function determined in this way has poles on the positive real axis of the ω plane and a cut along the negative ω axis. Therefore it can be constructed from the complete set of eigenfunctions of the BFKL operator in the usual way The spectrum of the eigenvalues ω n was found to be discrete for positive values of ω and continuous for negative value of ω. The complete set of eigenfunctions with positive and negative eigenvalues ω was found to satisfy the closure relation and the orthonormality condition. In addition, the Green Function converges rapidly so it was sufficient to use only O(10) discrete eigenfunctions (see the discussion below eq.(2.17)) to describe properly the gluon density, as compared to our previous work [3,4], where we needed more than 100 eigenfunctions.

Eigenvalues and eigenfunctions
In LO BFKL [14], with fixed QCD coupling constant α S , the eigenfunctions have a simple oscillatory behaviour in terms of the gluon transverse variable t, The frequency ν of these oscillations is connected to the eigenvalue ω by the characteristic equation With fixed α S the frequency ν is a one-to-one function of ω. However, when α S is running ν becomes a function of t, ν ω (t), in order to compensate the t variation of α S . For sufficiently large values of t there is no real solution for ν ω (t) of eq.(2.6). The transition from the real to imaginary values of ν ω (t) singles out a special value of t = t c for which For values of t below the critical point t c the behaviour of the eigenfunction remains oscillatory, but above it becomes exponentially attenuated. This fixes the phase of the eigenfunction at t = t c and together with some fixed non-perturbative phase η np leads to quantisation, i.e to a discrete set of eigenfunctions.
To analyse the behaviour of the BFKL equation in the neighbourhood of the turning point, t c , it is convenient to define first two related variables, s ω (t) and z(t). The variable s ω (t) gives the phase shift from the turning point t c to the point t and corresponds to the argument of the wave function of eq.(2.5). It is defined as (2.9) and the (ω dependent) variable z(t) is defined as . (2.10) Using these variables we have shown in [1] that the BFKL operator,Ω, can be related to the "generalized Airy operator" as . (2.11) In this derivation the diffusion approximation was used in the vicinity of the turning point and the semi-classical approximation far away from it. Using these approximations we have shown [1,2] that the most general solution to equation 2.11 is given by the Green Function with Bi(z(t)) = Bi(z(t)) + cot (φ(ω)) Ai(z(t)). (2.13) Here Ai(z) and Bi(z) denote the two independent Airy functions. The function φ(ω) is defined as with η np (ω) being a non-perturbative phase, fixed at some small t 0 . From 2.12 and 2.13 it follows, as discussed in ref. [1], that the BFKL Green function has poles when φ(ω) = nπ, n = 0, 1, 2, 3 ... . The equations 2.15 and 2.14 define the eigenvalues ω n , which are a function of the nonperturbative boundary condition η np (n).
Furthermore, in [1] we have shown that, in case of positive ω n , the eigenfunctions of the BFKL operator are given by N ωn (t)Ai(z(t)), (2.16) whith N ωn (t) being the normalisation factor, which is given by .

(2.17)
Here χ denotes the BFKL characteristic function which in LO is simply equal to χ 0 but is more complicated in NLO.
The above expression is similar to the eigenfunctions used in ref. [3,4] with the difference that the normalisation factor, N ωn , was not t dependent and was determined by numerical integration. In the first paper, in which we developed our new approach [1], we argued that this difference is not very important because the t dependence of the normalisation factor is very slow and would not sizeably change the shape of the eigenfunctions. Whereas this is correct for the shapes in the physical region, it is not true for the normalisation. The numerical integration, which determines the normalisation factor, extends to very large t regions (given by t c , see Fig. 3), much above the physical region. Therefore, enhanced t dependence in [3,4] has a substantial effect when integrated over large t regions. As explained in [2] the eigenfunction of eq.(2.16) converge as 1/n 2 , whereas these of ref. [3,4] converge on a much slower pace, as 1/n.
To understand the physical meaning of the function φ it is useful to asymptotically expand the Airy function of eq.(2.16), around t = t 0 (but far away from t c ), This means that the function φ is the difference between the perturbative and non-perturbative phases of the wave function, which should not depend on t 0 . 2 For negative values of ω eq.(2.8) has no solution, i.e. there is no critical point and no quantization of eigenvalues. The negative ω eigenfunctions were derived in [2] and are given by The eigenfunctions defined by eq.(2.16) and (2.19) fulfil the completeness relation and are orthonormal, as shown in [2].

NLO evaluation
To obtain the eigenfunctions of the BFKL equation in NLO we just need to replace eq.(2.6) by its NLO counterpart where χ 0 (ν) and χ 1 (ν) are the LO and NLO characteristic functions respectively. The NLO value of α s was fixed by measurement at Z 0 pole. In our numerical analysis, we modify χ 1 following the method of Salam [8] in which the collinear contributions are resummed, leaving a remnant which is accessible to a perturbative analysis. For the analysis of this paper we use Scheme 3 of ref. [8] (see Appendix A).
To create the eigenfunctions we have chosen the value of t 0 equivalent to k 0 = 1 GeV, close to Λ QCD but still in the perturbative region, with α s (k 0 ) = 0.50. To be able to describe the measured structure function F 2 , which has a changing slope λ, η np should vary with n and the value of the non-perturbative phase η np for the leading eigenfunctions should be close to zero (see the discussion in Sections 4.1 and 4.3). We have therefore adopted the convention that n in eq.(2.15) should be counted from 1 and η np should be confined to the interval between +π/4 and −3π/4. 3 The values of η n and the corresponding eigenfunctions, used later in the fit, are not limited to this interval. They are obtained from the periodicity of η n , i.e. by adding (or subtracting) multiples of π on both sides of eq.(2.14). In the following we will label the eigenvalues and eigenfunctions with n ≥ 1 and denote the n dependent phase η np (n) simply by η n .
In Fig. 2 we display the eigenvalues ω n obtained from eq.(3.1), using three different nonperturbative phases, η n = 0, π/4, −π/4. The dotted line shows, as an example the η n = 0 case, that the dependence of ω n values from n (for n > 1) can be simply parametrised by as noticed already in [6]. For η n = 0 we found in NLO, that A = 0.52223, B= 1.62001. Since we apply this parametrisation below to describe data we recall its derivation given in ref. [6].
implies that where a, b are constants independent of ω. This leads to the relation 3.2. In NLO this relation is satisfied already for n ≥ 2, since ν ω (t 0 ) is less dependent on ω than in LO. The relation 3.2 indicates also that for large n, t c = χ 0 (0)/β 0 ω n should grow almost linearly with n. This is also a feature of the NLO computation, see Fig.3. The value of t c is related to the value of the critical momenta k c by t c = ln k 2 c /Λ 2 QCD with Λ QCD = 275 MeV. In Fig. 4 we show as example the first three different eigenfunctions 1,2 and 3, computed from eqs.(2.16) and (3.1), at phases η n = 0, π/4, −π/4.

Application to data
To apply the BFKL Green Function to data, we express the low-x structure function of the proton, F 2 (x, Q 2 ), in terms of the discrete BFKL eigenfunctions by where xg x ζ , k denotes the unintegrated gluon density (4.2) and Φ p (k) denotes the impact factor that describes how proton couples to the BFKL amplitudes at zero momentum transfer. The impact factor, Φ γ (ζ, Q, k), which describes the coupling of the virtual photon to the eigenfunctions is given in [7]; the dependence on ζ reflects the fact that beyond the leading logarithm approximation, the longitudinal momentum fraction, x, of the gluon differs from the Bjorken-value, determined by Q 2 . Φ γ (ζ, Q, k) of Ref. [7] is determined taking into account kinematical constraints allowing for non-zero quark masses. The (k /k) ωn factor arises from a mismatch between the "rapidity", Y , of the forward gluon-gluon scattering amplitude used in the BFKL approach Y = ln s kk and the logarithm of Bjorken x, which is given by This ambiguity has no effect in LO but in NLO it can be compensated by replacing the LO characteristic function χ 0 (ν) by χ 0 (ω/2, ν), which modifies the NLO characteristic function The proton impact factor is determined by the confining forces. It is therefore barely known, besides the fact that it should be concentrated at the values of k < O(1) GeV. We use here a simple parametrisation in the form which vanishes as k 2 → 0, as a consequence of colour transparency and is everywhere positive. The value of b should be around 13 GeV −2 , i.e of the inverse square of Λ QCD = 275 MeV. This is much higher than the value of b determined from data for the proton form factor, b ≈ 4 GeV −2 . Since the range of the proton impact factor is much smaller than the oscillation period of the BFKL eigenfunctions we do not expect that the results should have substantial sensitivity to a value of b. Therefore we performed the investigation assuming two very different values of the impact factor, b = 10 and b = 20 GeV −2 , corresponding to Λ QCD ≈ 320 or 220 MeV. We also used, as a check, an extreme proton impact factor,

Properties of HERA data
The HERA F 2 data in the low x region can be simply parametrised by F 2 = c (1/x) λ , with the constants c and λ being functions of Q 2 , see e.g. [10]. As Q 2 increases from 4 GeV 2 to 100 GeV 2 λ changes from about 0.15 to 0.3. The BFKL evaluation of F 2 , which assumes that η n is independent of n, would predict that λ is a constant, independent of Q 2 with λ ≈ ω 1 , since it is the first pole which dominates F 2 , when the value of η n is fixed. Therefore, the only way that λ can depend on Q 2 is if the infrared phases, η n , depend on n. Otherwise, the predicted value of λ will be about 0.25, independent of Q 2 (see Fig.2) in clear contradiction with HERA data.
The fits utilize the highest precision HERA data [9] given in terms of reduced cross sections from which we extracted the F 2 values, using the assumption that F L is proportional to F 2 . We also limit the y range in order to avoid possible complications of a larger contribution from F L (see e.g. [10]). Since we are focusing on the comparison with the F 2 measurements, we only use the 920 GeV data set of [9]. We also limited the comparison with data to the region x < 0.001 and Q 2 > 6 GeV 2 since the BFKL equation is valid at very low x only. The Q 2 cut was chosen to be relatively high to avoid any complications due to possible saturation corrections [11]. The number of experimental points used for fits was then N p = 51. (It represents around 1/3 of the whole low x data sample, defined as x < 0.01, Q 2 > 3 GeV 2 ).
For this investigation we have taken the uncorrelated errors, obtained by adding in quadrature all the correlated errors of ref. [9]. From the data analysis of ref. [11] we know that the uncorrelated errors overestimate the error sizeably, so that the χ 2 /N df of a good fit should be around 0.7, instead of about 1 as in case of correlated errors (see also [13]).

Boundary condition
The major challenge in confronting the BFKL predictions with data is the determination of the infrared boundary condition. i.e. finding the relation between the infrared phases, η n and the eigenfunction number n which generates a precise description of the data. At the beginning we tried to parametrise η as a function of n, using polynomial or other functional dependences. This failed because we were not able to find any functional dependence which would lead to χ 2 < O(500). In the next step we tried to find a set of η n (with n = 1, 2, 3...10) values using only some assumptions of local continuity. This was essentially a 10 parameter fit, with some limitations. After a longer search, using permutational methods to avoid any pre-conceptional bias on the form of η − n relation, we found a set of 10 η n values which gave an acceptable χ 2 ≈ 40. Studying this set we noticed that it can be well parametrised by an ω − n relation, similar to eq.(3.2), with a value of C which is very small, but nevertheless non-zero. The η n values were then obtained from eqs.(2.14) and (2.15), by The parameters A, B and C, together with η neg , the phase of the negative omega contribution, were considered as free parameters of the fit, which we call in the following the ABC-Fit. In addition to these four parameters the overall normalisation was also fitted to data.
As we observed that the system was exhibiting a multitude of local optima, we used the Bayesian Analysis Toolkit (BAT) [12] to find the global optimum. BAT generates samples in parameter space via Markov chain Monte Carlo (MCMC), distributed according to the posterior probability of the parameters. The best fit value is the parameter set with the highest posterior probability, corresponding to the lowest χ 2 -value. Fig. 5 shows a marginalised distribution of the ABC-Fit, for the variables, B and C. The regions of higher probability are shown as coloured areas, with probability increasing as the colour changes from blue over green to yellow. The small circle shows the position of the best fit, given in Table 1. The complicated structure of the probability distribution is also seen as a function of A and B variables, (see Fig. 6) Fig.s 5 and 6 show that the distribution of probability has a complicated structure; there are several extended regions of higher probability, which are completely disconnected. In this situation the usual fitting methods, based on MINUIT, work poorly, since they assume a smooth increase in probability towards the real minimum.
Using the BAT together with the above parametrisation we found an excellent agreement with data, χ 2 /N df ≈ 33/46. We performed this fit for several specific values of the parameters b of Φ p and found that the χ 2 values were the same, within the computational precision of the fit, ∆χ 2 = ±1. For each value of b the values of the fit parameters, A, B, C and η neg , were somewhat different and compensated the change of b (see for example Table 1). The In spite of the fact that C is very small, it is impossible to put its value to zero without seriously deteriorating the quality of the ABC fit (to χ 2 ≈ 150). In standard QCD we should expect C to be zero so that ω n → 0 when n → ∞, as in the LO calculation discussed above. However, we noticed, that the parameter C can to be set to zero if we let η 1 , the phase of the first eigenfunction, to be a free parameter, instead of C. The fits obtained in this way are of the same quality as the ABC fits, they have however an unexpected property; the value of the η 1 parameter is always chosen such that the first eigenfunction decouples (or nearly decouples) from the proton. This means that its overlap with the proton form-factor becomes zero (or nearly zero), independent of the choice of b. Therefore, we determined the phase η 1 solely from the requirement that the first eigenfunction should be orthogonal to the proton impact factor (in this way the parameters A and B are correlated, for a given impact factor, with the value of the phase η 1 ). We call this fit the AB-Fit and give its results in Table 2, for two values of b as example. 4 In the AB-Fit the first eigenfunction is not used since it is decoupling from the proton. In addition, we note that an approximate decoupling happens also in the ABC-Fit, where the contribution of the first pole is much smaller than that of the second one, by more than a factor of 10. Finally we note that in fits of Table 1 and 2 we used 20 eigenfunctions, to see the convergence (see below).  Table 2: Results of the AB-Fit to 51 data points with x < 0.001 and Q 2 > 6 GeV 2 .
The assumption of the decoupling of the first eigenfunction, together with the AB-relation of eq.(3.2), leads to a much simpler probability structure, (see Fig. 7), with a steady increase of probability towards one minimum, i.e., without a multitude of local minima. In Fig. 8 we show the η − n relation as computed from the parameters A, B of the AB-Fit for two values of b. Note that η − n relation is visibly different in the two cases, although the parameters A, B differ by a fraction of per mill only. In Fig. 9 we show the same relation as computed from the parameters A, B, C of the ABC-Fit for the same two values of b. Note that the η − n relation is simpler in the AB-Fit than in the ABC-Fit.
In general, we observe that the AB and ABC parameterisations are characterised by a high sensitivity to the values of ω. The values of the parameters A, B for the case of constant η, given below eq.(3.2), differ only by about a percent from the values in Table 2, and yet produce a very different η − n relation. A fit to data with constant η would give χ 2 ≈ 3000!

Fit results
In Fig. 10 we show the comparison of the AB-Fit results with data (with b = 10 GeV 2 ).    The BFKL Green Function, determined in our approach, is able to describe the Q 2 dependence of the data, given by the F 2 values or by the slope λ, although neither the eigenvalues nor the AB(C)-parameters are Q 2 dependent. In Fig. 11 we show the comparison of the λ parameter obtained from the AB-Fit with data. The λ parameter was determined in the very low x < 0.001 region and in the Q 2 range between 6.5 and 35 GeV 2 . The Q 2 dependence enters indirectly because the eigenfunctions depend on the transverse momentum k, which in the convolution with the photon impact factor, leads to a Q 2 dependence .

Discussion of the phase tuning mechanism
The choice of the ω − n relation determines the set of phases η n which tune the contributions of the individual eigenfunctions to describe the data. To see how this happens we display in Fig. 12 the eigenfunctions 1, 2, 3, and as an example of subleading ones the eigenfunctions 7, 8, 9, as a function of k. The eigenfunctions are plotted with the η n phases, for n ≥ 2, given by the AB-Fit. The first eigenfunction has the phase η 1 , which suppresses its overlap with the proton impact factor. The figure shows that the leading eigenfunctions 2 and 3 have the values f n (k 0 ) ≈ 0, whereas the eigenfunctions 7, 8 and 9, have the values at k 0 which are substantially different from zero.
To see more precisely how the phases determine the overlaps, we display in Fig. 13 Figure 13: Eigenfunctions 1, 2, 3, in the k region close to k 0 . The eigenfunctions are plotted with the η n phases given by the AB fit. The first eigenfunction is plotted with the phase which gives zero overlap with the proton impact factor. The fits were performed with the b-value of the proton impact factor given by 10 GeV 2 (full lines) and 20 GeV 2 (dotted lines) (dotted lines) GeV 2 . We see that, in both cases, the eigenfunction 1 starts negative at k 0 = 1 GeV but then crosses zero at k 0 ≈ 1.05 and becomes positive. This small negative region is sufficient to suppress the overlap with the proton impact factor and effectively cancel its contribution to F 2 . The eigenfunction 2 and 3 do not cross zero, and in both cases the overlap with the proton and DIS impact factors have the same signs. They give, therefore, large contributions to F 2 . The contributions of the subleading eigenfunctions 7, 8 and 9 are also significant because η n values are substantially different than zero, η 7 = 0.23, η 8 = 0.32 and η 9 = 0.42. This leads to large overlaps with the proton and photon impact factor, but in this case they have have opposite signs. Their contributions to F 2 are therefore relatively large and have negative sign so that they can generate a Q 2 dependence in the slope λ.
With exception of the contributions of the second and of the continuous negative ω terms, the contributions of other eigenfunctions are displayed as a sum of two eigenfunctions, (3+4), (5+6), ... (19+20), to simplify the picture. The black full line shows the contribution of the second, leading eigenfunction, which is substantially larger than F 2 .
The contribution of the second eigenfunction, together with the contribution (3+4) and the contribution from the continuum with negative ω, is positive. The contributions of the eigenfunctions 5 to 20 are all negative. The negative contributions correct the positive one to reproduce precisely the measured F 2 . In this way the effective slope is also changed; the contribution of the dominating, second term, which has ω 2 = 0.144, is modified to λ = 0.176 at Q 2 = 6.5 GeV 2 and λ = 0.265 at Q 2 = 35 GeV 2 , in agreement with data. Note that the contributions from the subleading eigenfunctions are much larger at Q 2 = 35 GeV 2 than at Q 2 = 6.5 GeV 2 due to the increased overlap with the DIS impact factor. Note also that the variation of the non-perturbative phases leads to a slower convergence of the sub-leading terms than in the case of a constant η, studied in ref. [2]. This is expected because the contribution of the subleading terms has to be large enough to substantially correct the leading terms in order to reproduce the data. Nevertheless, we see from Fig. 14 that the contributions of eigenfunctions with n > 16 start to approach zero., i.e. show convergence.
Summarising we can confirm that an excellent description of data is achieved by a fine tune of the non-perturbative phases η n . This phase tune is a result of a simple ω − n relation which is well motivated in BFKL and is determined by only two or three parameters.

Decoupling of the first eigenfunction and its consequences
The decoupling or near decoupling of the first eigenfunction is an unexpected and puzzling feature of this investigation. The decoupling is not connected to a particular value of the proton impact factor or to its form. The fits of Tables 1 and 2 together with the example of Fig. 13 show that when we substantially change the value of the proton impact factor, from b = 10 GeV 2 to b = 20 GeV 2 , the values of the fit parameters are re-tuned such that the resulting phases, although sightly changed, reproduce the data very well and again lead to the decoupling of the leading eigenfunction. Note, that these re-tunes hardly change the physical properties of the solution, i.e. the position of the poles, owing to the interplay between the phases eigenfunctions and the parameters AB(C). A similar result is obtained when we choose a completely different impact factor, given by a delta function, Φ p = A δ(k − k 0 ). Although this is not a realistic impact factor, the results are similar; the fit selects the phase of the first eigenfunction such that f ω 1 (k 0 ) = 0. The other parameters are re-tuned so that the fit reproduce the data with a χ 2 value close to 33, as in the fits of Tables 1 and 2.
From the technical point of view this decoupling occurs because the position of the critical point of the first eigenfunction is relatively close to the physical region, k c (1) ≈ 50 GeV, whereas the critical point of the subsequent eigenfunctions is far away from it, k c (2) ≈ 3, 3 TeV, k c (3) ≈ 270 TeV, k c (4) ≈ 20000 TeV, etc. Therefore, the first eigenfunction varies more quickly near k 0 than the subsequent ones, so that a very small change in the phase, η 1 , leads to a large change of the first contribution.
We have also checked that the results do not depend on the number of eigenfunction used in the fit, provided this number exceeds 10. In Table 3 we show the results of fits made with the first 20, 16, 12 and 10 eigenfunctions. All the fits were made with the ABC relation and in all cases the fit has chosen a phase which decouples the first eigenfunction from the proton.  Table 3: Results of the ABC-Fit performed with different number of eigenfunctions, N ef . All fits were using the same 51 data points, with x < 0.001 and Q 2 > 6 GeV 2 . The value of the proton impact factor was b = 10 GeV 2 .
We conclude therefore that the decoupling or near decoupling of the first eigenfunction is a genuine property of this analysis, independent of the choice of the proton impact factor or the number of eigenfunctions used in the fit.
It is obvious, that this decoupling can only happen because the leading eigenfunction makes a transition from the negative to positive values in a region close to the starting point k 0 . Such a transition is an indication that the first eigenfunction, as chosen by the fit, cannot be a wave function of a ground state because the ground state has to be completely positive, see Appendix B. Therefore, the decoupling of the first eigenfunction should be interpreted as an indication that there exist an additional ground state, corresponding to n = 0.
Our computation gives us some hints about the properties of such a state. From the values of the turning points, t c (n), which grows almost linearly with n, Fig. 3, we can estimate the k c value of the ground state, n = 0, as being around 700 MeV 5 , just below our starting value of k 0 = 1 GeV. Such a state would have a high intercept, ω 0 ≈ 0.3, and would not have any oscillations above k 0 , it would just decay exponentially with increasing ln(k).
As example of such a state we show in Fig.15 the momentum distribution of a state which could be similar to the real ground state and which exists in our computation. 6 It has k c = 1.05 GeV, ω = 0.37 and η = −2.35. Indeed, the k c value of the additional ground state, of around 700 MeV, lays right in the middle of the saturation region [16][17][18][19][20][21][22][23], where multiple pomeron exchanges should dominate [24]. In our approach, these exchanges would almost entirely involve the interaction of the low k c ground state with itself, since its size is much larger than the size of higher eigenfunctions and the eigenfunctions are orthogonal to each other. This will lead to unitarisation (saturation) corrections which would substantially affect the properties of the ground state. The momentum distribution will be shifted towards the lower k values and therefore its overlap with the photon impact factor should diminish quickly with increasing Q 2 . In addition, the saturation correction will damp the effective exponent of the first eigenfunction, ω ≈ 0.3, to a value which is compatible with the non-perturbative pomeron state, λ ≈ 0.1. 7 It was already pointed out by Gribov [15], in the framework of the reggeon calculus, that the soft pomeron could be given by the renormalised, bare pomeron. The renormalisation procedure should take into account the corrections due to multiple interactions. This is somewhat similar to the picture emerging from our analysis. Of course, the soft pomeron discussed by Gribov, was essentially a non-perturbative state, determined mostly by nuclear forces. 8 In our case, the bare ground state is, however, a perturbative state and its multiple interaction are also of perturbative origin. Its properties are thought determined, to large extent, by the non-perturbative, nuclear forces which enter into our analysis through the choice of the non-perturbative phase η.

Q 2 dependence
In Table 4 we show the AB-Fit results for different Q 2 regions, Q 2 > 4, 6 and 9 GeV 2 , for b = 10 GeV −2 as an example. The fits with b = 20 GeV 2 and/or the ABC fits show very similar results. The fit with Q 2 > 4 GeV 2 of Table 4   than the one with Q 2 > 6 GeV 2 . Also the fit with Q 2 > 6 GeV 2 is significantly worse than the Q 2 > 9 GeV 2 one. Therefore, it is possible that the worsening of the fit quality with decreasing Q 2 cut is due to the presence of the hypothetical ground state discussed above.

Extrapolation to very low x
In Fig. 16 we show the extrapolation of the AB fit to very low x values, which can be possibly achieved in some future ep collider like VHEeP or LHeC. We see that at very large energies the increase of F 2 shows similar slopes at different Q 2 values, unlike at HERA. This is due to the dominance of the leading trajectory at very low x values.

Conclusions and Outlook
We have shown here that there exists an infrared boundary condition, which leads to a precise description of HERA F 2 data, for x < 0.001. We formulated it in terms of a relation between the eigenvalues, ω n , and the eigenfunction number, n. It has a simple form ω n = A/(B + n) or ω n = A/(B + n) + C, called here AB or ABC relations respectively. Both relations are well motivated in BFKL, for larger n. The ω − n relation determines, within the BFKL Green Function solution, the values of the phases of the eigenfunctions, η n , close to the non-perturbative region, at small k ∼ Λ QCD . The fits using both relations give an excellent description of data with similar χ 2 values.
The fits lead to the unexpected result that the first eigenfunction decouples or nearly decouples. This means that the overlap of the first eigenfunction with the proton impact factor is very small or even zero, due to the fact that the first eigenfunction has a transition region from negative to positive values, i.e. a node. Therefore, the first eigenfunction chosen by the fit cannot be a ground state. This suggests, as a consequence, the existence of a multiply interacting ground state which may have properties of the soft pomeron. The contributions of such a state would be rapidly attenuated as Q 2 increases. However, at low Q 2 , it should dominate the F 2 and diffractive processes. A particularly good place to study its effects should be the exclusive diffractive vector meson production, ρ, φ and J/ψ, because in this reaction the value of the Regge slope, α , is also measured. We may try to learn more about it in our forthcoming paper by focusing the investigation on the region closer to Λ QCD , by varying k 0 and, last but not least, using the complete information concerning the errors of HERA data [13].
The present BFKL fits to HERA data predict that in the very low x region, x << 10 −4 , Q 2 > 6 GeV 2 , F 2 should grow with a slope λ which is close to the eigenvalue of the second eigenfunction and which is Q 2 independent. This prediction is possible because there is no interference between the ground state and the second eigenfunction, since they are orthogonal to each other and have very different support. The ω value of the second eigenfunction could be easily measured on some future ep collider, such as VHEeP [26] or LHeC [25].
Finally, let us note that the AB(C) fits, should be affected by supersymmetry or other physics beyond the standard model (BSM), as discussed in our previous papers [3,4]. This is because (at least at LO) the constant A is proportional to the beta function, which changes its value drastically once the threshold for the production of gluimos or other BSM particles is crossed. The decoupling of the leading eigenfunctions makes the analysis of BSM physics simpler, especially on the future VHeP or LHeC colliders. This is because the ground state is now very well constrained, the value of the ω 2 can be directly measured, and the values of the higher intercepts, ω n>2 , can be parametrised reliably.

Appendix A
We rephrase here the original derivation of the BFKL resummation given in ref. [8]. It is convenient to write If a is small then up to order a we have We may write χ 1 (ν) (defining a quantity χ reg 1 (ν) ) as By a suitable choice of the constants A and B, we can arrange for χ reg 1 (ν) to be free of singularities as ν → ± i 2 . In this limit we have So that Therefore the constants A and B are selected to match the single and double poles respectively of the function χ 1 (ν) and in that way χ reg 1 is free from such singularities. in ref, [8] it is pointed out that the correction due to χ reg 1 is genuinely negligible and the entire large correction to the characteristic function come from the terms which are singular as ν → ±i/2. Now let us consider another functionω(ν) which is defined as the solution to the transcendental (implicit) equatioñ Solving to leading order in α s we havẽ Expandingω(ν) up to order α 2 s , and using (6.3) we obtaiñ Thus we see that up to order α 2 s , the quantities ω(ν) andω(ν) are identical so that up to that accuracy we may replace the usual perturbative expression given in (3.1) byω(ν).
On the other hand, the quantityω(ν) does not contain any singularities as ν → ± i 2 . The singularities we see in eq(6.10) are only present as a result of an expansion. They are therefore an artifact of this expansion and are not present for the entire function. Since it is these singular terms that give rise to the large NLO corrections found in χ 1 (ν) we may consider the quantityω(ν) to be the expression in which all of these large corrections have been resummed.
For the case of the third order pole, this has been established exactly, since we know what the origin of the triple pole is. In rev [8] it is explained that this arises from a mismatch between the "rapidity", Y , of the forward gluon-gluon scattering amplitude used in the BFKL approach Y = ln s kk For the resummation of the double and single poles, this is not known uniquely and there are an infinit enumber of possible resummation schemes, of which one is described here, and three oithers are discusswd in ref. [8]. All these reummation schemes have in common the fact that they resum all the collinear singularities (i.e. all poles as ν → ± i 2 and they are all equivalent to the ordinary pertubative expansion for ω up to order α 2 s . They, of course, differ, in the terms proportional to α 3 s and higher -but we have no reason to select one of these schemes above another in the absence of the NNLO calculation of the characteristic function. Scheme 3, which is the scheme considered here is the most convenient for our purposes.

Appendix B
Absence of nodes in the wave function of a ground state One can define the kinetic energy, T [ψ], as Integrating by parts we obtain provided the wave function ψ(x) is continuous and has continuous first derivatives. (The transition from (7.1) to (7.2) is not valid for the continuous wave functions which do not have a fully continuous first derivatives, like e.g. ψ(x) ∼ α|x| or ψ(x) ∼ exp(−α|x|).) In the following, we prefer to use for kinetic energy the expression (7.2) since, in contrast to (7.1), it is always positive.
Let us first consider the case of the one-dimensional Schrödinger equation and define the total energy as a functional Let us assume, that the ψ-function changes its sign, for example, ψ(x)| x→0 ∼ x, and prove, that there is a positive function χ(x) with χ(0) = 0, which has a smaller energy E. It would mean, that the wave function ψ with a node at x = 0 cannot be the wave function of the ground state.
Let us turn now to the BFKL equation with the running coupling constant. In the leading logarithmic approximation we have with where α s (t) = 1 β 0 t , t = ln |k ⊥ | 2 Λ 2 QCD . (7.9) Here E = −ω plays the role of the total energy in the Schrödinger equation. The operator ν denotes the momentum canonically conjugated to the coordinate t, [ν, t] = i . (7.10) As usual in QCD, one can use the perturbative hamiltonian H for large t > t 0 > 0 only. For t < t 0 it should be substituted by an hermitian non-perturbative hamiltonian H and the corresponding wave functions and their derivatives are matched at t = t 0 .
(7.11) Note, that for the BFKL hamiltonian, which has a non-linear dependence from ν 2 , it would be natural to introduce a trial function χ with continuous higher derivatives in the points t − t 1 = ± . But in the correction to the total energy, expressed in terms of the functional E = dt f (t) H f (t) , ||f || = 1 (7.12) with the substitution f (t) → |f (t)| → χ(t), the contribution from the region |t − t 1 | > will ramain unchanged. In the region |t−t 1 | < , the higher derivatives of the BFKL hamiltonian H, acting on the simple polynomial functions χ(t) and f (t), should be neglected. Note that this corresponds to the diffusion approximation, because only terms proportional to ν 2 , in the expansion of the hamiltonian H, should be taken into account.
As above, corrections to the normalisation condition and to the running coupling factors α s (t) are small. Thus, the main correction to the total energy of the trial function can be written as when → 0. Because this correction is negative we conclude that the ground state wave function for the BFKL pomeron cannot have nodes.