The Behaviour of the Green Function for the BFKL Pomeron with Running Coupling

We analyse here in LO the physical properties of the Green function solution for the BFKL equation. We show that the solution obeys the orthonormality conditions in the physical region and fulfills the completeness requirements. The unintegrated gluon density is shown to consists of a set of few poles with parameters which could be determined by comparison with the DIS data of high precision.


Introduction
The BFKL equation is particularly well suited for description of the behaviour of gluon density in the low-x region. Its application in this region is of major importance for the LHC and cosmic ray physics. In recent years we have therefore investigated [1,2] the solution of this equation using the discrete eigenfunction method, first proposed in [3]. The method is closely connected to the Green function approach which is, in our view, the most suitable since it does not require any cutoff on the BFKL dynamics.
The results reported in these papers where as interesting as they were puzzling. The proper description of data was only achieved by involving a large number of eigenfunction, O(100), contributing in a slowly convergent way. On the other hand, the third and higher eigenfunctions were already sensitive to physics at very large scales, much beyond the Q 2 region of data. If true, this property would offer a framework for investigations of Beyond the Standard Model (BSM) physics at large energy scales using relatively moderate beam energies combined with precision measurements.
The slow convergence of the procedure given in [1,2] suggests, in particular, that the discrete eigenfunctions might not by themselves form a complete set. Therefore, in a recent paper [4] we rederived the BFKL Green function using the complete BFKL spectrum and showed how the imposition of both UV and IR boundary conditions leads to a discrete set of poles. The purpose of this paper is to investigate the properties of the full Green function solution and to relate it to our previous work.
The paper is organized as follows: In section 2 we recapitulate the results of [4] which were obtained in the semi-classical approximation. To illustrate the role of approximations used in our method we discuss in section 3 the case in which the coupling runs without thresholds. In this case it is possible to obtain the Green function in analytic form, which reduces the solution to an expression in terms of Airy functions in the diffusion approximation in which the characteristic function, χ(ν), is simplified to a quadratic function of oscillation frequency, ν. In section 4 we discuss a restriction that is imposed on possible paths for the integration over the Mellin transform variable, ω, arising from thresholds in the running of the coupling, and discuss the agreement of numerical results using two substantially different ω-paths. In section 5 we discuss the orthonormality and completeness of the BFKL eigenfunctions obtained within the approximation that we are using and report that in order to obtain the required completeness relation it is necessary to include the continuum of states for which ω is negative. In section 6 we discuss the behaviour of the residues of the poles as the gluon transverse momentum increases and show that the leading pole is attenuated so that the subleading poles acquire an ever-increasing significance as the transverse momentum increases. In section 7 we present the results for the unintegrated gluon density in a model in which a very simple ansatz is used for the proton impact factor and for the phase of the oscillations of the wavefunctions in the infrared regime. In section 7.1 we discuss the convergence of the sum over poles and show that, in contrast to the results of ref. [1,2], in the formalism that we are using here this convergence is quite rapid. Section 7.2 discusses the role of the continuum for negative ω on the calculation of the unintegrated gluon density and section 7.3 shows the relation to the DGLAP calculation. In section 7.4 we discuss the behaviour of the gluon density as a function of x at fixed low virtualities, k 2 T , and show that qualitatively it has properties similar to data and display clear contributions from subleading poles. In section 8 we show how the presence of new physics albeit at very large energies affects the running of the coupling such that both the positions and the residues of the poles are altered and this in turn gives rise to a measurable change in the unintegrated gluon density. Section 9 contains our conclusions.

Green Function of the BFKL equation
We showed [4] that the (Mellin transform of the) Green function for the BFKL equation with running coupling can be solved in the semi-classical approximation in terms of Airy functions, Ai(z) and Bi(z), where z is a function of the transverse momentum of the QCD pomeron, k T , and the Mellin transform variable, ω.

