text:
# Discrete kink dynamics in hydrogen-bonded chains I: The one-component model

## I Introduction

Hydrogen bonds (H-bonds) play a crucial role in the structure and the dynamics in a whole variety of systems ranging from ferroelectrics to biomolecules. They are of central importance in biology, when reactions are considered at molecular level. In bioenergetics, they appear even more crucial because they enable transfers of protons from one molecule to another one in networks or chains formed via hydrogen bonding .

More specifically, H-bonds  are interactions linking two molecules or ions, for example, O, N, F, and Cl atoms, or in general any pair of hydroxyl groups, which may be denoted by X, via a hydrogen ion (proton) H+, forming a hydrogen-bonded (HB) bridge X–H$`\mathrm{}`$X as shown schematically in Fig. 1. The ion to which the proton in this bridge (H-bond) is more tightly linked is called the hydrogen donor, whereas the other ion is the hydrogen acceptor. More precisely, the proton is coupled to each X- ion through a pair ion-proton interaction potential of the standard type (Morse, Lennard-Jones, etc.) with an equilibrium distance $`r_0`$, which necessarily has a finite dissociation energy as the X$`\mathrm{}`$H distance tends to infinity.

As usual, the total potential for the HB proton is of a double-well shape, but this can occur only if the motion of the heavy ions along the H-bond is appropriately constrained, so that either (i) a sufficiently strong interaction (repulsion) between the X- ions that does not allow the ions to get closer each other a distance less or equal to $`2r_0`$ or (ii) a periodic substrate potential with period exceeding $`2r_0`$ is additionally involved. In this way, a one-dimensional network of hydrogen bonds can be formed as a diatomic chain of alternating heavy (ion) and light (proton) particles coupled nonlinearly (e.g., via a Morse-type potential, like in Ref. ), whereas the second-neighbor (ion-ion and proton-proton) interactions are involved in the harmonic approximation. Under certain conditions on the ion-ion coupling discussed in Ref. , the proton in each H-bond of the HB chain can be found in two equilibrium positions separated by a potential barrier, so that the two degenerate ground states of the chain $`\mathrm{}`$ X–H $`\mathrm{}`$ X–H $`\mathrm{}`$ X–H $`\mathrm{}`$ and $`\mathrm{}`$ H–X $`\mathrm{}`$ H–X $`\mathrm{}`$ H–X $`\mathrm{}`$ are assumed to exist. Another important property, more specific for biological systems, is that the height of the potential barrier crucially depends on the distance between adjacent X- ions. Using these properties as main features of HB chains, a number of one-dimensional two-sublattice models, whose dynamic behavior is governed by the soliton theory , has been suggested and studied extensively. These soliton-like theories are based on the well-known cooperativity  of the hydrogen bonding, simply defined through the coupling of protons in the nearest-neighbor hydrogen bridges of the chain.

Since the HB chain is a diatomic lattice, the mechanism of hydrogen bonding involves two types of particle displacements. Let $`Q_n`$ and $`q_n`$ be the displacements of the heavy ion and the proton in the $`n`$th unit cell of the lattice from their equilibrium positions, at which one of the two ground states of the chain is realized, respectively. These displacements are labeled according to the sequence $`\{\mathrm{},Q_{n1},q_{n1},Q_n,q_n,Q_{n+1},q_{n+1},\mathrm{}\}`$. Then the general and the most simple model for the proton transfers in such a diatomic chain can be given through the two-sublattice Hamiltonian, consisting of two parts :

$$H=H_0+H_{ion}.$$
(1)
The first part

$$H_0=\underset{n}{}\left[\frac{m_p}{2}\dot{q}_n^2+\frac{K_p}{2}(q_{n+1}q_n)^2+\epsilon _0V(u_n,\rho _n)\right],$$
(2)
with

$$u_n=q_n\frac{1}{2}(Q_n+Q_{n+1})\text{and}\rho _n=Q_{n+1}Q_n,$$
(3)
describes the proton kinetic energy, the nearest-neighbor proton-proton interaction, and the intrabond proton energy that depends on the displacements of the protons from the midpoints in the H-bonds $`u_n`$’s and the relative distances between the nearest-neighbor ions $`\rho _n`$’s, whereas the second (pure heavy-ion) part

$$H_{ion}=\underset{n}{}\left[\frac{1}{2}M\dot{Q}_n^2+\frac{1}{2}K_{ion}\rho _n^2+\frac{1}{2}K_{sub}Q_n^2\right],$$
(4)
describes the kinetic energy of the X- ions, the coupling energy between the nearest-neighbor ions, and the interaction energy of the X- ions with a possible substrate (e.g., formed by the walls of a pore crossing a membrane). Here the overdot denotes differentiation on time $`t`$. The proton and ion masses are denoted by $`m_p`$ and $`M`$, respectively. Similarly, $`K_p`$, $`K_{ion}`$, and $`K_{sub}`$ stand for the stiffness constants of the interaction between the nearest-neighbor protons, the nearest-neighbor ions, and the chain ions with the substrate, respectively. It is important that the intrabond proton energy is given in terms of a general double-well potential $`V(u,\rho )`$ as a function of two variables: $`u`$, the proton displacement from the middle of the hydrogen bond, and $`\rho `$, the relative ion displacement. If additionally this (dimensionless) function is normalized according to the relations $`V(0,0)=1`$ and $`V(\pm a,0)=0`$, where $`\pm a`$ are the positions of minima in this function, then $`\epsilon _0`$ is the barrier height of the proton potential in the H-bond. When the heavy ions are displaced from their equilibria, this potential is deformed, with its barrier top moved together with the ions.

There have been numerous studies  of soliton solutions to the equations of motion governed by the Hamiltonian (1)-(4) including also one-component models, where the heavy ions X- are assumed to be fixed . All these studies refer to the continuum limit, which presumes the existence of a sufficiently effective cooperativity of hydrogen bonding or, in other words, the inter-bond proton-proton coupling $`K_p`$ is required to be strong enough. However, according to the ab initio calculations of the proton-proton interaction in realistic HB chains by Godzik , $`K_p41`$ Kcal/mol Å2. This magnitude appears not to be sufficient for a free propagation of the ionic defects along the HB chain with realistic values of the potential barrier height $`\epsilon _0`$. The reduction of this barrier on the basis of the two-component modelling was also shown to be not enough to provide a free soliton regime and the ionic defects in HB chains with realistic parameter values  were shown to be very narrow objects .

On the other hand, a rapidly increasing number of publications over last years (since the pioneering work of Peyrard and Kruskal ) have demonstrated significant differences in the behavior of soliton solutions treated in the continuum limit and their spatially discrete relatives . The discrete versions of the partial differential equations brings about a number of critically important modifications to the dynamics. The moving kinks of the continuum theories become propagating structures that decelerate by emitting radiation as they traverse the lattice sites. This ultimately brakes the structures and brings them to rest, or “pins” them.

The above results has been obtained for the conventional models such as the discrete sine-Gordon and $`\varphi ^4`$ chains. For more general class of models, some interesting and intriguing results has been obtained. Thus, it has been shown that the shape of the on-site (in our case, the intrabond proton energy) potential is a factor of particular importance for modelling soliton motion in physical systems. To study this effect, Peyrard and Remoissenet  have introduced a modified sine-Gordon system, where the shape of the on-site potential differs sufficiently from the sin-function. They found that if the barrier between the potential wells is flat enough, the Peirels-Nabarro (PN) barrier does not decrease monotonically with the coupling constant, as in the ordinary discrete sine-Gordon chain. It decreases with oscillations, so that the PN barrier experiences dips, where it lowers by the order of magnitude. Later , it was found that if the PN barrier decreases nonmonotonically, there exist certain velocities, at which even very discrete kinks propagate with constant shape and velocity. Everywhere in between these velocities, there exist kinks with oscillatory asymptotics (nanopterons). Approaching the problem from another side, Schmidt in  has constructed a Klein-Gordon model that allows an exact moving kink solution of the form $`\mathrm{tanh}(nvt)`$ for some specific value of velocity $`v`$. Furthermore, Flach and coauthors  have shown that for this model the PN barrier is nonzero (when $`v0`$). It was also shown that kinks of the discrete sine-Gordon equation with topological charges greater then one exhibit some features, similar to those, described above, including free propagation at some selected velocities (see Refs. ). Note that if we step out from the Klein-Gordon class of discrete models (for instance, by introducing anharmonicity into the interparticle interaction), some new phenomena related to the kink mobility can appear, but this is out of the scope of the present paper. In this context, a few papers (see Refs.  and others therein) should also be mentioned.

Thus, owing to the importance of discreteness effects in the kink dynamics, it would be of big interest to apply these findings first for the one-component model of proton transport in HB chains and then for the two-component model given by the Hamiltonian (1)-(4). In this context, as found by Duan and Scheiner , a pair of Morse functions, placed tail-to-tail so as to allow for the approach of the proton towards the acceptor while it is departing from the donor (see Fig. 1), provides the best framework for reproducing their potentials obtained from ab initio calculations. It is important that the Morse-type functions contain parameters with clear physical meaning, which vary little from one H-bond to the next one in HB systems.

The aim of this paper is to investigate the properties of the one-dimensional Klein-Gordon chain with the on-site (intrabond) potential of the double-Morse type. We are going to find both stationary and dynamic (moving) kink solutions and to show that kinks can be mobile even if being very narrow.

The paper is organized as follows. In the next section, we present the Hamiltonian and the equations of motion for the one-component model. In Sec. III, we study the properties of the stationary kink solutions. In Sec. IV, the Peierls-Nabarro potential for the kinks is investigated. Section V is devoted to the studies of kink mobility. Conclusions are given in Sec.VI.

## II The double-Morse proton potential

In the limit, when the heavy ions are fixed (immobile) at a same distance $`l`$ forming a uniform lattice, we deal only with the first part of the Hamiltonian (2), where $`q_n=u_n`$. In what follows we adopt the dimensionless description, where for the dimensionless proton displacement $`u_n/l`$ we keep the same notation $`u_n`$ and use the time unit $`t_0=l/\sqrt{\epsilon _0/m_p}`$. In these dimensionless variables, the Hamiltonian (2) reads

$$=\underset{n}{}\left[\frac{1}{2}\dot{u}_n^2+\frac{\kappa }{2}(u_{n+1}u_n)^2+V(u_n)\right].$$
(5)
Here and in what follows, the overdot denotes the differentiation with respect to $`\tau =t/t_0`$ and $`\kappa =K_pl^2/\epsilon _0`$ is the dimensionless proton-proton coupling constant.

As described in Introduction and illustrated by Fig. 1, the intrabond proton potential $`V(u)`$ can be formed as a result of superposition of two pair ion-proton interaction potentials placed tail-to-tail. According to the ab initio studies of Duan and Scheiner , these potentials are preferred to be chosen of the Morse type. As a result, the potential $`V(u)`$ is a symmetric double-well function with minima at $`u=\pm a`$ and a maximum at $`u=0`$ :

$$V(u)=\left[\frac{\alpha \mathrm{cosh}(\beta u)}{\alpha 1}\right]^2,\alpha =\mathrm{cosh}(\beta a).$$
(6)
The inequality $`\alpha >1`$ ensures the double-well form of the function (6). Throughout the paper we take $`a=0.25`$. The potential (6) is normalized so that the barrier height always equals unity. Its shape for different values of $`\beta `$ is shown in Fig. 2. As can be seen from this figure, the parameter $`\beta `$ determines the curvature/flatness of the barrier and its shape strongly depends on this parameter. We assume its values to range over the whole half-axis $`0<\beta <\mathrm{}`$. For small $`\beta `$, the barrier is rather narrow being the limiting case of the $`\varphi ^4`$ model, i.e.,

$$\underset{\beta 0}{lim}V(u)=\left(1u^2/a^2\right)^2.$$
(7)
Increase of $`\beta `$ makes the barrier more flat and the wells more narrow, so that the other limit is

$$\underset{\beta \mathrm{}}{lim}V(u)=\{\begin{array}{cccc}\mathrm{},& & \mathrm{}<u<a,& \\ 0,& & u=\pm a,& \\ 1,& & a<u<a,& \\ \mathrm{},& & a<u<\mathrm{}.& \end{array}$$
(8)
The corresponding equation of motion is the well known discrete nonlinear Klein-Gordon equation:

$$\ddot{u}_n=\kappa (u_{n+1}2u_n+u_{n1})V^{}(u_n),n=0,\pm 1,\mathrm{}.$$
(9)
Here and in what follows the prime denotes differentiation of a function with respect to its argument.

Before embarking on a complete analysis of the discrete equation (9), we calculate the dispersion law of small-amplitude waves around one of the two ground states. As a result, this law is given by the following equation:

$`\omega ^2(q)`$ $`=`$ $`\omega _0^2+2\kappa (1\mathrm{cos}q),`$ (10)
$`\omega _0`$ $`=`$ $`\sqrt{2{\displaystyle \frac{\alpha +1}{\alpha 1}}}\beta ={\displaystyle \frac{\sqrt{2}\beta }{\mathrm{tanh}(\beta a/2)}}.`$ (11)
The gap of the spectrum depends on the parameter $`\beta `$; for large values of $`\beta `$, it increases linearly with $`\beta `$.

## III Kink states, their stability and bifurcations

In this section, we start from the classification of possible stationary (anti)kink states. To compute these states, we have used both the conjugate gradients method for minimization of the stationary part ($`\dot{u}_n0`$) in the Hamiltonian (5) and the Newton iteration method for solving the time independent nonlinear set of equations that originates from Eqs. (9):

$$\kappa (u_{n+1}2u_n+u_{n1})=V^{}(u_n),n=0,\pm 1,\mathrm{}.$$
(12)
These equations can be rewritten as a two-dimensional map:

$`p_{n+1}`$ $`=`$ $`u_n,`$ (13)
$`u_{n+1}`$ $`=`$ $`\kappa ^1V^{}(u_n)+2u_np_n,n=0,1,\mathrm{}.`$ (14)
In general, maps of this type are chaotic. However, an on-site potential, for which the map (14) is integrable, has been found in Ref. . The previous knowledge on the kink solutions in the most popular cases of the discrete nonlinear Klein-Gordon lattices such as the sine-Gordon or $`\varphi ^4`$ chains implies the existence of only two stationary kink states. These states possess the inversion symmetry with respect to the center of the kink, being monotonic functions on the lattice. They connect two hyperbolic fixed points of the map $`(a,a)`$ and $`(a,a)`$, which are the ground states of the chain.

Adapted to our case of a HB chain with the numbering of ions and protons according to the sequence $`\{\mathrm{},Q_{n1},q_{n1},Q_n,q_n,Q_{n+1},q_{n+1},\mathrm{}\}`$, one of these stationary states, which has its center positioned at a heavy ion (call it an ion-centered kink/antikink), say with a number $`n_0`$, is dynamically stable, whereas the other, with its center positioned in the middle of a $`n_0`$th H-bond, i.e., in between the $`n_0`$th and the $`(n_0+1)`$th ions (call it a bond- or proton-centered kink/antikink), is dynamically unstable. The symmetry of the ion-centered (on the $`n_0`$th ion, in between H-bonds $`n_01`$ and $`n_0`$) kink/antikink is defined by the relations

$$u_{n_0n}=u_{n_0+n1},n=0,\pm 1,\mathrm{},$$
(15)
whereas for the proton-centered (in the middle of the $`n_0`$th H-bond, in between ions $`n_0`$ and $`n_0+1`$) kink/antikink, the symmetry relation is given by

$$u_{n_0n}=u_{n_0+n},n=0,\pm 1,\mathrm{}.$$
(16)
The solutions of these types certainly exist also in our model, as illustrated by Fig. 3. However, their stability properties appear to be much more complicated and they are discussed below.

In general, a deformation of the barrier shape in the potential $`V(u)`$ leads to a more rich family of kink/antikink solutions. This has been observed in the previous studies . The first feature caused by the deformation of the proton potential (6) with increase of $`\beta `$ is the phenomenon of stability switching, according to which the two types of kink solutions with inversion symmetry defined by Eqs. (15) and (16) switch their stability, while varying the system parameters. The second one is the appearance of new types of kink solutions. For the first of these types the symmetry relations (15) and (16) are not valid anymore, whereas the other one is symmetric with a zigzag-like profile at their center, but still has monotonic asymptotics as $`|n|\mathrm{}`$.

### A An exactly solvable limit

To understand better the effect of stability switching, it is instructive to consider the limit $`\beta \mathrm{}`$, resulting in a similar potential behavior as studied in Ref. . In this limit, the double-Morse potential takes the form \[see Eq. (8)\], for which the system of equations of motion (12) becomes exactly solvable because the particles (protons) can appear only either in the wells or on the flat region of the proton potential $`V(u)`$. There exists an infinite, but a countable set of the stationary kink solutions with an arbitrary number of particles, $`m=0,1,\mathrm{}`$, lying on the barrier. This number uniquely defines a kink/antikink solution. The set of these solutions can be written as

$$u_n=\{\begin{array}{ccc}a,& & \mathrm{}<nn_0\frac{m}{2}1,\\ 2a\frac{nn_0+1/2}{m+1},& & n_0\frac{m}{2}1<n<n_0+\frac{m}{2},\\ a,& & n_0+\frac{m}{2}n<\mathrm{},\end{array}$$
(17)
for the kink centered on the $`n_0`$th ion ($`m=0,2,\mathrm{}`$) and

$$u_n=\{\begin{array}{ccc}a,& & \mathrm{}<nn_0\frac{m+1}{2},\\ 2a\frac{nn_0}{m+1},& & n_0\frac{m+1}{2}<n<n_0+\frac{m+1}{2},\\ a,& & n_0+\frac{m+1}{2}n<\mathrm{},\end{array}$$
(18)
for the kink centered on the $`n_0`$th H-bond ($`m=1,3,\mathrm{}`$). Similarly, the analytical expressions for the antikink solutions can be obtained. These stationary solutions can be found, using either energy arguments based on the Hamiltonian (5) or directly from the equations of motion (12), where

$$\underset{\beta \mathrm{}}{lim}V^{}(u)=\{\begin{array}{ccc}\mathrm{},& & \mathrm{}<u<a,\\ 0,& & aua,\\ \mathrm{},& & a<u<\mathrm{}.\end{array}$$
(19)
The energy of both these (anti)kink solutions in the limit $`\beta \mathrm{}`$ is easily calculated and for any integer $`m=0,1,\mathrm{},`$ it reads

$$E_m=E_m(\kappa )=m+2\kappa a^2/(m+1).$$
(20)
As illustrated by Fig. 4, where the linear dependences of the energy $`E_m`$ on $`\kappa `$ are plotted for different $`m`$’s that the crossings of these dependencies occur at some values of the coupling parameter $`\kappa `$.

The most interesting points in Fig. 4 are the crossings for the states with the energies that correspond to adjacent $`m`$’s, i.e., $`m=0`$ and 1, $`m=1`$ and 2, and so on because they occur at lowest energies. Thus, the crossings when the kink state with $`m`$ particles on the barrier is transformed into the state with $`m+1`$ particles on it, occur at the following values of $`\kappa `$:

$$\kappa _{m+1}^{(c)}=(m+1)(m+2)/2a^2,m=0,1,\mathrm{}.$$
(21)
Therefore, depending on the strength of the proton-proton coupling $`\kappa `$, the proton-centered or the ion-centered kink can reach a global minimum of the energy $`E_m>0`$ (except for the ground states, when $`E=0`$). In addition, it is interesting to notice that at the values of $`\kappa `$ given by Eqs. (21), the inter-bond and the intrabond energies are equal exactly each other.

### B Numerical results for finite $`\beta `$’s

Now let us investigate how the properties of the kink solutions found in the exactly solvable limit $`\beta \mathrm{}`$ change, when $`\beta `$ take finite values. Thus, changing $`\beta `$ allows us to explore the whole set of scenarios, starting from the $`\varphi ^4`$ limit and finishing with the limit $`\beta \mathrm{}`$. We start to compute the kink solutions from the anticontinuous limit ($`\kappa =0)`$, taking the solutions of the exactly solvable limit as an initial guess for the Newton iteration method. Then, we increase $`\kappa `$ and check how the kink profiles behave. In Fig. 5, the energy dependence on the coupling parameter $`\kappa `$ is plotted for $`\beta =10`$. First, let us focus on the behavior of the solutions with the lowest energies (namely those which correspond to $`m=0`$ and $`m=1`$ in the exactly solvable limit).

In Fig. 5, curve 1 corresponds to the ion-centered kink with $`m=0`$, the symmetry of which is given by Eqs. (15), whereas curve 2 to the proton-centered kink with $`m=1`$, the symmetry of which is given by Eqs. (16). Contrary to the previous knowledge on the stability properties of the kink solutions for the discrete sine-Gordon and $`\varphi ^4`$ models, where the (anti)kinks of the symmetry (15) are always dynamically stable, while the (anti)kinks of the symmetry (16) are always dynamically unstable, an interchange of stability is observed for the proton potential (6) with finite $`\beta `$’s as $`\kappa `$ varies. Thus, one can see that at a certain value of the coupling parameter, $`\kappa 24.5`$, the energies of the both types of symmetric kinks coincide and after passing this critical point, the proton-centered kink appears to be stable, while the ion-centered kink is unstable. For higher $`\kappa `$’s, several more interchanges of stability take place, with the energy difference between the sequential kink states that decreases with the growth of $`\kappa `$. These transitions of stability take place smoothly all the way up to the continuum limit. We refer to these transitions as to stability switchings. Thus, the solutions with $`m>1`$, which were clearly separated from the $`m=0`$ and the $`m=1`$ solutions in the exactly solvable limit $`\beta \mathrm{}`$, appear to be smoothly connected with them. In other words, while the coupling $`\kappa `$ increases, the particles slowly “climb” on the barrier, so that there is no abrupt transition from the state with $`m=0`$ to the state with $`m=2`$ or, further, to the states with $`m=4,6,\mathrm{}`$. The same can be concluded about the kinks with $`m`$ being odd. In this respect, the system still shows similarity with the $`\varphi ^4`$ model.

If we focus more carefully on the behavior of the system in the vicinity of the points, where the energies of proton-centered and ion-centered kinks become equal, we find that a new type of kinks appears. These kinks shown in Fig. 6 do not exhibit any of the symmetries defined by Eqs. (15) and (16), but they are doubly degenerate related to each other by the inversion with respect to the crossing point between the line $`u_n=0`$ and the line connecting two central particles of the kink \[$`n=48`$ and $`n=49`$ in Fig. 6(a)\]. The energy of these kink states with broken symmetry is always larger than the energy of the symmetric (ion-centered and proton-centered) kinks. The asymmetric (anti)kinks are always linearly unstable. In Fig. 5, their energy is shown by curve 3.

The kinks of the last type are shown in Fig. 5 by curves 4 and 5. This is what happens to the solutions with $`m>1`$ obtained in the limit $`\beta \mathrm{}`$, when we are moving from the anticontinuous limit ($`\kappa =0`$). Instead of attaining a regular monotonic form, the kink profiles with several particles on the barrier develop a zigzag-like structure in their centers, as demonstrated by Fig. 7. For finite $`\beta `$’s the zigzag-like solutions are linearly unstable and therefore we do not study them here in more detail.

Instead, we focus on the behavior of the monotonic kink solutions with increase of the coupling $`\kappa `$ in the vicinity of the stability switchings. To understand this effect better, in Fig. 8, we have plotted the position of the $`(N/2)`$th proton ($`N`$ is the total number of H-bonds in the chain) as a function of the coupling $`\kappa `$. Here the sequence of pitchfork bifurcations is clearly seen. Curve 1 corresponds to the kink centered in the middle of the $`(N/2)`$th bond, i.e., when the $`(N/2)`$th proton is always a central particle of the kink with $`u_{N/2}0`$. The displacements of the $`(N/2)`$th proton for the ion-centered kinks positioned on the $`(N/2+1)`$th and the $`(N/2)`$th ions are shown by curves 2 and 3, respectively. When increasing $`\kappa `$, the $`(N/2)`$th proton moves slowly out of the well. At a certain value of $`\kappa `$, more specifically, at $`\kappa 22.5`$, the pitchfork bifurcation of the proton-centered kink takes place. This configuration retains its stability and two new solutions (both linearly unstable) appear. These are precisely those asymmetric kinks (see curves 4 and 5 in Fig. 8), with their shape shown in Fig. 6. At the beginning, they look like slightly distorted proton-centered kinks, but with the growth of $`\kappa `$, they change more and more towards the ion-centered configuration. Eventually, the second pitchfork bifurcation takes place at $`\kappa 26.2`$. The asymmetric kinks join the ion-centered kinks (junction of curves 3 and 4, and curves 2 and 5) and the ion-centered configuration loses its stability. Now one can clearly see that two identical kinks, shifted by one lattice spacing with respect to each other, are connected via such a bifurcation sequence. Thus, this pitchfork bifurcation is nothing but a transition of the kink from one position to another position, one lattice period forward or backwards. This cascade of bifurcations can be continued further up or down in $`u_n`$’s or, in other words, two or more sites backwards or forward. Similar bifurcation scenario for discrete breathers in the ac-driven and damped Klein-Gordon lattice has been reported in .

For higher values of $`\beta `$, the effect of stability switchings exists, being more pronounced because the switchings start at smaller $`\kappa `$ and take place more frequently (see Fig. 9). Another feature that appears from the exactly solvable limit $`\beta \mathrm{}`$ is as follows. Curve 1 of Fig. 9 corresponds to the ion-centered kink and curve 2 to the proton-centered kink. They cross each other at $`\kappa 17.4`$, where the ion-centered kink loses its stability, disappearing shortly at $`\kappa 28.6`$. Here the same pitchfork bifurcation scenario takes place, but the distance between the first and the second bifurcations is much larger than in the $`\beta =10`$ case. In the meanwhile, a bit earlier, at $`\kappa 23.4`$, a new family of ion-centered kinks appears (curve 4). This curve corresponds to the ion-centered kinks with two protons on the barrier. Thus, one can observe the coexistence of two different kink solutions with the same type of symmetry (for more details see the upper inset of Fig. 9). This coexistence takes place on a rather narrow interval of $`\kappa `$ and one of the coexisting kinks is unstable, but still this phenomenon clearly originates from the limit $`\beta \mathrm{}`$. In Fig. 10, we show the shape of two coexisting kinks. One of them (more narrow) corresponds to curve 1 of Fig. 9 and the second one, which is more broad, corresponds to curve 4 of this figure. These kinks have zero and two protons on the barrier, respectively. When $`\kappa `$ increases further, the stability switchings occur between curves 2 and 4 under the same scenario as before for $`\beta =10`$ (see for details the lower inset of Fig. 9). Zigzag-like kinks are also presented in this case, as shown by curves from 6 to 8.

When increasing $`\beta `$ even more, the coexistence of different kinks with the same symmetry is even more pronounced. We have checked the case of $`\beta =50`$ and discovered for this value that several cases of this coexistence for different kinks take place and they are much more pronounced. Thus, this case is rather close to the exactly solvable limit $`\beta \mathrm{}`$.

### C Analytical approximation for finding kink solutions

One can use a simplification of the potential (6) in order to obtain analytically an exact (anti)kink solution. To this end, we approximate both the wells of the potential by parabolas connected by the constant equal to the barrier height as follows:

$$V(u)\{\begin{array}{ccc}(\omega _0^2/2)(u+a)^2,& & \mathrm{}<u<b,\\ 1,& & bub,\\ (\omega _0^2/2)(ua)^2,& & b<u<\mathrm{},\end{array}$$
(22)
where $`b=a\sqrt{2}/\omega _0`$. Schematic description of this “parabola-constant” approximation is presented in Fig. 11 by thick solid lines. The thin solid lines show the original potential (6) with $`\beta =20`$. The approximation is expected to work well, when the barrier is flat enough, i.e., when $`\omega _0\sqrt{2}/a`$. Within this approximation, we are able to solve the problem of finding stationary kink solutions analytically.

Let $`m=0,1,\mathrm{},`$ be the number of protons on the barrier of the potential (22). Then the discrete kink profiles are given by

$$u_n=\{\begin{array}{cccccc}a+A_me^{\lambda (nn_0+m/2+1)}& & & & & \\ \text{if}\mathrm{}<nn_0m/21,& & & & & \\ (nn_0+1/2)D_m& & & & & \\ \text{if}n_0m/21<n<n_0+m/2,& & & & & \\ aA_me^{\lambda (nn_0m/2)}& & & & & \\ \text{if}n_0+m/2n<\mathrm{},& & & & & \end{array}$$
(23)
for the kink centered on the $`n_0`$th ion ($`m=0,2,\mathrm{}`$) and

$$u_n=\{\begin{array}{cccccc}a+A_me^{\lambda [nn_0+(m+1)/2]}& & & & & \\ \text{if}\mathrm{}<nn_0(m+1)/2,& & & & & \\ (nn_0)D_m& & & & & \\ \text{if}n_0(m+1)/2<n<n_0+(m+1)/2,& & & & & \\ aA_me^{\lambda [nn_0(m+1)/2]}& & & & & \\ \text{if}n_0+(m+1)/2n<\mathrm{},& & & & & \end{array}$$
(24)
for the kink solution centered on the $`n_0`$th H-bond ($`m=1,3,\mathrm{}`$). Here $`\lambda `$ is a “localization” parameter that measures the transition width in the (anti)kink profile between the uniform distribution of protons on the barrier and the (anti)kink asymptotics $`u_n\pm a`$. It is given by a positive root of the equation

$$\mathrm{cosh}\lambda =1+\omega _0^2/2\kappa .$$
(25)
The other two parameters, the amplitude $`A`$ and the uniform distance between the nearest-neighbor protons on the flat $`D`$ can be expressed through the localization parameter $`\lambda `$ as

$`A_m`$ $`=`$ $`{\displaystyle \frac{2a}{(m+1)e^\lambda m+1}},`$ (26)
$`D_m`$ $`=`$ $`{\displaystyle \frac{2ae^\lambda }{(m+1)e^\lambda m+1}}.`$ (27)
As follows from Eqs. (11) and (25), the limit $`\lambda \mathrm{}`$ (when $`\omega _0^2/\kappa \mathrm{}`$) is more general than the limit $`\beta \mathrm{}`$ because $`\lambda `$ contains both $`\beta `$ and $`\kappa `$. Therefore one can check that Eqs. (23)-(27) are reduced to the stationary kink solution given by Eqs. (17) and (18) as $`\lambda \mathrm{}`$. In particular, the amplitude $`A_m`$ and distance between the protons on the barrier $`D_m`$ tend to zero and $`2a/(m+1)`$, respectively.

Using Eqs. (23)-(27) one can easily compute the energy of both the kink configurations:

$$E_m=m+2\kappa a^2\frac{\mathrm{tanh}(\lambda /2)}{1+m\mathrm{tanh}(\lambda /2)}.$$
(28)
Similarly, this expression is also transformed to the energy (20) as $`\lambda \mathrm{}`$.

Now we investigate the behavior of the energy difference

$$\mathrm{\Delta }E_{m+1}(\kappa ,\beta )=E_{m+1}E_m,m=0,1,\mathrm{}.$$
(29)
We find that this difference as a function of $`\kappa `$ has a number of zeroes and these zeroes depend on the parameter $`\beta `$. Thus, within our approximation one can predict the effect of switching of the stable and unstable kink configurations. In Table I, we show the values of the coupling parameter $`\kappa `$, for which the first \[i.e., when $`\kappa =\kappa _1`$ and $`m=0`$; see also Eq. (21)\] switching of the kink stability takes place.

We see that the approximation works fairly well when $`\beta `$ is rather large and the barrier between the wells is close to being completely flat. It is improving with increase of $`\beta `$ and in the limit $`\beta \mathrm{}`$, our approximation coincides with the exact result, shown above. Thus, switchings of the stability of kink states with different symmetries is a generic effect that does not depend on a specific model, but on the properties of the on-site potential. Another piecewise approximation of the proton potential $`V(u)`$ for a HB chain has been constructed earlier by Weiner and Askar , using alternatively inversed parabolas. Our potential approximation (22) seems to be more appropriate for the studies in detail of flatness effects of the on-site (intrabond) potential and more close to the realistic double-Morse potential (6) if $`\beta `$ is not so small.

Note that the zigzag-like kink profiles obtained above numerically shown, e.g., in Fig. 7 can also be given analytically within the approximation (22). Indeed, the ion-centered kink shown in Fig. 7(a) is described by

$$u_n=\{\begin{array}{ccc}a+B_0e^{\lambda (n+2)},n=2,3,\mathrm{},& & \\ u_0=u_1=\xi _0,& & \\ aB_0e^{\lambda (n1)},n=1,2,\mathrm{},& & \end{array}$$
(30)
whereas the proton-centered profile illustrated by Fig. 7(b) is given by

$$u_n=\{\begin{array}{ccc}a+B_1e^{\lambda (n+2)},n=2,3,\mathrm{},& & \\ u_0=0,u_1=u_1=\xi _1,& & \\ aB_1e^{\lambda (n1)},n=1,2,\mathrm{},& & \end{array}$$
(31)
where $`\lambda `$ is given by Eq. (25) and

$`B_m`$ $`=`$ $`2a{\displaystyle \frac{e^\lambda +e^\lambda m/21}{e^\lambda \left(e^\lambda m+1\right)}},`$ (32)
$`\xi _m`$ $`=`$ $`a{\displaystyle \frac{3e^\lambda 2e^\lambda }{e^\lambda m+1}},m=0,1.`$ (33)
Using the last equations, the kink energy can be calculated in a similar way as the energy (20). In the limit $`\lambda \mathrm{}`$, the energy of the zigzag-like kinks, profiles of which are illustrated by Fig. 7, become

$$E_0=6\kappa a^2\text{and}E_1=1+5\kappa a^2.$$
(34)
These expressions clearly show that for this type of kinks their energies also can coincide for certain values of coupling $`\kappa `$ if $`\beta `$ (and, consequently, $`\lambda `$) is large enough.

### D Elementary excitations on the kink background

Let $`u_n^{(0)}`$ be a stationary kink solution of Eq. (12), i.e., a fixed point of the map (14). We are interested in the properties of small-amplitude excitations on the kink background. Linearizing Eq. (9) around the stationary kink solution according to

$$u_n=u_n^{(0)}+A_ne^{i\mathrm{\Omega }\tau },$$
(35)
we arrive at the eigenvalue problem

$$\widehat{L}𝐀=\mathrm{\Lambda }𝐀,𝐀=\{\mathrm{},A_{n1},A_n,A_{n+1},\mathrm{},\}.$$
(36)
Here the operator $`\widehat{L}`$ \[the Hessian of the Hamiltonian Eq. (5)\] acts on a vector $`𝐀`$ as

$$(\widehat{L}𝐀)_n=\kappa (A_{n+1}2A_n+A_{n1})+V_nA_n,$$
(37)
where $`V_n=V^{\prime \prime }[u_n^{(0)}]`$. The operator $`\widehat{L}`$ is a symmetric (so all eigenvalues are real) tridiagonal matrix and the spectral parameter is $`\mathrm{\Lambda }\mathrm{\Omega }^2`$. This eigenvalue problem can be treated as a quantum-mechanical problem of a particle, trapped in a single-well spatially discrete potential formed by the kink. Its depth depends on the curvature/flatness of the proton potential in the middle of the H-bond, $`V^{\prime \prime }(0)`$, and tends to $`\omega _0^2`$ as $`n\pm \mathrm{}`$.

The eigenvalues of the problem also give information about the linear stability of the kink solution. If there exists at least one eigenvalue $`\mathrm{\Lambda }=\mathrm{\Omega }^2<0`$, the linear excitation on the kink grows exponentially in time and the corresponding kink solution is linearly unstable. Otherwise, it is linearly stable. The stability of stationary kink solutions is determined by the system parameters (in our case, by the curvature/flatness parameter $`\beta `$ and the coupling $`\kappa `$).

In Fig. 12, we depict the dependence of the eigenfrequencies $`\mathrm{\Omega }_n`$’s on the coupling parameter $`\kappa `$ for different curvatures $`\beta `$. The spectrum consists of two parts. The first one describes the eigenfrequencies that lie above the potential $`V_n`$ and therefore determine delocalized eigenvectors corresponding to the linear (phonon) spectrum of the lattice. The second part consists of eigenfrequencies that lie below $`\omega _0^2`$ corresponding to spatially localized eigenvectors. The lowest eigenfrequency is the Goldstone mode, which is universally presented in all the Klein-Gordon models.

The second mode corresponds to small-amplitude oscillations of the kink core, being often referred to as the Rice mode. In the continuum $`\varphi ^4`$ model, there are only two localized modes: the Goldstone mode and the Rice mode. The properties of the internal modes can be significantly altered due to different factors such as change of the nature of the interparticle interaction  or change of the shape of the on-site potential . If $`\beta `$ is not so large, the behavior of the eigenfrequencies and eigenvectors (see Fig. 13) of our system is reminiscent to that of the $`\varphi ^4`$ model, as shown in panels (a) and (b) of Fig. 12. However, several differences occur. Thus, the Goldstone mode collides with the zero axis and becomes unstable \[see Fig. 12(a)\] for the ion-centered kink. Meanwhile, for the proton-centered kink \[see Fig. 12(b)\], the Goldstone mode was initially unstable and became stable later on. The stability switchings are caused by the pitchfork bifurcations as described at the beginning of this section. Another difference is the appearance of the new eigenfrequency for the proton-centered kink with the eigenvector which has two nodes \[see Fig. 13(e)\].

Thus, one can see that curvatures play an important role in the properties of the eigenfrequencies of our system. Increasing $`\beta `$ further, we observe more pitchfork bifurcations of the Goldstone mode and the appearance of more localized modes. In panels (c) and (d), we show how the eigenfrequencies behave for $`\beta =10`$. We observe more localized internal modes, some of them surviving in the continuum limit and some of them disappearing there. The shape of the corresponding eigenvectors (see Fig. 14) behaves accordingly to the wavefunction shape of the bound states of the quantum-mechanical Schrödinger equation.

## IV Properties of the Peierls-Nabarro barrier

In general, the propagation of topological solitons (kinks and antikinks) in lattices are subject to their discreteness. The discreteness effects can be described by a spatially periodic potential with the period coinciding with the lattice spacing $`l`$, known as a Peierls-Nabarro (PN) potential. The PN potential $`E_{PN}(n_c)`$ is a function of the kink center position

$$n_c=\underset{n}{}n\frac{u_{n+1}u_{n1}}{2(u_{\mathrm{}}u_{\mathrm{}})}=\frac{1}{4a}\underset{n}{}n(u_{n+1}u_{n1}),$$
(38)
and the height of the PN barrier equals

$$\delta E=E(n_{c,max})E(n_{c,min}).$$
(39)
This difference measures the activation energy for the kink propagation though one lattice spacing. Here $`n_{c,max}`$ and $`n_{c,min}`$ are the positions of the kink, where the PN potential has its maxima or minima, respectively. The minimum position $`n_{c,min}`$ coincides with the kink position, which is a global minimum of the kink energy. Thus, normally $`n_{c,min}`$ is an integer or a half-integer.

It is a well established fact that for the $`\varphi ^4`$ and sine-Gordon models, the height of the PN barrier coincides with the difference between the energies of the proton-centered and ion-centered kinks. As demonstrated above, the deformation of the proton (intrabond) potential $`V(u)`$ leads to the switching of stability of these two states. This does not mean, however, that the PN barrier disappears, when the energy of these states coincides because the site-centered and proton-centered states do not represent the states with the highest and the lowest energy of the PN potential .

In Fig. 15, some examples of variation of the energy $`E(n_c)`$ for different values of $`\kappa `$ are presented. Panel (a) corresponds to the situation, when the proton-centered kink is an energy minimum and the ion-centered one reaches a maximum. The second panel (b) demonstrates the case, when the first pitchfork bifurcation occurs and the asymmetric kink appears. Panel (c) corresponds to the case, when the stability switching happens and both the proton- and ion-centered kink states get minima, and the asymmetric kinks are maxima. The period of the PN potential is decreased by one half. After that, as shown in panel (d), the proton-centered kink becomes unstable and the ion-centered one stable. The last panel (e) is obtained after the second pitchfork bifurcation has occured and the asymmetric kinks disappear.

Dependence of the PN barrier is non-monotonic not only in $`\kappa `$ but also on $`\beta `$. In Fig. 16, we show the dependence of the PN barrier on the parameter $`\beta `$ for a fixed value of $`\kappa `$. We observe a minimum at $`\beta 8.3`$, which is obviously caused by the stability switching. Indeed, as shown by the upper inset, the minimum of the PN potential moves from a half-integer value to an integer one. Thus, the stable kink configuration changes from the ion-centered to the proton-centered kink. Note, that the height of the PN barrier does not attain zero (see the lower inset), but decreases by two order of magnitude.

## V Discrete kink mobility

Another interesting consequence of the deformation of the on-site potential barrier is the possibility of the existence of non-oscillating travelling kinks. In the continuum Klein-Gordon models which admit moving topological solitons (kinks and antikinks), the domain of admissible soliton velocities $`s`$ is the interval $`0s<c_0`$, where $`c_0=\sqrt{\kappa }`$ is the characteristic velocity. Thus, kinks in these continuum models are a one-parametric family of solutions with the kink velocity $`s`$ as a parameter. In a general case of the discrete Klein-Gordon model, this family is reduced to a discrete set of travelling kink solutions with some velocities $`s_0`$, $`s_1,\mathrm{},`$ $`s_k`$, $`s_0=0`$, $`s_k<c_0`$. In the $`\varphi ^4`$ and the sine-Gordon models, there exists only $`s_0`$. In general, everywhere in between $`s_n`$’s, there exist moving kinks with oscillating asymptotics known as nanopterons. The existence of velocities $`s_n0`$ has been shown both numerically  and analytically  for several models. This effect is due to the properties of the on-site (intrabond) potential $`V(u)`$ and, more precisely, the flatness of its barrier which, in its turn, causes the stability switchings described in the previous sections.

For finding non-oscillating kink solutions, we have used a pseudospectral method . The method allows us to find the travelling-wave kink solutions of the type

$$u_n(t)=u(ns\tau )u(z),$$
(40)
solving the differential equation with advanced and delay terms:

$$s^2u^{\prime \prime }(z)=\kappa [u(z+1)2u(z)+u(z1)]V^{}[u(z)].$$
(41)
Dependence of the kink velocity on the coupling parameter $`\kappa `$ is shown in Fig. 17. As follows from this figure, it starts (see curve 1) at the value of $`\kappa `$, which is close to that when the pitchfork bifurcation takes place (see Fig. 8). Then the velocity grows with $`\kappa `$. The second velocity dependence starts at $`\kappa `$ being close to the point of the second pitchfork bifurcation.

Dependence of the kink velocity on the parameter $`\beta `$ has similar behavior, as demonstrated by Fig. 18. It is interesting to note that the value of the velocity of the moving kink does not decrease down to zero as $`\kappa `$ or $`\beta `$ decreases, but stops at some finite $`s`$ instead. This means that the kinks should have sufficient kinetic energy to overcome the pinning effects.

In Fig. 19, we plot the profiles $`u(z)`$ of the moving kinks that correspond to both two curves of Fig. 17. The profile in panel (a) corresponds to the moving kink from curve 1, whereas the profile in panel (b) to the kink from curve 2. In both the cases, the coupling constant was the same: $`\kappa =100`$. The second kink appears to be wider and this can be explained by the fact that its velocity is more than two times smaller. The insets of the figures show the velocity profiles.

In Fig. 22, we investigate the behavior of moving kinks for different $`\beta `$’s. Increase of the barrier flatness is reflected in change of the shape at the kink center. We have considered the kinks corresponding to curve 2 of Fig. 18, when the coupling is fixed: $`\kappa =120`$. Deformation of the kink profile can easily be seen on the velocity profile shown in the inset. The kink profile experiences deformation of its slope part, which is seen as a dip in the velocity profile. The dip grows with increase of $`\beta `$. More flat the barrier becomes, more possibilities of the kink’s profile deformation.

These results demonstrate that the existence of a finite set of velocities, at which the kink can move with constant shape is a generic effect. Its appearance is due to the shape of the on-site (intrabond proton) potential. If the proton potential has a barrier which is flat enough to allow the symmetry switchings (accompanied by the pitchfork bifurcations), the Peierls-Nabarro barrier experiences lowering for specific values of the system parameters. Thus, it becomes possible that the kinetic energy of the kink is sufficient to overcome the pinning forces of the lattice.

## VI Conclusions

We have studied the dynamics of the one-dimensional Klein-Gordon lattice with the on-site potential of the double-Morse type. This is a physically motivated model, which is the simplest one for the proton transport in a hydrogen-bonded chain, where the on-site potential plays the role of the potential for proton transfers in the hydrogen bond. Therefore throughout this paper we call it an it intrabond proton potential. A Morse-type function was found to offer the best combination of accuracy in reproducing quantum-mechanically computed potentials . The model has two parameters, the proton-proton coupling $`\kappa `$ and the anharmonicity of the Morse potential, the curvature parameter $`\beta `$. The anharmonicity parameter is responsible for the shape of the intrabond potential $`V(u)`$, especially on the convexity of its barrier. For larger $`\beta `$, the barrier becomes more flat and the wells become more narrow. Changing this parameter, one can explore the variety of possible intrabond potentials, starting from the $`\varphi ^4`$ model (as $`\beta 0`$) and finishing with the exactly solvable limit $`\beta \mathrm{}`$.

Deformation of the barrier of the intrabond potential becomes crucial for the properties of the stationary kink solutions. While the $`\varphi ^4`$ limit allows us only two types of stationary kinks: ion-centered and proton-centered with their stability properties being constant for any coupling $`\kappa `$, the opposite limit shows the existence of an infinite countable set of stationary kink solutions. For some values of the coupling $`\kappa `$ forming an infinite countable set, the states with different symmetries can have the same energy. In between these two limiting cases, some of the kink properties survive from the exactly solvable limit $`\beta \mathrm{}`$. First of all, the stability properties of the kink solutions depend drastically on the coupling parameter $`\kappa `$. With increase of $`\kappa `$, the initially stable ion-centered kink becomes unstable, while the proton-centered kink gets stable. With further increase of $`\kappa `$, these stability switchings which are, in fact, pitchfork bifurcations, go on several times. Another result can be seen for rather high $`\beta `$, at least, for $`\beta >20`$, in our calculations. It is the coexistence of several kink solutions of the same symmetry for the same $`\kappa `$. This is a left-over from the exactly solvable limit and it disappears with lowering $`\beta `$. An analytical approximation has been constructed to show the effect of the symmetry switchings analytically and to confirm that the effect is not confined to a specific model, but has a universal nature.

The stability switchings contribute to the non-monotonic behavior of the Peierls-Nabarro (PN) barrier. The barrier decays with increase of $`\kappa `$, however, experiences local minima at the stability switching points, where its value decreases by one order of magnitude. The same happens for the dependence of the PN barrier on $`\beta `$. This phenomenon, in its turn, assists the kink mobility and leads to the appearance of a finite set of velocities for which the propagation of very narrow kinks is possible. This is a generic effect, attributed to the shape of the on-site (intrabond) potential, and also is not confined to a specific model. Note that the PN barrier is not required to vanish completely (see Ref. ). Simply, the barrier is low enough for the kinetic energy of the kink to carry it over the barrier.

Thus, we have concluded that highly mobile kinks are possible in our model of the proton transport in hydrogen-bonded chains even for those proton-proton nearest-neighbor interactions, when the proton kinks are very narrow. This gives us a reason to believe that the soliton mechanism of proton transfers can work for physically reasonable values of the proton-proton coupling $`\kappa `$.

This work has benefited from discussions with M. Peyrard. We acknowledge financial support from the RTN Project No. LOCNET HPRN-CT-1999-00163, the INTAS Grant No. 97-0368, and the Danish Research Agency.


format: text

text:
# Modeling generalized parton distributions to describe deeply virtual Compton scattering data

## I Introduction

Generalized parton distributions (GPDs) have been studied extensively in recent years mrgdh ; ji ; rad1 ; die ; mfort ; vgg ; jcaf ; ffgs ; ffs ; gb1 ; gb2 ; bemuothers . This interest was spurred by the realization that these distributions are not only the basic, non-perturbative ingredients in hard, exclusive processes such as deeply virtual Compton scattering (DVCS), or exclusive vector meson production, but that they are generalizations of the well known parton distribution functions (PDFs) from inclusive reactions. GPDs incorporate both a partonic and distributional amplitude behavior and hence contain more information about the hadronic degrees of freedom than PDFs. In fact, GPDs are true two-parton correlation functions, allowing access to the highly non-trivial parton correlations inside hadrons foot .

GPDs can be broadly characterized by the following features:

* They depend on two momentum fraction variables, a partonic variable defined with respect either to the incoming or the average of the incoming and outgoing proton momentum and the skewedness (which is the difference between the momentum fractions of two adjacent partons in the parton ladder).
* For fixed skewedness, they are continuous functions of the dependent variable and span two distinct regions, the DGLAP region and the ERBL region, in which their evolution in scale obeys generalized versions of the DGLAP and ERBL evolution equations, respectively, and in which their behavior is qualitatively different.
* They are even functions of the skewedness variable and the singlet, non-singlet and gluon distributions are either symmetric or anti-symmetric about the center point of the ERBL region (the symmetry obeyed depends on the precise definitions used).
* The Lorentz structure of their definitions implies a polynomiality condition ji ; rad1 ; gb2 ; poly : their ($`N1`$)-th moments are polynomials in square of the skewedness of degree no greater than $`N/2`$.
* They reduce to the ordinary PDFs in the limit of zero skewedness (the ‘forward limit’).

All of the above features have to be preserved under evolution in scale.

Any suggested model of GPDs should adhere to these mathematical features. In rad2  such a model, based on double distributions (DDs), was suggested for the GPD input distributions (see also rad3 ). In poly  it was realized that an additional term, the so-called D-term, was required in the ERBL region for the unpolarized quark singlet and gluon distributions in order to satisfy polynomiality for even $`N`$. The use of factorized factor  DDs augmented with a D-term has become a popular phenomenological model. Unfortunately, when this type of model for input GPDs was used in its current form to calculate deeply virtual Compton scattering at both leading (LO) and next-to-leading order (NLO), the results were not in agreement with the H1 data h1  on the DVCS photon level cross section, $`\sigma (\gamma ^{}p\gamma p)`$, and the HERMES and CLAS data hermclas  on the DVCS single spin asymmetry or charge asymmetry afmmshort ; afmmlong ; bemuk1 .

Another popular model for input GPDs, inspired by the aligned jet model (AJM) ajm  and its QCD extension fs88  is based on the observation that at a scale $`Q^212\text{GeV}^2`$ and a wide range of $`x_{bj}`$, soft physics gives the dominant contribution to the parton densities. As a result the effect of skewedness at small $`x_{bj}`$ should be rather small and hence at the input scale, it is a good approximation to set the GPDs equal to the forward PDFs at the same parton fraction, $`X`$, defined with respect to the incoming proton ffgs  (for any skewedness). This has the advantage that it automatically satisfies the requirements of polynomiality for the first two moments, however one encounters infinities in the quark singlet GPD in the middle of the ERBL region.

Another ‘forward model’, which may be considered to be an extreme case of a DD model, was adopted in gb2  where one assumes that the GPD is equal to the forward PDF at the same parton momentum fraction, $`v`$, with respect to the average of the incoming and outgoing proton momentum $`\overline{P}=(p+p^{})/2`$ which implicitly contains the skewedness. This translates to an $`X`$ which is shifted to lower values by an amount controlled by the skewedness. This ansatz works fine for the DGLAP region. Unfortunately in the ERBL region it also involves sampling the forward PDFs right the way down to zero in momentum fraction where they have not yet been measured (this is especially problematic for singular quark distributions). In this paper we construct an alternative, finite ‘forward model’ for the input GPDs, using the forward input PDFs in the DGLAP region and imposing a simple form in the ERBL region that has the correct symmetries and ensures polynomiality is respected in the first two moments (see gb2  for alternative ways of dealing with this problem). As we will demonstrate this forward model reproduces the available data on DVCS reasonably well.

Throughout this paper we will use the off-diagonal representation of GPDs, $`^i(X,\zeta )`$, defined by Golec-Biernat and Martin gb1  and used, for example, in the numerical solution of the GPD evolution equations in afmmgpd  (see gb1 ; bemuevolve ; radmus  for other approaches to numerical evolution). They depend on the momentum fraction $`X[0,1]`$ of the incoming proton’s momentum, $`p`$, and the skewedness variable $`\zeta =\mathrm{\Delta }^+/p^+`$ (so that $`\zeta =x_{bj}`$ for DVCS). For the quark case, the relationship of the quark and anti-quark distributions, $`^q(X,\zeta ),^{\overline{q}}(X,\zeta )`$, to Ji’s GPD $`H^q(v,\xi )`$ is shown in Fig. 1 with $`\xi =\mathrm{\Delta }^+/2\overline{P}^+`$. More explicitly, for $`v[\xi ,1]`$:

$$^{q,a}\left(X=\frac{v+\xi }{1+\xi },\zeta \right)=\frac{H^{q,a}(v,\xi )}{1\zeta /2},$$
(1)
and for $`v[1,\xi ]`$

$$^{\overline{q},a}\left(X=\frac{\xi v}{1+\xi },\zeta \right)=\frac{H^{q,a}(v,\xi )}{1\zeta /2}.$$
(2)
The two distinct transformations between $`v`$ and $`X`$ for the quark and anti-quark cases are shown explicitly on the left hand side of Eqs.(1, 2) . There are two distinct regions: the DGLAP region, $`X>\zeta `$ ($`|v|>\xi `$), and the ERBL region, $`X<\zeta `$ ($`|v|<\xi `$). In the ERBL region, due to the fermion symmetry, $`^q`$ and $`^{\overline{q}}`$ are not independent. In fact $`^q(X,\zeta )=^{\overline{q}}(\zeta X,\zeta )`$, which leads to an anti-symmetry of the unpolarized quark singlet distributions (summed over flavor $`a`$), $`^S=\mathrm{\Sigma }_a^{q,a}+^{\overline{q},a}`$, about the point $`\zeta /2`$ (the non-singlet and the gluon, $`^g`$, which is built from $`vH_{Ji}^g(v,\xi )`$, are symmetric about this point). $`^{NS,g}`$ are constructed analogously to $`^S`$ from $`H^{NS,g}`$.

The rest of the paper is organized as follows. Section II contains a detailed explanation of why DD-based models in their current form cannot describe the data. In Section III we construct our alternative forward model for input GPDs, which is motivated by the AJM ajm ; fs88  and describes the data well. In Section IV we propose a phenomenological model for the slope of the $`t`$-dependence in which the slope parameter is allowed to change with photon virtuality, $`q^2=Q^2`$. The model improves the theoretical description of the $`Q^2`$-dependence of the HERA data, relative to using a constant slope. Finally we summarize our findings in Section V.

## II The problem with double distribution models

In this section we discuss factorized double distribution (DD) based models for GPDs and explain in detail why the sampling of the forward PDF at extremely small $`x`$ in constructing the GPD leads to a problem in the quark singlet GPD.

Symmetric DDs, $`F_{DD}(x,y,t,Q^2)`$, were introduced in rad1  with plus momentum fractions, $`x,y`$, of the outgoing and returning partons defined as shown in the left hand plot of Fig. 2. They exist on the diamond-shaped domain shown to the right of Fig. 2. The outgoing parton lines of course only have a single plus momentum relative to any particular external momenta, so the GPDs are related to these DDs via a reduction integral, involving $`\delta (vx\xi y)`$, along the off-vertical lines in the diamond (the dotted line corresponds to $`v=\xi `$):

$`H(v,\xi )={\displaystyle _1^1}𝑑x^{}{\displaystyle _{1+|x^{}|}^{1|x^{}|}}𝑑y^{}\delta (x^{}+\xi y^{}v)F_{DD}(x^{},y^{}).`$ (3)
In rad2 ; rad3  a model for $`F_{DD}(x,y,t,\mu ^2)`$, at the input scale $`\mu ^2=Q_0^2`$, was introduced in which the functional form is factorized in $`x`$ and $`y`$ and in the $`t`$-dependence, an assumption borne out of convenience rather than physical necessity:

$$F_{DD}^i(x,y,\mu ^2,t)=\pi ^i(x,y)f^i(x,\mu ^2)r^i(t).$$
(4)
Here $`f^i`$ and $`r^i`$ are the standard PDF and form factor for the parton distribution of general type $`i`$. Since $`t`$-dependence is assumed to factorize and thus of no importance to the following, we will suppress it from now on. The profile functions, $`\pi ^i(x,y)`$, are asymptotic shape functions rad1  for quarks and gluons of the general form

$$\pi (x,y)=\frac{\mathrm{\Gamma }(2b+2)}{2^{2b+1}\mathrm{\Gamma }^2(b+1)}\frac{[(1|x|)^2y^2]^b}{(1|x|)^{2b+1}},$$
(5)
and normalized such that

$`{\displaystyle _{1+|x|}^{1|x|}}𝑑y\pi (x,y)=1.`$ (6)
Note that $`\pi `$ is an even function of both of its arguments footpi . The power $`b`$ controls the size of the skewing effects in the input GPD. Usually $`b=1`$ is chosen for the quarks, corresponding to maximum skewedness, whereas $`b=2`$ is chosen for the gluons. In the limit $`b\mathrm{}`$ there is no external skewedness effect, however, since in this limit, $`H(v,\xi )=q(v)`$, this translates into internal dependence on the skewedness $`(X,\zeta )=q\left(\frac{X\zeta /2}{1\zeta /2}\right)/(1\zeta /2)`$. Note that a consequence of the above model is a ratio of GPD/PDF for the quark singlet in the DGLAP region at $`X=\zeta =x_{bj}`$, which is substantially larger than $`1`$ afmmgpd ; radmus  for all experimentally relevant values of $`x_{bj}`$.

Concerning the description of the data using the above model, it was shown in afmmlong  that maximal skewing ($`b=1`$) at ‘conventional’ input scales ($`Q_0=1,2\text{GeV}`$) overshoots the H1 data by a factor of $`610`$. It was also demonstrated that one can describe, in LO only, the H1 data without including skewedness effects at the input scale, $`Q_0=2\text{GeV}`$, if one neglects evolution bemuk1 . This simplification, however, is not warranted since we know that the effects of skewed evolution are much stronger in the region of $`X\zeta `$ compared to forward evolution ffgs ; afmmgpd . Note that this region strongly influences the cross section at small $`x_{bj}`$ and some asymmetries both at small and large $`x_{bj}`$ foot3 , if they are dominated by the imaginary part of the amplitude which is proportional to $`^S(\zeta ,\zeta )`$ at LO (see for example afmmlong  as well as bemu2 ; afmmamp ; bemus ).

One may wonder whether one can come closer to the data by choosing a very low input scale and valence-like partons as in the GRV scenario grv , generating the rise of the parton distributions entirely through evolution. It turns out that choosing the canonical value of $`b_q=1`$ and GRV98 input distributions ($`Q_0=0.51,0.63`$ GeV in LO and NLO, respectively) the curves still overshoot the data considerably. Indeed, even if one tries to minimize the effect of the enhancement due to skewedness at the input scale, by choosing a large value of $`b_q=100`$, evolution still drives the prediction above the data by at least a factor of $`4`$. Since, for example, $`d\sigma _{DVCS}|^S(\zeta ,\zeta )|^2`$, the large enhancement of the quark singlet GPD at $`X=\zeta `$ is the root of the problem in this model. As we will explain below the origin of the enhancement stems from sampling singular forward sea distributions at extremely small $`x`$ in the DD-based model. To understand the last statement, one has to first establish the regions in $`x`$ in which the PDFs are sampled in the double distribution model, particularly at small $`x`$.

First, having defined the model for the factorized double distribution in Eq. (4) one may then perform the $`y^{}`$-integration in Eq. (3) using the delta function. This then modifies the limits on the $`x^{}`$-integration according to the region concerned: for the DGLAP region $`X>\zeta (v>\xi )`$ one has for the quark GPD

$`^{q,a}(X,\zeta )=`$
$`{\displaystyle \frac{2}{\zeta }}{\displaystyle _{\frac{X\zeta }{1\zeta }}^X}𝑑x^{}\pi ^q(x^{},{\displaystyle \frac{2}{\zeta }}(Xx^{})+x^{}1)q^a(x^{}).`$ (7)
For the anti-quark GPD in the DGLAP region $`X>\zeta (v<\xi )`$ one has

$`^{\overline{q},a}(X,\zeta )=`$
$`{\displaystyle \frac{2}{\zeta }}{\displaystyle _X^{\frac{X+\zeta }{1\zeta }}}𝑑x^{}\pi ^q(x^{},{\displaystyle \frac{2}{\zeta }}(X+x^{})+x^{}+1)\overline{q}^a(|x^{}|),`$ (8)
changing variables from $`xx`$ and exploiting the fact that the profile functions are even in both arguments one arrives at

$`^{\overline{q},a}(X,\zeta )=`$
$`{\displaystyle \frac{2}{\zeta }}{\displaystyle _{\frac{X\zeta }{1\zeta }}^X}𝑑x^{}\pi ^q(x^{},{\displaystyle \frac{2}{\zeta }}(Xx^{})+x^{}1)\overline{q}^a(|x^{}|),`$ (9)
so that the singlet and non-singlet quark distributions are given by

$`^S(X,\zeta )={\displaystyle \underset{a}{}}^{q,a}+^{\overline{q},a}`$
$`{\displaystyle \underset{a}{}}{\displaystyle \frac{2}{\zeta }}{\displaystyle _{\frac{X\zeta }{1\zeta }}^X}𝑑x^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))[q^a(x^{})+\overline{q}^a(x^{})],`$
$`^{NS,a}(X,\zeta )=^{q,a}^{\overline{q},a}`$
$`{\displaystyle \frac{2}{\zeta }}{\displaystyle _{\frac{X\zeta }{1\zeta }}^X}𝑑x^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))[q^a(x^{})\overline{q}^a(x^{})],`$ (10)
where $`\stackrel{~}{y}(x^{})=2(Xx^{})/\zeta +x^{}1`$.

In the ERBL region, $`X<\zeta (|v|<\xi `$) integration over $`y`$ leads to:

$`^{q,a}(X,\zeta )=`$
$`{\displaystyle \frac{2}{\zeta }}[{\displaystyle _0^X}dx^{}\pi ^q(x^{},{\displaystyle \frac{2}{\zeta }}(Xx^{})+x^{}1)q^a(x^{})`$
$`{\displaystyle _{X\zeta }^0}dx^{}\pi ^q(x^{},{\displaystyle \frac{2}{\zeta }}(Xx^{})+x^{}1)\overline{q}^a(|x^{}|)],`$
$`^{\overline{q},a}(X,\zeta )=`$
$`{\displaystyle \frac{2}{\zeta }}[{\displaystyle _0^{\zeta X}}dx^{}\pi ^q(x^{},{\displaystyle \frac{2}{\zeta }}(X+x^{})+x^{}+1)q^a(x^{})`$
$`{\displaystyle _X^0}dx^{}\pi ^q(x^{},{\displaystyle \frac{2}{\zeta }}(X+x^{})+x^{}+1)\overline{q}^a(|x^{}|)].`$ (11)
Again, using $`xx`$ and $`\pi (|x^{}|,|y^{}|)`$, one gets

$`^{q,a}(X,\zeta )`$ $`={\displaystyle \frac{2}{\zeta }}[{\displaystyle _0^X}dx^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))q^a(x^{})`$
$`{\displaystyle _0^{\zeta X}}dx^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))\overline{q}^a(x^{})],`$
$`^{\overline{q},a}(X,\zeta )`$ $`={\displaystyle \frac{2}{\zeta }}[{\displaystyle _0^{\zeta X}}dx^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))q^a(x^{})`$
$`{\displaystyle _0^X}dx^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))\overline{q}^a(x^{})].`$ (12)
Hence, for the singlet and non-singlet combinations one has

$`^S(X,\zeta )={\displaystyle \underset{a}{}}^{q,a}+^{\overline{q},a}`$
$`{\displaystyle \underset{a}{}}{\displaystyle \frac{2}{\zeta }}\left[{\displaystyle _0^X}𝑑x^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))[q^a(x^{})+\overline{q}^a(x^{})]\right]`$
$`{\displaystyle _0^{\zeta X}}dx^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))[q^a(x^{})+\overline{q}^a(x^{})]],`$
$`^{NS,a}(X,\zeta )=^{q,a}^{\overline{q},a}`$
$`{\displaystyle \frac{2}{\zeta }}[{\displaystyle _0^X}dx^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))[q^a(x^{})\overline{q}^a(x^{})]`$
$`+{\displaystyle _0^{\zeta X}}dx^{}\pi ^q(x^{},\stackrel{~}{y}(x^{}))[q^a(x^{})\overline{q}^a(x^{})]].`$ (13)
These expressions clearly satisfy the correct symmetries properties, i.e., $`F^S(\zeta X,\zeta )=F^S(X,\zeta )`$, $`F^{NS,a}(\zeta X,\zeta )=F^{NS,a}(X,\zeta )`$ (n.b. $`\stackrel{~}{y}(x^{})\stackrel{~}{y}(x^{})`$ when $`X\zeta X`$ ). Analogously for the gluon one obtains

$`^g(X,\zeta )={\displaystyle \frac{2}{\zeta }}{\displaystyle _{\frac{X\zeta }{1\zeta }}^X}𝑑x^{}\pi ^g(x^{},\stackrel{~}{y}(x^{}))x^{}g(x^{}),`$ (14)
for the DGLAP region, and

$`^g(X,\zeta )`$ $`={\displaystyle \frac{2}{\zeta }}[{\displaystyle _0^X}dx^{}\pi ^g(x^{},\stackrel{~}{y}(x^{}))x^{}g(x^{})`$
$`+{\displaystyle _0^{\zeta X}}dx^{}\pi ^g(x^{},\stackrel{~}{y}(x^{}))x^{}g(x^{})],`$ (15)
for the ERBL region (which is symmetric under $`X\zeta X`$).

Inspection of the integration limits in Eqs. (10, 13) highlights the main problem. In the limit $`X\zeta `$, as a result of the lower limits of the integrals the forward PDF is sampled closer and closer to $`x^{}=0`$, where it has not yet been measured. This will be irrelevant providing the integrand is sufficiently non-singular in $`x^{}`$ in this region i.e. it can happen if the profile functions, $`\pi ^i`$, provide a strong suppression of this region, or if the PDFs themselves are sufficiently non-singular. However, we know that phenomenological quark and sometimes even gluon input distributions are singular in the small $`x`$ region. In the quark case this problem is made worse by the fact that we sample the number distribution $`q(x^{})`$ rather than the momentum distribution $`x^{}q(x^{})`$ (so that a non-singular momentum distribution $`xq(x)x^a`$ for $`a[0,1]`$ will give a singular number distribution $`q(x)x^{a1}`$).

It turns out that for realistic quark distributions the region close is $`x^{}=0`$ is very significantly sampled for small $`\zeta =x_{bj}`$. This leads to two serious problems. Firstly, the forward distributions are unknown here so one must extrapolate the ‘known’ analytic forms downwards in $`x^{}`$. Secondly, and much more importantly, it leads to a very significant enhancement of the quark singlet GPDs relative to the PDFs for $`X\zeta `$, i.e., the region most relevant for DVCS. Though of paramount importance for DVCS, this region is but a small region of phase space where the current factorized DD models fail.

We illustrate this using a series of three figures relating to the formation of the quark singlet GPD in the DGLAP region close to $`X=\zeta `$. Fig. 3 shows the integrand, $`I(x^{})=\pi ^u(x^{},\stackrel{~}{y})u^S(x^{})`$, of Eq.(10) for the up quark singlet (multiplied by $`\zeta `$) as a function of $`x^{}/\zeta `$ for two values of $`\zeta =0.1,0.0001`$ and two values of $`X\zeta =0.1,0.001`$. Clearly as $`X`$ approaches $`\zeta `$ the PDF is sampled at progressively smaller values of $`x^{}\zeta `$, where for small $`\zeta `$ it is unknown. Fig. 4 shows the average value of $`x^{}`$ sampled in this integral (divided by $`\zeta `$) as a function of $`\zeta `$ for several values of $`X\zeta `$. For very small values of $`X\zeta `$ the average value of $`x^{}`$ settles down to about $`\zeta /4`$, for small $`\zeta `$. Finally in Fig. 5 we show, for MRST input PDFs, the ratio of the quark singlet GPD to PDF at $`\zeta =0.0001,0.1`$, for the canonical value of the power, $`b_q=1`$, in Eq. (5). Note the large enhancement of the GPD at $`X\zeta `$, particularly for small $`\zeta `$ in the upper plot. We emphasize that this enhancement, which leads to an overshoot of the DVCS data, is built in right at the start in the modeling of the quark singlet GPD at the input scale. One also sees that for the gluon, which uses $`x^{}g(x^{})`$, and $`b_g=2`$, the ratio remains close to unity.

The most important enhancement effect in the valence region, $`\zeta \text{ }>0.1`$, originates from the relative shift of the parton momentum fraction $`X`$ to smaller values close to $`X=\zeta `$ (although the enhancement from small $`x^{^{}}`$ is still significant). As we will show in the next section, the assumption that $`H(v,\xi )=(1\zeta /2)(X,\zeta )=q(v)`$ with $`v=(X\zeta /2)(1\zeta /2)`$ gives a good description of the data at both small and large $`x_{bj}`$. As stated before, this corresponds to a factorized DD-model with $`b=\mathrm{}`$, i.e., with no external skewedness. However, in terms of a comparison of GPD to forward PDF, there is a residual effect of skewedness since one now has to compare $`q(v)`$ with $`q(X)`$. Since we are comparing number distributions which are more singular than momentum distributions, any shift in the momentum fraction to smaller values will lead to a quite a large enhancement of $`q(v)`$ relative to $`q(X)`$. For CTEQ6M for example the enhancement at $`X=0.1`$ and $`\zeta =0.1`$ is about $`1.7`$ for the quark singlet, which increases further if more skewedness is added by decreasing $`b`$.

However, as we will demonstrate in Section III, the available data allows little room for further enhancement due to skewedness at the input scale since the LO result, at least, is already close to the upper bound of the experimental errors. Therefore, only the extremal “$`b=\mathrm{}`$” version of the current factorized DD-model can be used to describe the data. An obvious solution to this is to modify the quark singlet profile functions in Eq.(4) in such a way as to suppress the region of very small $`x^{}`$. However, one must find a new functional form which achieves this without spoiling the known mathematical features of GPDs discussed above. Alternatively, one abandons the attempt to model a DD using a factorized form $`\pi (x,y)f(x)`$, though appealing due to its simplicity but possibly too simplistic in its form. This remains an open problem and has to be addressed by those who wish to use the double distribution framework, factorized or not, to model GPDs.

## III The forward input model and the aligned jet model

In this section we revisit the logic of setting the GPDs equal to the forward PDFs by proposing an alternative forward model to that suggested in gb2 , with suitably-symmetrized input GPDs in the ERBL region constructed so as to satisfy the requirements of polynomiality for the first two moments.

In ffs  DVCS was predicted to be measurable at DESY-HERA and, allowing for the freedom associated with choosing the slope parameter, $`B`$, the predictions successfully describe both the H1 data h1  and the recent ZEUS result zeus ; zeus2  on the photon-level DVCS cross section. This was achieved by modeling the imaginary part of the DVCS amplitude at the input scale using the aligned jet model (AJM) fs88 . This was then compared to the imaginary part of the DIS amplitude, calculated within the same framework, which was found to be smaller by a factor of about two. The comparison enabled the normalization of the DVCS amplitude at the input scale to be set using $`F_2`$ structure function data. The DVCS amplitude was then evolved to higher scales using LO skewed evolution in perturbative QCD.

The basic relation between the DVCS and DIS amplitudes, using the AJM, is given by ffs

$`R={\displaystyle \frac{\text{Im}𝒯_{\text{DVCS}}}{\text{Im}𝒯_{\text{DIS}}}}=\mathrm{ln}\left(1+{\displaystyle \frac{Q^2}{M_0^2}}\right)\left(1+{\displaystyle \frac{M_0^2}{Q^2}}\right)1.52.5,`$ (16)
where $`Q^2`$ for the AJM is typically $`13`$ GeV<sup>2</sup> and $`M_0`$ is a hadronic scale which roughly corresponds to the lowest allowed, excited intermediate state in the $`s`$-channel. Therefore, $`M_0^20.40.6\text{GeV}^2`$, or about $`m_\rho ^2`$. The AJM neglects the contribution of quarks with large transverse momenta in the quark loop attached to the photons in the handbag diagram. Since the contribution of small transverse momenta is more symmetric than the one at large transverse momenta, the AJM may somewhat overestimate the effect of skewedness at the input scale. Eq. (16) can be generalized to demonstrate how the forward limit $`\text{Im}𝒯_{\text{DVCS}}=\text{Im}𝒯_{\text{DIS}}`$ is achieved, i.e., how the skewedness effect is reduced by giving the outgoing photon a space-like virtuality, $`q^{}_{}{}^{}2=Q^{}_{}{}^{}2`$:

$`R={\displaystyle \frac{\text{Im}𝒯_{\text{DVCS}}}{\text{Im}𝒯_{\text{DIS}}}}=\mathrm{ln}\left({\displaystyle \frac{1+\frac{Q^2}{M_0^2}}{1+\frac{Q^{}_{}{}^{}2}{M_0^2}}}\right){\displaystyle \frac{1+\frac{M_0^2}{Q^2}}{1\frac{Q^{}_{}{}^{}2}{Q^2}}}.`$ (17)
This procedure allows one to derive a very important relation between the relative momentum fractions of the outgoing and returning partons, $`X`$ and $`(X\zeta )`$ of the quark singlet GPD, and the virtualities of the incoming and outgoing photons

$$\lambda =\frac{X\zeta }{X}=\frac{Q^{}_{}{}^{}2}{Q^2}.$$
(18)
We illustrate this point in Fig. 6 by plotting $`R`$ as a function of $`\lambda `$ for several values of $`Q^2`$ and two values of $`M_0^2`$ to demonstrate the relative insensitivity of $`R`$ (within $`2030\%`$) to $`M_0`$.

The plot shows that as $`X`$ increases relative to $`\zeta `$ in the DGLAP region the ratio drops rapidly to its forward limit. For example, at $`\lambda =1/2`$, i.e., $`X=2\zeta `$, the curves are very flat and there is only a modest enhancement of $`2040\%`$. One also encounters this behavior in the DD model if one investigates the ratio of the GPD to the PDF in the DGLAP region (see for example afmmgpd ; radmus ). It would therefore be advantageous to be able to directly relate $`R`$ to a ratio of GPD to PDF. Trusting that perturbative QCD is applicable at the AJM input scale one can, in LO at least where the coefficient function is trivial, directly translate the ratio in amplitudes for a particular $`\lambda `$ into a ratio of GPD to PDF:

$`R(\lambda )`$ $`={\displaystyle \frac{\text{Im}𝒯_{\text{DVCS}}}{\text{Im}𝒯_{\text{DIS}}}}={\displaystyle \frac{H^S(v=\xi \frac{1+\lambda }{1\lambda },\xi )}{q^S(X)}},`$
$`={\displaystyle \frac{^S\left(X=\zeta /(1\lambda ),\zeta \right)}{(1\zeta /2)q^S(X)}},`$ (19)
i.e.,

$$^S(\zeta /(1\lambda ),\zeta )=(1\zeta /2)R(\lambda )q^S(X).$$
(20)
There are several comments in order at this point: $`\lambda `$ is now bounded from above through $`\lambda 1\zeta `$. This implies that the relationship between the ratios in Eq. (20) is only strictly true for $`\lambda =0`$ (i.e., for DVCS for which $`Q^{}_{}{}^{}2=0`$). The case $`\lambda 0`$ should be viewed as follows: for $`Q^22\text{GeV}^2`$, there is still the possibility of having more than one rung in the partonic ladder. Probing the uppermost rung with a virtuality $`Q^2`$ reveals the distribution in momentum fractions, in this case $`X=\zeta `$ i.e., $`X\zeta =0`$ corresponding to $`\lambda =0`$. The next rung and its distribution in momentum fractions can be accessed by ‘emitting’ a photon with space-like virtuality (i.e., $`Q^{}_{}{}^{}2>0`$). As $`Q^{}_{}{}^{}2`$ increases and one goes further down the ladder to where $`X\zeta `$, one approaches the forward limit. If one keeps the interpretation of the $`s`$-channel cut as being equal to the imaginary part of the ‘scattering’ amplitude for $`Q^{}_{}{}^{}20`$, which, in LO, is directly proportional to GPD/PDF at $`X\zeta `$ rather than at $`X=\zeta `$, then the ratio $`R`$ is a direct measure of quark singlet GPD to the quark singlet PDF for $`X\zeta `$. However, this logic is valid only at LO. At NLO, the situation radically changes since, first, the gluon directly enters into the amplitude, secondly, the convolution of the coefficient function with the GPD is no longer as trivial as in LO and thirdly, the value of $`\alpha _s`$ at low $`Q^2`$ is quite different in LO and NLO. Therefore such a simple relation as in Eq. (20) should and can no longer be valid. In order to keep the ansatz as simple as possible, we will only require that the model GPDs, at least in LO, produce a ratio, $`R`$, which is in broad agreement with the $`R`$ values obtained in the AJM from Eq. (20).

If one chooses the forward model ansatz where the GPD equals the PDF at $`v`$ in both LO and NLO, due to a lack of a better Ansatz, (see e.g. gb2 ):

$$H^S(v,\xi )=q^S(v)q^S\left(\frac{X\zeta /2}{1\zeta /2}\right),$$
(21)
which corresponds to the $`b\mathrm{}`$ limit of the DD model, one obtains a ratio of GPD to PDF at $`X=\zeta `$ of $`2.3`$ for the quark singlet, in agreement with the AJM prediction though for slightly different values of $`Q_0^2`$ (directly compare the upper line in the upper plot of Fig. 6 with the quark singlet in Fig. 7 keeping in mind that $`\lambda =0X/\zeta =1`$ and $`\lambda =0.9X/\zeta =10`$). In consequence, our model ansatz corresponds to an AJM with maximal skewedness.

For our forward model in the DGLAP region we, therefore, choose for simplicity, the ansatz of Eq. (21) for the quark singlet, the non-singlet (i.e., the valence) and the gluon. The ratios of GPD to PDF at the input scale for MRST01 quark singlet and gluon distributions mrst01  at LO and NLO at $`Q_0=1\text{GeV}`$ are shown in Fig. 7. For a function which falls as $`x`$ decreases, such as the valence quark at small $`x`$ or the MRST gluon at LO, this ansatz leads to a suppression of the point $`X=\zeta `$ relative to the forward case (see the lower, dotted line in Fig. 7). Note that in NLO the MRST gluon actually goes negative at small $`x`$, so a ratio of GPD/PDF $`>1`$ close to $`X=\zeta `$ in this case leads to a suppression of the DVCS cross section from the gluon contribution, relative to using PDFs. However, the DVCS cross section is rather insensitive to the behavior of this ratio close to $`X=\zeta `$ for the gluon, since it only enters in NLO, and is completely insensitive to it for the non-singlet quark case since this distribution only enters into the evolution.

The above reasoning indicates that the physics of the AJM model provides a guide for modeling input GPDs in the DGLAP region. If one compares the NLO imaginary part of the DVCS amplitude from the above model to the NLO imaginary part of the DIS amplitude as extracted from a recent H1 $`F_2`$ fit H1f2fit , we find for MRST2001, $`R3.8`$ and for CTEQ6M, $`R2.7`$ for $`Q^2=3\text{GeV}^2`$ and $`x=0.0005`$. the enhancement effect of evolution, the AJM result for $`R`$ is basically reproduced for CTEQ6M but not for MRST2001 using eq. (21). Given the widely different parameterizations at NLO, this seems acceptable to us, at present. As we will see below the enhancement effect generated through the shift is too strong for both LO and NLO at low values of $`Q^2`$ near the input scale.

The prescription in eq. (21) does not dictate what to do in the ERBL region, which does not have a forward analog. Naturally the GPDs should be continuous through the point $`X=\zeta `$ and should have the correct symmetries around the midpoint of the ERBL region. They are also required to satisfy the requirements of polynomiality:

$`M_N`$ $`={\displaystyle _1^1}𝑑vv^{N1}[H^q(v,\xi )+H^g(v,\xi )]`$
$`={\displaystyle \underset{k=0}{\overset{N/2}{}}}\xi ^{2k}C_{2k,N}.`$ (22)
At this point we choose to model the ERBL region with these natural features in mind. We demand that the resultant GPDs reproduce the first moment $`M_13`$ and the second moment $`M_21`$ footmom  (with the D-term set to zero to remove the quadratic piece in Eq. (24)) where

$`M_1`$ $`=`$ $`{\displaystyle _0^1}𝑑vH^{q,NS}(v,\xi )=3,`$ (23)
$`M_2`$ $`=`$ $`{\displaystyle _0^1}𝑑vv[H^{q,S}(v,\xi )+H^g(v,\xi )]`$ (24)
$`=`$ $`1+C\xi ^2`$
and $`C`$ in Eq. (24) was computed in the chiral-quark-soliton model vanderhagen1 . This reasoning suggests the following simple analytical form for the ERBL region ($`X<\zeta `$):

$`^{g,NS}(X,\zeta )=^{g,NS}(\zeta )\left[1+A^{g,NS}(\zeta )C^{g,NS}(X,\zeta )\right],`$
$`^S(X,\zeta )=^S(\zeta )\left({\displaystyle \frac{X\zeta /2}{\zeta /2}}\right)\left[1+A^S(\zeta )C^S(X,\zeta )\right],`$ (25)
where the functions

$`C^{g,NS}(X,\zeta )={\displaystyle \frac{3}{2}}{\displaystyle \frac{2\zeta }{\zeta }}\left(1\left({\displaystyle \frac{X\zeta /2}{\zeta /2}}\right)^2\right),`$
$`C^S(X,\zeta )={\displaystyle \frac{15}{2}}\left({\displaystyle \frac{2\zeta }{\zeta }}\right)^2\left(1\left({\displaystyle \frac{X\zeta /2}{\zeta /2}}\right)^2\right),`$ (26)
vanish at $`X=\zeta `$ to guarantee continuity of the GPDs. The $`A^i(\zeta )`$ are then calculated for each $`\zeta `$ by demanding that the first two moments of the GPDs are explicitly satisfied (remembering to include the D-term in the ERBL region which only provides the quadratic term in $`\xi `$ in Eq. (23)). For the second moment what we do in practice is to set the D-term to zero and demand that for each flavor the whole integral over the GPD is equal to the whole integral over the forward PDF for the input distribution concerned (due to the inherent small errors on the PDFs, the sum of such integrals will be close to, but not precisely equal to, unity). Note that the modeling of the ERBL region is unimportant at small $`\zeta `$ since the unknown subtraction constant in the dispersion relation between the real part of the amplitude which formally depends on both the ERBL and the DGLAP region, and the imaginary part which formally depends only on the DGLAP region, is proportional to $`\zeta `$ and therefore inconsequential at small $`\zeta `$.

It would be straightforward to extend this algorithm to satisfy polynomiality to arbitrary accuracy by writing the $`A^i(\zeta )`$ explicitly as a polynomial in $`\zeta `$ where the first few coefficients are set by the first two moments and the other coefficients are then randomly chosen since nothing is known about them. The above algorithm is extremely flexible both in its implementation and adaption to either other forward PDFs or other functional forms in the ERBL region. Therefore it can be easily incorporated into a fitting procedure.

In Fig. 8 we show the shape of the resulting input GPDs for two characteristic values of $`\zeta =0.001,0.1`$. The upper plot in this figure explicitly shows the antisymmetry of the singlet GPD and the symmetry of the gluon GPD about the point $`X=\zeta /2`$.

The photon level cross section results from this model, using MRST01 mrst01  and CTEQ6 cteq6  input distributions at LO and NLO, are compared in Figs. 9,10 to the H1 h1  and ZEUS zeus2  data at their average kinematic points, respectively. In these curves we chose to use an $`x`$ and $`Q^2`$-independent slope parameter of $`B=6.5`$ GeV<sup>-2</sup>, but realistically there is a $`3040\%`$ uncertainty associated with the value of this unknown parameter. The figures illustrate that within the framework of the forward input model for GPDs the DVCS cross section remains rather sensitive to the choice of input PDF and to the accuracy with which the calculation is performed (i.e., LO or NLO). It should be noted that the description at lowest H1 value of $`Q^2`$ is bad. However, the enhancement we chose using eq. (21) at $`Q_0^2=11.69\text{GeV}^2`$ corresponds more to the AJM at $`Q^2=3\text{GeV}^2`$ (compare Figs. 6,7). Hence it is not surprising that the description at low $`Q^2`$ is not good suggesting that the shift in $`X`$ in eq. (21) should be less at lower values of $`Q_0`$ . When increasing the input scale of CTEQ6M, as done below, we find an appropriate reduction in the cross section at low $`Q^2`$ much more in line with the low $`Q^2`$ data and the AJM value.

It is important to note that the preliminary ZEUS data lies systematically above the H1 data (see Fig. 11 of zeus2 ). Overall NLO seems to be doing better than LO, particularly on the slope of the energy dependence. It is fair to say that all of the theory curves appear to have a $`Q^2`$-dependence that is too steep to describe all of the data. We will return to this point in the next section.

The difference between the MRST and CTEQ curves at LO and NLO reflects the relative size of the quark singlet and gluon distributions for each set. It is possible that more precise data on DVCS may eventually allow a discrimination between various input scenarios using NLO QCD. For this to be realistic one would first need to pin down the uncertainty associated with the slope by explicitly measuring the $`t`$-dependence.

We also investigated the effect on the cross section of increasing the input scale for skewed evolution using CTEQ input distributions, from the starting scale $`Q_0=1.3`$ GeV to $`Q_0=2.0`$ GeV. We then use the forward PDFs at the new scale in our model for the GPDs. Fig. 11 shows that the reduced lever arm for skewed evolution starting at the higher scale leads to a smaller cross section at LO and NLO, as expected, and that, in LO at least, the effect of this change is rather large. In fact the CTEQ and MRST collaborations only advocate the use of their forward PDFs above $`Q^234`$ GeV<sup>2</sup> (they start evolution at a lower scale $`Q^212`$ GeV<sup>2</sup> due to technicalities associated with a consistent implementation of charm). Hence, it is not completely clear where one should start skewed evolution, and this constitutes an additional uncertainty in the theoretical predictions.

Having compared to small $`x_{bj}`$ data we will now test the AJM ansatz for large $`x_{bj}`$ by comparing to data on the single-spin asymmetry, SSA, (HERMES and CLAS hermclas ) and the charge asymmetry, CA, (HERMES only), defined by

$`SSA={\displaystyle \frac{2_0^{2\pi }𝑑\varphi \mathrm{sin}(\varphi )(d\sigma ^{}d\sigma ^{})}{_0^{2\pi }𝑑\varphi (d\sigma ^{}+d\sigma ^{})}},`$
$`CA={\displaystyle \frac{2_0^{2\pi }𝑑\varphi \mathrm{cos}(\varphi )(d\sigma ^+d\sigma ^{})}{_0^{2\pi }𝑑\varphi (d\sigma ^++d\sigma ^{})}}.`$ (27)
Here $`d\sigma ^{}`$, $`d\sigma ^{}`$ refer to the differential cross sections with the lepton polarized along or against its direction of motion, respectively; $`d\sigma ^+,d\sigma ^{}`$ are the unpolarized differential cross sections for positrons and electrons, respectively.

Such a comparison of QCD models with the available high $`x_{bj}`$ data may be viewed with some scepticism, especially in the case of CLAS which has such a low $`Q^212\text{GeV}^2`$ (HERMES is only slightly better with a typical $`Q^2`$ of $`24\text{GeV}^2`$). Firstly, it is a priori not clear that perturbation theory is applicable at such low $`Q^2`$ values, (in particular, higher twist corrections may be expected to become important in this region and our approximations correspond to the DVCS cross section being divergent as $`Q^20`$.). Secondly, the previously neglected GPDs $`\stackrel{~}{H},E`$ and $`\stackrel{~}{E}`$ become increasingly important as $`x`$ increases afmmlong ; bemu3 . In the following we will include the dominant twist-3 contributions bemu4 , which are entirely kinematic in origin, in our calculation of the differential cross section, neglecting the sub-dominant twist-3 effects. We use the same input models for $`\stackrel{~}{H},E`$ and $`\stackrel{~}{E}`$ as well as $`t`$-dependence for the various amplitudes i.e. various dipole form factors which are about equivalent with an exponential at small $`t`$ with a slope $`B`$ between $`5`$ and $`8`$, as in afmmlong . For HERMES we perform a full LO and NLO QCD analysis, whereas for CLAS we are restricting ourselves to LO, i.e., we are testing handbag dominance with no or little evolution. Furthermore, we shall restrict ourselves to MRST01 input PDF for simplicity.

It transpires that the average kinematics of HERMES is such that $`H`$ is still the leading GPD and within our model assumptions $`\stackrel{~}{H},E`$ and $`\stackrel{~}{E}`$ could be set to zero for those values, with negligible difference to the final answer foot7b . Within the above caveats, we find for average HERMES kinematics ($`x=0.11,Q^2=2.56\text{GeV}^2,t=0.265\text{GeV}^2`$)

SSA $`=0.28\text{(LO)},0.23\text{(NLO)}`$
CA $`=0.12\text{(LO)},0.09\text{(NLO)}`$ (28)
compared to the quoted experimental results hermnew

SSA $`=0.21\pm 0.08`$
CA $`=0.11\pm 0.07.`$ (29)
For the average CLAS kinematics ($`x=0.19,Q^2=1.31\text{GeV}^2,t=0.19\text{GeV}^2`$) we find

$$\text{SSA}=0.2\text{(LO)}$$
(30)
compared to the experimental value (second reference of hermclas )

$$\text{SSA}=0.202\pm 0.041.$$
(31)
This demonstrates that the AJM ansatz works surprisingly well even at large $`x_{bj}`$ giving us confidence in the AJM-based model and suggesting that a fit to the available data should be possible without tuning too many input parameters.

## IV A simple model for the slope parameter

It was pointed out in ffs  that the $`t`$-slope of the DVCS cross section at small $`x`$ should depend strongly on $`Q^2`$ in the transitional region from $`Q^2`$ of a few GeV<sup>2</sup> to large $`Q^2`$. At $`Q^22`$ GeV<sup>2</sup> it is natural to expect that the slope will be pretty close to that for exclusive $`\rho `$-meson production: $`B8`$ GeV<sup>-2</sup> brho . For large $`Q^2`$ the dominant contribution is governed by evolution trajectories which, at the resolution $`Q_0^22`$ GeV<sup>2</sup>, originate from the gluon field. Hence we expect that in this case the slope will be given by the square of the two-gluon form factor of the nucleon at $`X,X\zeta x_{bj}`$. Recently fs  it was demonstrated that for $`x_{bj}0.05`$ this $`t`$-dependence can be approximated in a wide range of $`t`$ as $`1/(1t/m_{2g}^2)^4`$ with $`m_{2g}^21.1`$ GeV<sup>2</sup>. This corresponds to a $`t`$-slope of $`B3`$ GeV<sup>-2</sup> for exponential fits foot8 . At smaller $`x`$ an increase of the slope is expected which could originate from several effects, including Gribov diffusion. Hence for the highest $`Q^2`$ point of ZEUS of about $`90`$ GeV<sup>2</sup> we expect $`B=3.5\pm 0.5`$ GeV<sup>-2</sup>. The recent H1 and ZEUS $`\rho `$-meson production data, for a $`W`$-range similar to the DVCS experiments, indicates that the slope of $`\rho `$-production both for $`\sigma _L`$ and $`\sigma _T`$ drops rather rapidly with increasing $`Q^2`$ reaching $`B5`$ at $`Q^210`$ GeV<sup>2</sup> brho2 .

A simple parameterization which reflects the discussed constraints for the range of $`2Q^2100`$ GeV<sup>2</sup> is

$$B(Q^2)=B_0\left(1C\mathrm{ln}\left(\frac{Q^2}{Q_0^2}\right)\right)$$
(32)
with $`B_0=8`$ GeV<sup>-2</sup>, $`Q_0=2`$ GeV<sup>2</sup>, $`C=0.15`$ being reasonable values for the various parameters. This gives $`B(Q^2=9.6)=6.1`$ GeV<sup>-2</sup> and $`B(Q^2=4.5)=7.0`$ GeV<sup>-2</sup> at the average $`Q^2`$ values of the ZEUS and H1 data, respectively (in broad agreement with our chosen constant value of $`B=6.5`$ GeV<sup>-2</sup>). Fig. 12 illustrates the effect of introducing this simple model on the description of the $`Q^2`$-dependence of the ZEUS data.

This modification of the $`B`$ slope gives a great improvement in comparison with the data and shows how important an experimental determination of the $`B`$ slope is, since it constitutes a large theoretical uncertainty at this point.

## V Conclusions

Using sound phenomenological guidelines such as the aligned jet model, we have constructed a model for generalized parton distributions at the input scale. Within certain theoretical uncertainties (such as the exact shape, the input scale and the functional form of $`B`$ in $`W`$ and $`Q^2`$) this model can be used in a NLO QCD analysis to describe the recent DVCS data from H1, ZEUS, HERMES and CLAS within their experimental errors. In constructing this model we have given a simple and flexible algorithm which can be easily incorporated into a fitting procedure.

We have also demonstrated and explained the failure of the most widely used model for generalized parton distributions, the factorized double distribution based model, to describe the available DVCS data, when rigorously applied in a LO or NLO QCD analysis in its current form.

The modeling of the input GPDs is now sufficiently advanced to justify attempting to fit some of the input parameters directly to the available data. A basic analysis of the data would seem to favor a $`t`$-dependence with a slope parameter, $`B`$, that depends on $`Q^2`$. Hence, an accurate measurement of this slope is of crucial importance for further progress of the comparison of theory and experiment.

We would like to thank H. Abramowicz, C. Weiss, M. Diehl, L. Frankfurt, K. Golec-Biernat, D. Müller and A. Schäfer for helpful discussions. A.F. was supported by the DFG Emmi-Noether Program, M.M. by PPARC and M.S. by the DOE. This work is dedicated to Dina, Joanne and Patricia.