More precisely the equation for the Green function is
where t = ln(k 2 T /Λ 2 ). We note that we have introduced the running of the coupling in such a way that the hermiticity of the kernel is preserved so that its eigenfunctions form a complete orthonormal set. Eq.(2.1) has a solution, in the semi-classical approximation, of the form (2. 2) The argument z(t) of the Airy functions is given by and χ(ν) are the eigenvalues of the kernel K. The parameter t c is the (ω-dependent) value of t for which ν ω (t c ) = 0. N ω (t) is a normalization factor given by This Green function is analytic in the entire ω-plane with the exception of an essential singularity at ω = 0.
The expression on the RHS of eq.(2.2) has the desired ultraviolet behaviour, namely it is exponentially attenuated as t → ∞, but the infrared behaviour (where t is small) is not properly determined. To obtain the most general solution to (2.1) this Green function should be generalized by adding to Bi the solution of the homogeneous BFKL equation, i.e. the transformation Bi(z(t)) → Bi(z(t)) ≡ Bi(z(t)) + cot (φ(ω)) Ai(z(t))., (2.6) and eq.(2.2) becomes The transformation (2.6) plays the dual role of providing a set of poles at the values of ω for which the function φ(ω) = nπ and fixing the phase of the oscillatory behaviour of the Green function for very small t (or t ), thereby providing the connection between the determination of the infrared behaviour of the Green function and the position of the poles, originally suggested in [3].
It was pointed out in [4] that upon inverting the Mellin transform of the amplitude, by integrating along a suitable path in the ω-plane, the saddle-point approximation, used to match the BFKL solution with the result of a DGLAP analysis [5,7,8], could be valid provided the saddle-point was to the right of all the poles in the ω-plane. On the other hand, if this is not the case then the selected contour for the integral over ω must surround the poles to the right of the saddle-point and will provide significant supplementary contributions to the unintegrated gluon density which are not matched to the DGLAP result.
In this paper we report on a numerical analysis of the above-mentioned semi-classical solution and discuss the behaviour of the eigenfunctions of the kernel. We confine ourselves to the leading order BFKL kernel. The effects of the large components of the NLO BFKL kernel will be discussed in a subsequent paper.

Explicit Solution in the Absence of Thresholds
We begin by considering the simplified case in which the running coupling is given bȳ Consider the eigenfunctions, f ω (t) of the kernel with running coupling The eigenfunctions may be written in integral form as Provided g ω (ν) is square-integrable, the contour C may be taken to be along the real axis in the ν-plane and we have (3.5) These eigenfunctions obey the orthonormality relation and completeness relation For the leading order BFKL equation with running coupling, g ω is given by The integral over ν in eq.(3.3) can be approximated by a Gaussian integral around the saddle-point, ν ω (t), given by χ (ν ω (t)) =β 0 ωt (3.9) For sufficiently small t the eigenfunctions have an oscillatory behaviour In terms of the function g ω (ν), the Green function may be written (see [6]) as (3.11) where c(ω) is an arbitrary function of ω and this second term reflects the fact that one can add to a Green function any solution to the homogeneous part of the equation for the Green function. If we write c(ω) ≡ cot (φ(ω)) , (3.12) then in the semi-classical limit, we can interpret φ(ω) as being the phase of the oscillations of the BFKL eigenfunctions at some small value of t.
That (3.11) is indeed the solution to the equation for the Green function can be seen by applying the (hermitian) operator ×g ω (ν )e itν e it ν (3.14) Integrating over ν by parts and using g ω (−ν) = (g ω (ν)) −1 yields Thus we see that it is the factor (ν + ν ) inside the integration over ν and ν that generates the required inhomogeneous term in the Green-function equation.
For small t (t < t ) this Green function has an oscillatory t behaviour The phase of these oscillations fixed by the boundary conditions at some small value, t 0 , of t determines the function φ(ω) and hence the positions of the poles at ω = ω n where φ(ω) = nπ.
The Green function has a spectral representation in terms of these poles, namely and we see from the completeness relation that In the limit t → −∞ (keeping t fixed), the integral over ν may be approximated by its contributions in the regions of the two saddle points at where χ(ν ω (t)) =β 0 ωt. (3.19) In this limit, |ν ω (t )| > |ν| and so the function (ν + ν ) is replaced by ±1 at the saddle points ±ν ω (t ) respectively.
We define Performing the gaussian integrals over ν and ν around the two saddle-points we obtain Exploiting the symmetry under t ↔ t , we arrive at the semi-classical approximation for the Green funcion in the region × cos s ω (t) + π 4 + cot(φ(ω)) sin s ω (t) + π 4 + t ↔ t , (3.23) This is the approximation to the Green function given in eq.(2.7) when t, t t c . A similar argument, using the semi-classical approximation, can be used to show that the expression (2.7) matches the semi-classical approximation when t or t t c However, the semi-classical approximation breaks down if t or t ≈ t c , since ν ω (t) becomes very small and the curvature at the saddle-points becomes small. Nevertheless, as we have shown in [4], in this limit the characteristic function may be approximated by a function which is quadratic in ν such that the homogeneous part of the equation for the Green function reduces to Airy's equation and the particular combination of Airy functions given in eq.(2.7) generates the appropriate inhomogeneous term, so that the expression (2.7) is a good approximation to the Green function in all regions of t and t .
Henceforth we take the running coupling to be given by whereβ 0 (t) has steps at the heavy particle thresholds. With this more realistic function for the running coupling, it is no longer possible to solve the Green function analytically, even in integral form, and we confine ourselves to a numerical analysis within the semi-classical approximation for which the integral over the variable ν has a saddle-point at ν ω (t).

Numerical Solution Using Different Paths on the ω-Plane
The unintegrated gluon densityġ(x, t) is given by the inverse Mellin transform of the Green function byġ where Φ P (t) is the proton impact factor. We take the Green function to be given by eq.(2.7) and the running coupling given by eq.(3.24). The function φ(ω) is given by with t c given by the relation .
This means that at some small value, t 0 , of t the phase of the eigenfunction with eigenvalue ω is given by a non-perturbative phase, η N P . 1 For a numerical evaluation of the integral over ω, it would be most efficient to identify the saddle-point of the integrand and select a path for ω which passes through that saddle point and follow the path of steepest descent. In the case where the saddle point lies to the left of any of the identified poles of the Green-function the integral must be supplemented by the integral over a contour surrounding all poles to the right of the saddle-point. 2 .
Unfortunately, there are restrictions on the permitted paths in the case where the running coupling encounters thresholds. In order to consider complex values of ω, we require complex values of α s (t) and hence complex values of t. In the presence of fermion masses m i , the running coupling to leading order is given in terms of its measured value at t = t by [ where the function F , given by 1 η N P can be a function of ω but in this paper we take it to be constant, although for a realistic fit to data we would expect it to possess some ω dependence. Note that φ(ω) does not depend on the choice of the infrared scale, t 0 , but the infrared phase η N P does. 2 A full discussion of this is found in [4] is multi-valued in the complex plane. The running coupling is therefore only uniquely specified in terms of t corresponding to k T covering the complex plane once so that the imaginary part of t is restricted to −π ≤ m{t} < π, which restricts the imaginary part ofᾱ s and consequently the imaginary part of ω. In particular, the calculation of the argument of the Airy functions requires the identification of α s (t c ), so that the restriction on the range of the imaginary part ofᾱ s leads to a corresponding restriction on the imaginary part of ω. If the real part of ω is small then the real part of α(t c ) is small, i.e. the real part of t c is large. The restriction on the allowed range of the imaginary part of t c therefore implies that the imaginary part of ω must also be small -i.e. we need to select paths which are very close to the real axis in this region. Furthermore, the restriction on the imaginary part pushes the possible paths towards the essential singularity at very small ω so that we could perform the integration only to some small value of ω min of around 0.05.
We have selected several paths whose imaginary part differ substantially for large ω. Two of them are shown as an example in Fig.1. We have performed the integral of eq.(4.1) along all of these paths and find negligible difference over a large range of t for t = 1 to t = 17 (corresponding to transverse momentum k T ≈ 2 TeV). 3 In principle, this result should be expected, since there are no singularities of the Green function off the real axis and therefore one can deform the contour into any contour that surrounds the real axis and crosses the real axis to the right of all poles. However, our contours are not fully closed although the results, at lower t values 4 , were independent of the ω min value. This means that the missing piece of the paths gave a negligible contribution, in this t region. The independence of the results from ω min was a first sign that the full Green function of eq.(2.1) and eq.(2.6) behaves differently from the slowly converging sum of eigenfunctions constructed in ref. [1,2]. In Section 7.1 we will explain in details why this happens.
In addition, the good agreement of the integrals over the different paths shows the numerical consistency of our Mellin transform integration.

Properties of the Eigenfunctions
The Green function of the BFKL equation has poles at ω = ω n where φ(ω n ) = nπ. If the infrared phase η N P is set to a constant (ω-independent) value of −π/4, the first few eigenvalues are shown in Table 1.   Except for the first three of these, the eigenvalues are well approximated by , consistent with the estimate found in [3] and [1].
These correspond to a discrete set of eigenvalues of the BFKL kernel with running coupling whose normalized eigenfunctions are given by These functions look superficially like the eigenfunctions of ref. [1,2]; the Airy functions were defined exactly in the standard way and the normalization coefficient was previously determined numerically (from the requirement that every eigenfunction should be normalized to unity). Here the normalization coefficient, N ωn , is known analytically and is given by eq.(2.5); the first factor in eq.(5.2) arises from the conversion between the continuum and discrete eigenfunctions. The continuum ones are normalized to a δ-function in ω whereas the discrete ones are normalized to a Kronecker δ-function in the eigenfunction number, n, This conversion factor can be obtained from the relation where ∆ω n is the separation of the n th and (n + 1) th eigenvalues. For large n the eigenvalues are closely packed and the separation to the eigenvalues may be approximated by The eigenfunctions start by oscillating for small t with n turning points for t < t c and are evanescent for t > t c . The first three such eigenfunctions are plotted in Fig.2  In Table 2 we show the numerical results for the orthonormality relation (5.3) evaluated for the first 7 eigenfunctions of eq.(5.2).
We see that the use of the semi-classical approximation has only had a very small effect on the orthonormality relation (5.3). Furthermore, in order to preserve the validity of the perturbative expansion, we cannot integrate over all values of t, but only for t > t 0 . However, the eigenfunctions are small at small values of t. To see this we note that the additional semi-classical factor, N ω (t), given by eq.(2.5) is constant for t ≈ t c , where the Airy functions alone is a good approximation to the eigenfunction, but for sufficiently small t, the variable z and χ (ν ω (t)) both become very insensitive to t so that N ω (t) then has a t dependence ∼ 1/ ᾱ s (t) -i.e. it decreases in the infrared region asᾱ s (t) grows. The t dependence of these normalization factors for the first three eigenfunctions is shown in Fig.3.
From Table 2 we find that the orthonormality condition is obeyed to a very high accuracy. We consider this as a strong indication that we are using a consistent set of physical approximations; in particular the semi-classical approximation and the running of α s are preserving the hermiticity of the Hamiltonian.  For the completeness condition we consider the sum For k T = 100 GeV we plot this sum in Fig.4 for the N = 5, 10 and 20 eigenfunctions.
We see that the sum has converged 5 after 10 -20 eigenfunctions to a distribution on k T which is peaked at k T = k T . However, we note that the distribution is very broad. This tells us that the discrete eigenfunctions do not form a complete set by themselves. Rather the completeness requires that the sum over the discrete eigenfunctions must be supplemented by the integral over the continuum of states for which ω takes negative values.
For negative ω there is no critical transverse momentum, t c , beyond which the eigenfunctions diminish, but have oscillatory behaviour for all t. We can write the negative omega eigenfunctions as Since the ultraviolet boundary condition does not impose a specific ultraviolet phase, the spectrum is continuous. 6 Although there is no ultraviolet phase-fixing condition, there can be an infrared boundary condition which determines the phase of the oscillations at small t. For small positive ω the eigenfunctions are very closely spaced and become indistinguishable from a continuum. For small negative ω, the non-perturbative phase should match its value for small positive ω in order to ensure a smooth function as ω changes sign. Note that for large negative ω it may be the case that the infrared phase is not defined. An example of a mechanism in which this happens is where the infrared behaviour of QCD is simulated by an effective gluon mass [11]. Here it is found that at some negative ω = −ω 1 , there is a phase transition below which the infrared phase is not determined. Other possible sources of such phase transitions could arise from the restoration of conformal invariance at some high-energy scale. 5 the fast convergence of this sum will be discussed below 6 This is analogous to the fact that particle only form bound-states for negative energy. Here the analogue of energy is −ω. For sufficiently large (negative) values of ω, ν ω (t) is given by and from which we can see that for large t these negative ω eigenfunctions have oscillations whose amplitude and frequency increase rapidly as t increases. An example for which ω = −1 is shown in Fig. 5. With the inclusion of the continuum states, the completeness relation becomes In Fig. 6 we plot this quantity for ω min = −1 and ω min = −2 for t = 10. We see that as ω min becomes more negative the LHS of eq.(5.8) becomes more sharply peaked -tending to the required δ-function in the asymptotic limit ω min → −∞.

Transverse Momentum Dependence of the Residues of the Poles
As t increases from t = t 0 , (where the infrared phase is set), the eigenfunction f n oscillates through (n − 1) nodes before the value of t c (n) for that eigenfunction, with (in leading order) t c (n) = 4 ln 2 ω n (6.1) whereas for values of t > t c (n) it decays exponentially. This means that for small t, t ≤ t c (1) the Green function as a function of ω has a series of poles at ω = ω n with residues that oscillate with amplitudes that decrease with increasing n, reflecting a convergence of the sum over pole contributions. This is shown in the left-hand graph of Fig.7, where we have taken t = 2 (and ω taken close to the real axis).
In the region t c (1) < t < t c (2), the residue of the leading pole is attenuated and the leading behaviour is now given by the first sub-leading pole, which still has an oscillatory residue. This is shown in the right-hand graph of Fig.7 for t = 12. We note that for this value of t the residue of the leading pole has significantly diminished. Note also that the sign of the residue of the sub-leading pole is opposite to the case where t = 2, reflecting the oscillatory nature of the residues of the poles for t < t c (n).
As t increases further, the residues of more and more of the sub-leading poles start to decay and the inverse Mellin transform of the Green-function is dominated by the contribution from smaller and smaller values of ω n . This is in sharp contrast to the situation in which the coupling is kept fixed and for which the inverse Mellin (for large rapidity) transform is always dominated by the region close to 4ᾱ s ln (2). For the running coupling this is not the case, but the value of t at which the particular sub-leading poles dominates depend on the value of the rapidity, Y (or ln(1/x) in the case of deep-inelastic scattering, leading to an effective pomeron intercept which depends on Y .

Properties of the Unintegrated Gluon Density
The Green function also has a spectral representation given by eq.(3.17) and so we should be able to obtain a good approximation to the unintegrated gluon density by summing over the pole contributions from n = 1 to n = n max .
As a qualitative demonstration of the unintegrated gluon density that can be obtained from the BFKL Green function, we consider a very simple model in which the non-perturbative phase, η N P is set to the constant value 7 of − π 4 . The proton impact factor has to be positive everywhere and concentrated at the values of k T < O(1) GeV. It is usually assumed to be of the form as discussed in ref. [1,2]. The form (7.1) vanishes as k 2 T for small k T , as required by colour transparency and the coefficient b has the interpretation of the average inverse square transverse momentum of partons inside the proton and is therefore of the order of 10 GeV −2 . The overlap integral of the proton impact factor with the eigenfunction for t > t 0 (t 0 corresponds to k T = 1 GeV) is therefore determined by the value of the amplitude at t = t 0 . This is due to the fact that the proton impact factor falls for t > t 0 at a rate which is much faster than the oscillation frequency of the eigenfunctions in the region t ∼ t 0 . In ref [1,2] also other forms of the proton impact factor were investigated, e.g. with different powers of k 2 in the prefactor and/or the exponent. It was found, however, that the fit to data has no sensitivity to such alternatives due, again, to a small range of the impact factor in comparison with the rate of change of the eigenfunction amplitudes. Therefore, for the purpose of this paper it is sufficient to take an impact factor which have support only at t = t 0 . With these parameters we plot the unintegrated gluon density by performing the inverse Mellin transform given by eq.(4.1) over one of the contours shown in Fig. 1 and compare it with the result obtained from the summation over the first 7 poles. We plot this in Fig. 8 for two different values of x, namely x = 10 −2 and x = 10 −3 . The pole result was increased by a factor of 1.03 for visibility. Without this increase both results would be indistinguishable. This perfect agreement was obtained because the path enclosed all the poles used in the sum. The integration over the path of Fig. 1 was performed down to ω min = 0.065, which is between the ω min values of the 7th and 8th pole. We checked this agreement for other values of ω min and the corresponding sum of poles and obtained an equally good agreement.
This perfect agreement is, of course, due to Cauchy's theorem. This agreement is however non-trivial because the path integration is not closed what means that the missing piece of the path gives a negligible contribution, which is a first sign of a very good convergence of both the path and pole computations. In addition we note that both computations are numerically very different so it is therefore a very good check of the computational accuracy.
Another reason for this very good agreement is the fact that we limited the comparison to the region of relatively small virtualities, k T < 100 GeV. At x = 10 −2 , this region is experimentally relatively well accessible at LHC and possible future colliders. At x = 10 −3 and for smaller x values the experimentally accessible k T region diminishes substantially. Notwithstanding this, the subleading poles turn out to have a measurable effect (as we shall discuss in subsection 7.4) and this will be significant for the prospect, discussed in section 8 of the identification of new physics from the shift in positions and residues of sub-leading poles.
In the low k T region, as first indicated in Fig. 8, we have a fast convergence of the gluon density as a function of the number of poles used in the summation. In Fig. 9 we show the gluon density computed with different number of poles. The dashed-dotted line shows the leading pole contribution, the dashed line shows the sum of 5 poles, the solid line the sum of 10 poles and the dotted line the sum over 15 poles. We see that x = 10 −2 the first 10 poles almost converges in the whole k T region, for x = 10 −3 the 10 and 15 pole summation are completely indistinguishable. We note also that for transverse momenta below around 10 GeV for the case x = 10 −2 and around 20 GeV for smaller values of x, the leading pole provides a reasonable approximation.

Convergence of the Sum Over Poles
We can understand how the pole sum converges by considering the behaviour of the normalized eigenfunctions for large n. In this region the eigenvalues are very small (ω n ∼ 1/n) and the critical momentum t c is very large (proportional to n). This means that for accessible values of t the RHS of eq.(2.4) is very small and eigenfunctions oscillate with approximately a fixed frequency, ν 0 , given by χ(ν 0 ) = 0.
In this region we may use eq.(3.16) and the fact that the phase of the oscillation is η N P at t = t 0 to write φ(ω) as The normalization factor N ω (t) given by eq.(2.5) has a numerator factor |z(t)| 1/4 which cancels an identical factor in the denominator of the Airy function Ai for t t c and the remaining factor is approximately n independent as we replace ν ω (t) by ν 0 in the argument of χ .
The upshot of this is that the two factors of 1/ φ (ω) (see eq.(5.2)) give rise to a convergence of the eigenfunction series at small t (relative to t c ) like ∼ 1/n 2 . Since t c increases quickly with the eigenfunction number n this fast convergence always happens in the experimentally accessible region of t.

The Role of the Continuum for Negative ω
We have seen in section 5 that the contributions from eigenfunctions with negative eigenvalues (i.e. negative ω) are essential in order to provide a complete set of eigenfunctions obeying the closure relation (5.8). In this section we discuss the effect that the Green function with negative ω has on the unintegrated gluon density.
The contribution to this expression for the unintegrated gluon density from the (positive ω) poles vanishes asymptotically with t, reflecting the behaviour described in section 6 whereby as t increases the residues of more and more of the poles pass from an oscillatory behaviour to an exponentially damped behaviour.
The contribution, ∆ ω<0 to g(x, k) from negative ω is given by with f −|ω| (t) given by eq.(5.5) (and t = 2 ln(k T /Λ)). At first sight, one would expect this contribution to have a negligible effect on the unintegrated gluon density owing to the factor of x ω . However, as can be seen from Fig. 10, for values of ω just below ω = 0 the integrand of the RHS of eq.(7.4) is still quite large and rapidly oscillating, although we can also see that the integral converges to a fairly small value by ω = −1.
We have computed the gluon density including the contribution from negative ω. The contribution of positive ω's was given by the sum of the first 10 poles and the contribution between ω = 0 and ω 10 (where the poles are densely packed) was treated as a continuum. The negative omega contribution was evaluated from ω = −2 to ω = 0 using eq. (7.4). Fig.  11 shows the unintegrated gluon density for x = 10 −2 , x = 10 −3 and x = 10 −4 including the contribution from negative ω. The dashed line shows the pole contribution (computed from the sum of 10 poles), the solid line shows the same pole gluon density with added contribution of negative ω's. The negative ω contribution was computed assuming that the infrared phase for positive and negative ω match at ω = 0. Fig. 11 shows that the contribution of negative ω may be substantial at x = 10 −2 if the negative ω infrared phase is fixed and substantially different from 0. However, this contribution diminishes very fast with decreasing x as can be seen in the same figure.

Comparison with DGLAP
We see from Fig.8 that for x = 10 −2 above k T of around 30 GeV, the unintegrated gluon density ceases to rise, whereas for smaller values of x, e.g. x = 10 −3 , the unintegrated gluon density continues to rise up to transverse momenta of 100 GeV. This can be understood from the t-dependence of the residues of the poles as seen for example in Fig.7. For t = 12 (corresponding to k T of around 100 GeV), we see that the residue of the sub-leading pole is of opposite sign and slightly larger in magnitude that that of the leading pole (whose residue is evanescent since t is above t c for that pole). The contribution of the sub-leading pole is suppressed by a factor of For x = 10 −2 the contribution from the sub-leading pole is sufficient at such values of t to halt the rise in the unintegrated gluon density, whereas for x = 10 −3 it is insufficient. Nevertheless, at sufficiently large t the unintegrated gluon density for x = 10 −3 will also cease to rise and will eventually display oscillations. Above some large t, outside the range of t ( k T ) range considered in this paper, these oscillation are certainly unphysical because they lead to negative gluon density.
There is no reason a priori why the BFKL amplitude should not display oscillations. The inversion of the Mellin transform consists of an integral over ω which has greatest support at the saddle-point ω s . For values of t below t c for this value of ω the amplitude displays oscillatory behaviour and it is only when t exceeds this critical value that the oscillations halt. The unphysical oscillatory behaviour indicates that the solution to the BFKL equation is being applied to deep-inelastic scattering outside the kinematic range for which it was intended. The application of the gluon scattering amplitude in the Regge regime to the determination of the unintegrated gluon density identifies the rapidity with ln(1/x), which in LO is only valid provided the rapidity significantly exceeds (t − t ). Therefore for x = 10 −2 and t confined to the low t region where the proton impact factor has non-negligible support, we would expect the BFKL determination of the unintegrated gluon density to become invalid if t is substantially larger than ∼ 5 (corresponding to k T of order 10 GeV).
We would expect this limitation on the allowed range of t to be less stringent when NLO effects are taken into account. Indeed, as pointed out by Salam [13], the LO treatment ignores the discrepancy between the rapidity variable used in a BFKL analysis (which is symmetric in t and t ) and the variable x. This introduces a factor of e ω(t−t )/2 , whose absorption generates the largest part of the NLO contribution to the characteristic function, χ(ν). We would therefore expect the BFKL amplitude computed at NLO to be less sensitive to the difference (t − t ) than a purely LO analysis.
It is known that at sufficiently large t and sufficient small x, the double logarithm limit (DLL) of a BFKL analysis matches that of a DGLAP analysis. In Mellin space, the region in which this approximation is valid is given by In terms of x this translates into the limits α s (t) ln 1 x 1, ln 1 x 1 Moreover the match between a DGLAP analysis and a BFKL analysis can only be valid if t exceeds t c at ω ≈ ω s , where ω s is the saddle-point for the inversion of the Mellin transform, i.e. the region around ω s is where the integrand has its maximum support (assuming that this saddle-point lies to the right of all poles). We have seen that we need to have values of x smaller than x = 10 −3 in order to avoid a signal from the oscillatory part of the BFKL eigenfunctions. This means that in the case of a DGLAP analysis, the DLL limit can only be reached for very small values ofᾱ s (t), i.e. values of t way beyond any reasonable experimentally accessible region. For the case of BFKL with running coupling we need to go to even smaller values of x before the DLL becomes a reasonable approximation. This is because for running coupling the contributions from leading poles are attenuated at large t and we need to be at sufficently large rapidity to ensure that these leading poles dominate the unintegrated gluon density. The input for a DGLAP fit is the structure function at some reference photon invariant mass, Q 0 . In the case of the discrete BFKL pomeron the input would be the proton impact factor and the infrared phase η N P . As discussed in section 7 this impact factor is expected to be dominated by its value at t = t 0 so that the only free parameter associated with this impact factor would be the overall normalization. The only other parameter which can be substantially varied is the infrared phase, η N P ,which should be a function of the eigenvalue ω n . The infrared phases η n are generated, as the eigenvalues ω n , by the quasi-bound states of gluons inside the proton and therefore should be described by a simple parameterization. Our previous experience indicate that two parameters may be sufficient in order to generate the η − ω dependence required to fit data.
Because of these very different parameterizations, it is quite likely that there exists an overlap kinematic region at low-x for which data can be equally well described by a conventional DGLAP fit or a fit to the discrete BFKL pomeron. A detailed comparison of the discrete BFKL pomeron with data will be performed in the next paper, after the NLO corrections are introduced. In the same time we will discuss the comparison with the DGLAP fit. It is well known from HERA data that the x dependence of F 2 , which is directly connected to the gluon density, exhibits a striking Regge type behaviour, i.e. ∼ (1/x) λ . The parameter λ is not a constant, it increases logarithmically with Q 2 (which we set here equal to k 2 T ), see e.g. [14]. Such a behaviour is also a feature of the gluon density obtained from the Green function solution of BFKL investigated here, see Fig. 12. The figure shows the unintegrated gluon density as a function of x, determined from the pole contribution only, for various values of k 2 T . In this log log plot the function (1/x) λ is a straight line, so we see immediately that the gluon densities exhibit the same striking linearity as data. The slope λ increases slightly with k 2 T which can be seen by comparison of the slope of gluon density (full line) with the k 2 T -independent slope of the leading pole contribution (dashed line). The leading pole contribution is, for each k 2 T , normalized to the values of the gluon density at x = 10 −2 . In Fig. 12 we display the gluon density in the validity region of our solution to BFKL, ∆t < log(1/x), which means that k T should be smaller than order of 10 GeV at x = 10 −2 . At larger k T , not shown in this figure, we observe a clear deviation from linearity in the region of x between 10 −2 and 10 −3 . A close inspection of the highest k 2 T line in Fig. 12 shows that this effect sets in already at a relatively low k T of about 30 GeV.
In Fig. 13 we show the x dependence of the gluon density which includes the negative ω contribution (full line) for various k 2 T 's. The dashed line shows for comparison the same leading pole contribution as in Fig. 12. We see clearly from this plot that the negative ω contribution substantially affects the linearity in the region of x > 10 −3 . Since this nonnegligible contribution from negative ω occurs at the relatively low values of k 2 T considered here whereas data strongly indicates linear trajectories in x up to x ∼ 10 −2 , we take this as a strong hint that the infrared phase, η N P for ω < 0 and the precise form of the proton impact factor are such that the overlap of the proton impact factor with the negative ω part of the Green-function is very small. In view of the fact that the proton impact factor is expected only to have significant support near t 0 where the infrared phase, η N P , is fixed, such suppression of the overlap would occur if the infrared phase were very small in the region ω ∼ 0.
In this paper we do not attempt to make any detailed comparison with data because we are working here in the LO only and it is well known that the NLO and LO results differ substantially in BFKL. We note, however, that the difference between the full gluon density and the leading pole contribution seen in Fig. 12 is due to the subleading poles. A similar change of slope λ with Q 2 8 , even more prominent then in our LO calculation, is also seen in the data [14]. Therefore it is highly possible that the properties of the subleading poles can be well determined from the measurements of the slope λ since in our solution there is only one free parameter per pole, η N P , and the number of contributing poles is small, due to the fast convergence of their sum.

The Effect of New Physics
In previous publications [1,2], we have pointed out that both the positions of the poles in ω and their residues are sensitive to any new physics which affects the running of the coupling and hence the value of the critical momentum, t c , provided t c is above the threshold for such new physics. This unique effect is due to the fact that an eigenfunction of the BFKL kernel, (5.2), is a (quasi) bound state of gluons with very different virtualities, ranging from t 0 to t c . Even if such a state is probed at low t value, much below the Beyond Standard Model (BSM) threshold, the result is sensitive to the properties of the whole state since the eigenvalue and the residue at the probed t is determined by all the gluons of the state 9 .
Formally speaking the sensitivity to BSM thresholds emerge from the fact that the function φ(ω), eq. (4.2), which defines the eigenvalues ω n by the requirement φ(ω) = nπ, contains an integral over the frequency ν ω (t) which ranges from t 0 to t c . This frequency, defined by eq. (3.19), is strongly sensitive to a supersymmetry (SUSY) threshold because the value of β 0 changes substantially in the SUSY region.
In Fig. 14 we show an example of this in which we plot the Green function as a function of ω on a path close to the real axis, for a typical low t values of 9, 6.7 and 4.4 corresponding to k T of 30, 10 and 3 GeV. We are comparing the Standard Model (SM) with the MSSM with a supersymmetry (SUSY) threshold of 3 TeV. We see that the position of the first pole is unaffected because the corresponding t c lies much below the SUSY threshold of 3 TeV. However the positions of the sub-leading poles are shifted to the right because their t c 's are either close to the 3 TeV threshold, as in the case of the second pole, or much above it for the rest. It is interesting to observe that the residues of the non-leading poles oscillates strongly in the displayed t region. This may help to disentangle their contribution since this t region is well accessible to high precision measurements.

Conclusions and Outlook
We have investigated here the properties of the complete Green function solution to BFKL equation in LO approximation, using a mixed technique of analytic and numerical analyses. We have shown that this solution fulfills the completeness requirement and leads to a set of eigenfunctions which are properly normalized and are orthogonal to each other. These We have taken t = 9, 6.7, 4.4, corresponding to k T of 30, 10, 3 GeV, and t = 2.
mathematical properties are fulfilled with a high numerical precision, which is not trivial in view of the fact that the eigenfunction states are defined only for t values above some small, but perturbative t 0 , and not in the whole region −∞ < t < ∞, as would be mathematically required. To achieve the completeness it is mandatory to take into account the contribution from the states of the negative ω continuum.
The unintegrated gluon density is defined, in this paper, as the inverse Mellin transform for the Green function of the BFKL equation over a suitable path in the complex ω plane. We show that at low t the integration over an ω path is exactly equivalent to the sum of poles supplemented by a possible contribution of the negative ω continuum. The sum of poles is dominated by a relatively few terms, which is in contrast to the situation found in the discrete eigenfunction solution of [1][2][3]. The fast convergence of the complete Green function solution is due to the presence of the t dependent normalization factor which was missing in the pure eigenfunction approach of [1][2][3].
We have investigated the region of validity of the unintegrated gluon density obtained in this paper and found that it is limited to the region ∆t = t − t < log(1/x). If ∆t is sizably larger the unintegrated gluon density starts to exhibit an oscillatory behaviour, which eventually leads to an unphysical, i.e. negative gluon density. Since in DIS t is limited by the proton factor, which confines it to very small values, the t values cannot be too large. At x = 10 −2 they corresponds to k T values of the order of 10 GeV. At smaller x, like x = 10 −3 , the virtuality k T would be an order or more of magnitude larger, however at HERA or LHC the region of accessible virtuality decreases with x substantially, so that the region of applicability of the BFKL equation is limited effectively to k T of the order of 10 GeV.
At low k T the gluon density is dominated by the leading pole, which leads to a linear behaviour in the logarithmic x dependence with a slope which varies with k T . This variation is due to the contribution of the subleading poles. We show here that the ω values of the subleading poles are sensitive to the Beyond Standard Model (BSM) effects due to the same mechanism of threshold sensitivity as investigated in [1,2]. Therefore, the deviations from the leading pole behaviour are sensitive to the BSM effects and could be measurable due to the high precision of the measured slopes of the logarithmic x distribution, especially in future projects, see e.g. [15,16]. However, in this paper we did not attempt to perform any data analysis because we are working here in the LO order approximation only and it is well known that the ω values differ substantially between LO and NLO approximation of BFKL. The LO analysis presented here allows the full understanding of the qualitative properties of the BFKL solution in the low k T region which we expect to be the same as in the NLO analysis. The quantitative results in the NLO approximation will be presented in our next paper in which we will also perform a comparison with the high precision DIS data.