Sum rule for the Compton amplitude and implications for the proton-neutron mass difference

The Cottingham formula expresses the leading contribution of the electromagnetic interaction to the proton-neutron mass difference as an integral over the forward Compton amplitude. Since quarks and gluons Reggeize, the dispersive representation of this amplitude requires a subtraction. We assume that the asymptotic behaviour is dominated by Reggeon exchange. This leads to a sum rule that expresses the subtraction function in terms of measurable quantities. As shown in our previous paper, the evaluation of this sum rule leads to $m_{QED}^{p-n}=0.58\pm 0.16\,\mbox{MeV}$. The present article describes the underlying analysis in detail.


Introduction
In the framework of the Standard Model, the fact that proton and neutron have nearly the same mass is explained as consequence of an approximate symmetry: isospin [2]. The symmetry is broken explicitly because the two lightest quarks neither have the same charge nor the same mass. The violation of the symmetry is very weak, because the e.m. coupling e as well as the difference between m u and m d are small. The weak interaction provides the neutron mass with an imaginary part and generates a shift of the real part as well, but these effects are tiny and will be neglected. The Standard Model then reduces to QED+QCD. In that framework, the expansion of the mass difference between proton and neutron in powers of e starts with where m QCD is what remains if e is turned off and is proportional to m u − m d , while m QED stands for the term of order e 2 . It is well-known that the splitting of the physical masses into an electromagnetic and a strong part is not unique. In our analysis, the ambiguity shows up through the scale of the logarithm occurring in the e.m. renormalization of the quark masses and will be discussed in detail.
As shown by Cottingham [3], the leading contribution of the e.m. interaction to the mass of a particle is given by an integral over the spin averaged forward Compton scattering amplitude, This amplitude is determined by QCD. If the mass of a particle is expanded in powers of the e.m. coupling constant e, the explicit expression for the term of order e 2 is formally given by 1 There are two problems with this formula: (i) The short distance properties of QCD imply that the amplitude T µ µ (p, q) does not fall off rapidly enough at large values of q for the integral to converge. (ii) T µ µ (p, q) does not obey an unsubtracted dispersion relation -causality alone determines the Compton amplitude through the structure functions of lepton-nucleon scattering only up to a subtraction function. The asymptotic behaviour in the deep inelastic region is now fully understood on the basis of asymptotic freedom, but the properties of the subtraction function are still under debate.
Elitzur and Harari [4] pointed out that if the exchange of Reggeons correctly describes the asymptotic behaviour in the limit ν → ∞ at fixed q 2 -an assumption we refer to as Reggeon dominance -then the subtraction function obeys a sum rule which fully determines it through the cross section of lepton-nucleon scattering. Their paper appeared in 1970, at a time when the origin of the ∆I = 1 mass differences within the isospin multiplets was totally mysterious: evaluations of the Cottingham formula invariably led to the conclusion that the proton should be heavier than the neutron and hence unstable.
In 1975 Gasser and Leutwyler [5] then showed that the mystery disappears if the popular conviction, according to which the strong interaction conserves isospin, is dismissed. They showed that a coherent picture of isospin breaking can be reached within the Quark Model, provided the masses of the two lightest quarks are not only very small but also very different. 1 We fix the normalization of the one-particle states with p ′ , s ′ |p, s = 2p 0 (2π) 3 δ 3 (p ′ −p)δ s ′ s and m is the mass of the particle. The spin averaged matrix element of the operator O is abbreviated with p|O|p ≡ 1 2 s p, s|O|p, s .
At that time, the experimental results on deep inelastic scattering were consistent with the scaling laws of Bjorken [6]. The implications of Reggeon dominance were worked out in this framework, using models to substitute the lack of experimental information in part of phase space, with the result m QED = 0.7±0.3 MeV [5].
Lattice calculations of the proton-neutron mass difference are very demanding and became feasible only in the 21st century. Early calculations were consistent with the result obtained from the Cottingham formula, but more recent evaluations indicate higher values for m QED -we will compare the available results with the outcome of our calculation in section 23. In 2012, Walker-Loud, Carlson and Miller [7] performed a new evaluation of the Cottingham formula. They claimed that the analysis in [5] is inconsistent and replaced our sum rule by a model where the subtraction function T 1 (0, q 2 ) is parametrized with a simple algebraic formula. This paper triggered renewed interest and several authors investigated the matter [8][9][10][11]. We will discuss these works in section 24. A critical examination of some of the claims made in [7] can be found in appendix E of [12] and in [13,14].
The discovery of QCD and asymptotic freedom led to a fully transparent picture for the properties of the Compton amplitude in the region where both ν and q 2 are large and where the divergence of the Cottingham formula arises [15][16][17][18]. In the letter [1], we showed that the formal relation (3) can be rewritten in such a way that the divergences are under full theoretical control, exclusively concern the contribution from the subtraction function and are absorbed in the e.m. renormalization of quark masses and QCD coupling constant. The aim of the present paper is to describe the analysis underlying these statements in detail.
The presentation is organized as follows. In a first part, sections 2-7, we discuss the mathematical underpinnings: decomposition of the Compton amplitude, dispersion relations, sum rule for the subtraction function, Wick rotation, mass formulae. The second part, sections 8-13, deals with the operator product expansion, which governs the behaviour of the amplitudes at large momenta. The renormalization of the mass difference is discussed in sections 14 and 15, whereas the data concerning the structure functions used in our work and the numerical determination of the subtraction function and of the mass difference are described in sections [16][17][18][19][20][21][22]. Sections 23 and 24 compare the outcome of our analysis with results obtained on the lattice and with other recent evaluations of the Cottingham formula. A summary and conclusions are provided in section 25. The Appendices contain material concerning the operator product expansion as well as a detailed derivation of the sum rule for the subtraction function that plays a central role in the present work.

Lorentz invariance, kinematic zeros
Causality ensures that the time-ordered amplitude is unique up to contact terms and the ambiguity can be fixed in such a manner that T µν (p, q) is Lorentz covariant. 2 Together with symmetry under space reflections, this property implies that the Compton amplitude can be decomposed as T µν (p, q) = A g µν + B p µ p ν + Cp µ q ν + C ′ p ν q µ + D q µ q ν , where A, B, C, C ′ , D only depend on the two variables q 2 and ν = p · q/m. Current conservation imposes the constraints A + mν C + q 2 D = 0 , mν B + q 2 C = 0 , C ′ = C . (5) Since the physical spectrum of QCD does not contain massless particles, the amplitude T µν (p, q) cannot have a pole at q 2 = 0. Hence the second relation shows that B vanishes for q 2 = 0 and can therefore be represented as B = −q 2 T 2 /m 2 . Setting T 1 = D and solving the constraints (5) for A and C, this leads to the decomposition T µν (p, q) = T 1 (ν, q 2 )K µν 1 + T 2 (ν, q 2 )K µν 2 , K µν 1 = q µ q ν − g µν q 2 , K µν 2 = 1 m 2 {(p µ q ν + p ν q µ )p · q − g µν (p · q) 2 − p µ p ν q 2 }. Crossing symmetry, T µν (p, q) = T νµ (p, −q), implies that T 1 and T 2 are even in ν.
A popular alternative decomposition identifies the two independent amplitudes instead withT 1 = −A and T 2 = m 2 B. It is related to the one specified above bŷ The problem with this choice is that, in contrast to T 1 , T 2 , the amplitudesT 1 ,T 2 contain kinematic zeros. This makes it difficult to determine their asymptotic behaviour. That is important because analytic functions are fully determined by their singularities only if the asymptotic behaviour is known. In dispersion theory, theoretical constraints are needed to determine the asymptotic properties of the amplitudes.
To illustrate the problems encountered when working with amplitudes that are not free of kinematic zeros, consider the "Born terms", i.e. the poles generated by the one-particle intermediate states. Their residues are determined by the elastic form factors of the nucleon. The Cauchy formula implies that an analytic function of the variable z is determined uniquely by its singularities (poles, cuts) and by its behaviour for z → ∞. The amplitudes T 1 and T 2 are analytic in ν at fixed q 2 . The Born terms concern the contributions from the nucleon poles at ν = ± q 2 /2m. They are fixed uniquely by the requirement that they disappear for ν → ∞ [12]: .
For notation, in particular also for the definition of the Sachs form factors G E and G M , we refer to [12]. For the alternative decomposition (7), the elastic part ofT 1 does not disappear when ν → ∞. In terms of Regge poles, the elastic part ofT 1 contains a fixed pole at α = 0, with a residue that is determined by the nucleon form factors:T 1 picks up asymptotic contributions that do not have anything to do with the phenomena that dominate the high energy behaviour of the Compton amplitude -they merely reflect the fact that the amplitudeT 1 contains kinematic zeros. It is not advisable to work with such amplitudes -for further discussion of the problems encountered in the presence of kinematic zeros, we refer to [14,19,20].
As pointed out in the letter [1], the operator product expansion shows that, up to normalization, the leading spin 2 contributions in T 1 and T 2 are the same: in the combination these contributions drop out. For this reason, the analysis of the asymptotic behaviour simplifies considerably if the pair T 1 , T 2 is replaced by the pairT , T 2 , which is also free of kinematic zeros.

Dispersion relations
The dispersion relations express the Compton amplitude in terms of the structure functions. These represent the Fourier transform of the current commutator: The structure functions are experimentally accessible only for q 2 ≤ 0 and it is customary to replace q 2 by Q 2 ≡ −q 2 . In the standard notation, where the structure functions are denoted by F 1 (x, Q 2 ), F 2 (x, Q 2 ) with x = Q 2 /2mν, V 1 and V 2 are given by: ForT , the structure functionV = V 1 + 1 2 V 2 is relevant: We assume that the Compton amplitude exhibits Regge behaviour for ν → ∞:T ∝ ν α , T 2 ∝ ν α−2 . Accordingly, the dispersion relation forT requires a subtraction while T 2 obeys an unsubtracted dispersion relation: S(q 2 ) represents the subtraction function, ν 2 0 is the subtraction point in the variable ν 2 and the lower limit corresponds to the threshold for inelastic reactions, ν th = (2mM π + M 2 π − q 2 )/(2m). The elastic part of T is given byT el = T el 1 + 1 2 T el 2 . As such, the choice of the subtraction point is arbitrary (provided that ν 2 0 < ν 2 th ), but as pointed out in [1], it is convenient to set ν 2 0 = − 1 4 Q 2 rather than to subtract at ν 0 = 0. As will be seen below, this choice simplifies the asymptotic behaviour of the subtraction function for Q 2 → ∞. Replacing the variable of integration ν ′ by x = Q 2 /(2mν ′ ), the dispersive representation then takes the form ,

Reggeon dominance
While T 2 is fully determined by the form factors and the structure functions because it obeys an unsubtracted dispersion relation, the representation forT involves a subtraction function, which causality alone leaves undetermined. This illustrates a venerable theorem which concerns the implications of causality for the structure functions [21][22][23]. The theorem states that the values ofV (ν, q 2 ), V 2 (ν, q 2 ) in the space-like region q 2 ≤ 0 determine these functions in the time-like region, up to a polynomial in the variable ν. The implications for the dispersive analysis of the Compton amplitude are discussed in [24]. In Regge language, integer powers of ν are called fixed poles: the continuation from the space-like to the time-like region is unique up to fixed poles. Regge asymptotics excludes such contributions in V 2 , but the continuation ofV into the time-like region is unique only up to a term that depends on ν exclusively through the step function: where σ(s) vanishes for s < 0. InT , the ambiguity shows up in the form which is independent of ν and thus only affects the subtraction function. Since the ambiguity amounts to a superposition of free propagators, it is evident that the Fourier transform ofV fp vanishes outside the light-cone. The nontrivial part of the theorem is that the values of the structure functions in the space-like region, where they can be measured, fully determine the amplitude in the time-like region up to contributions of this particular type.
In QED, the electrons Reggeize [25], but the photon remains elementary [26]. QCD, however, does satisfy the criteria for Reggeization formulated in [26]: the gluons as well as the quarks have this property [27,28]. In the meantime, the graphs that need to be summed up to study the high energy properties of the scattering amplitudes within QCD perturbation theory have been identified and Reggeon field theory has been developed for the analysis of exchanges of more than one Reggeon, in particular also of the cuts generated by these [29][30][31][32].
It is generally assumed that the asymptotic behaviour of the Compton amplitude is indeed governed by Reggeon exchange. The exchange of Reggeons generates contributions which at high energies are of the form where s = m 2 + 2mν + q 2 and u = m 2 − 2mν + q 2 represent the square of the centre of mass energy in the s-and u-channels, respectively. In general, the power α depends on t: α(t) moves on a Regge trajectory. In our context, however, only the forward scattering amplitude is of interest, so that only the intercept α = α(0) is relevant. For the Compton amplitude of the proton or the neutron, the Pomeron yields the dominant contribution; it involves a superposition of terms of the above form with intercepts in the vicinity of α = 1. The Reggeon with the quantum numbers of the f 2 and an intercept in the vicinity of α = 1 2 represents the most important non-leading contribution. In the difference between the amplitudes relevant for proton and neutron, the Pomeron drops out -the Reggeon with the quantum numbers of the f 2 then represents the leading term.
We assume that the Reggeons dominate the asymptotic behaviour [5]: This amounts to the assumption that the amplitudeT does not contain a fixed pole at α = 0. We do not know of a physical phenomenon that could produce a fixed pole at α = 0 inT . Neither causality nor the shortdistance singularities nor the Reggeons generate terms of this sort. The constraints imposed on the subtraction function by causality and unitarity have been analyzed within the alternative dispersive framework set up in [33]. Model-independent bounds for the subtraction function S 1 (q 2 ) are derived and it is shown that the results obtained in [12] from Reggeon dominance at low values of Q 2 are consistent with these. An extension of this work to the higher values of Q 2 investigated in the present paper would be most welcome, as it would allow to subject Reggeon dominance to a further test. Modelindependent bounds on the subtraction functionS(q 2 ) would be particularly interesting, because the operator product expansion of this quantity is free of the short distance singularities associated with operators of spin 2 -asymptotic freedom fully determines the asymptotic behaviour ofS.

Sum rule for the subtraction function
As pointed out in [5], Reggeon dominance determines the subtraction function in terms of the cross sections of inelastic scattering. In [12], a sum rule for the subtraction function S 1 (−Q 2 ) ≡ T 1 (0, −Q 2 ) was derived that represents the inelastic part of this quantity in terms of integrals over the cross sections. We now derive an analogous sum rule that expresses the subtraction func-tionS(−Q 2 ) as an integral over the structure function F (x, Q 2 ) -an immediate consequence of the dispersion relation (15) and Reggeon dominance (19). The derivation is not trivial, however, because the Reggeons generate singularities at x = 0. The leading singularity inF is of the form: For this reason, taking the limit ν → ∞ in the dispersion relation (15) requires some care: the limit cannot simply be exchanged with the dispersion integral because, in the vicinity of x = 0, the term 4m 2 x 2 ν 2 fails to dominate over Q 4 . Since the elastic part of the amplitude tends to zero for ν → ∞, the Reggeon dominance hypothesis (19) amounts to the requirement that the subtraction function cancels the limiting value of the difference between the dispersion integral and the asymptotic representation: The limit is worked out in Appendix C. The result takes the form of a sum rule that expresses the subtraction functionS in terms of the structure functionF A Finite Energy Sum Rule variant of this relation was proposed by Elitzur and Harari [4]. The above formulation shows that the sum rule is perfectly consistent with the scaling violations required by QCD -contrary to statements made in [18].

Wick rotation
We now return to the mass formula (3). In the decomposition introduced above, the trace of the Compton amplitude is given by The integration in (3) runs over all q 2 , space-like as well as time-like. For the integral to converge, it needs to be regularized, for instance by replacing the photon propagator 1/q 2 with Λ 2 /(Λ 2 − q 2 )/q 2 .
In the rest frame of the particle, ν coincides with the component q 0 of the photon momentum. Cottingham [3] observed that the time-ordered amplitude is analytic in q 0 and that the path of integration over this variable may be rotated into the imaginary axis, at fixed three momentum q -without crossing any singularities of the integrand in (3) (Wick rotation).
Setting q 0 = iQ 4 and identifying Q 1 , Q 2 , Q 3 with the space components of the physical momentum, m γ takes the form of a euclidean integral extending over the four-vector Q µ : The result for the renormalized mass difference is independent of the form used for the regularization. It is customary to use a cutoff in momentum space: restrict the integration to the euclidean sphere Q 2 ≤ Λ 2 and write the regularized Cottingham formula as In this formula, the amplitudesT and T 2 are to be evaluated at ν = iQ 4 , q 2 = −Q 2 .

Decomposition of the mass shift
In the framework of QCD+QED, the mass of a particle is determined by the bare parameters that occur in the Lagrangian and the cutoff used to regularize the theory. If the electromagnetic interaction is turned off, only the QCD coupling constant, the quark masses and the cutoff are relevant. To order e 2 , the e.m. interaction changes the mass not only by the above integral, but in addition by the contribution ∆m Λ , which arises from the change in the bare parameters needed for the mass of the particle to stay finite when the cutoff Λ is removed -the bare quantities depend on Λ as well as on e. The e.m. part of the mass is obtained by adding this contribution to the integral in equation (26): Inserting the dispersion relations (15) in formula (26), we obtain a representation of the e.m. part of the mass as a sum of four terms: In each one of these, the integrals over the direction of the vector Q µ can be done explicitly.
In the first term, which collects the contributions fromT el and T el 2 , this leads to a set of integrals over the form factors of the nucleon, which are known very accurately -an explicit expression for m el is given, for instance, in [12].
The second and third term arise from the dispersion integrals over the structure functionsF and F 2 , respectively. With the above choice of the subtraction point, the integrands are proportional to the factor Q 2 − 4Q 2 4 , in both cases. Taken by itself, the angular integral of this factor over the directions of the vector Q µ vanishes. Moreover, when Q 2 becomes large, the remainder of the integrand becomes independent of Q 4 . Hence the angular integration suppresses the contributions arising from large values of Q 2 -these integrals approach finite limits when the cutoff is removed: The normalization constant is given by the variable y stands for y = Q 2 /(4m 2 x 2 ) and the explicit expression for the function f (y) reads The suppression of the angular integrals relevant for mF and m F2 manifests itself in the fact that the function f (y) rapidly falls off when y becomes large: The angular integral can be done explicitly in the fourth term as well, but there, it does not suppress the contributions from large values of Q 2 , so that the cutoff must be retained: Together with the sum rule (23), the above formulae fully specify the e.m. part of the mass difference between proton and neutron, in terms of measurable quantities. The next four sections concern asymptotic properties of the Compton amplitude that are not of direct relevance for the Cottingham formula -the contributions generated by short distance singularities of spin 2, for instance. If the reader is more interested in the numerical outcome of our analysis for the mass difference, he or she may go directly to section 12.

Operator product expansion
The behaviour of the amplitudesT and T 2 at large momenta is determined by the short-distance properties of the matrix element p|T j µ (x)j ν (0)|p , which can be analyzed by means of the operator product expansion (OPE) [34]. The asymptotic freedom of QCD implies that perturbation theory can be used to work out the leading terms of this expansion [15][16][17][18]. For the timeordered product of two currents, the behaviour at short distances z = x − y is of the form where O n enumerates the renormalized gauge invariant local operators of QCD and X = 1 2 (x + y). The expansion starts with the operators of lowest dimension. The Wilson coefficientsC n µν (z) vanish unless O n has the same flavour quantum numbers as the product of two e.m. currents. The symmetry of QCD under P, T and C also prevents some operators from contributing to the expansion. The coefficients depend in a nontrivial manner on z only through z 2 : they are polynomials in the components of the vector z, with coefficients that depend on z 2 .
In momentum space, the OPE governs the behaviour at large momenta. We denote the Fourier transform with respect to z byT µν : The limit z = λz, λ → 0 in coordinate space corresponds to the limit q = λq, λ → ∞ in momentum space. We refer to this limit as q → ∞. In this notation, we havẽ The coefficients C n µν (q) represent the Fourier transforms of those in equation (35) and are polynomials in the components of q, with coefficients that depend on q 2 . This immediately implies that the coefficients occurring in the expansion of the invariant amplitudes T 1 (ν, q 2 ), T 2 (ν, q 2 ) are polynomials in ν. The expansion thus also holds for imaginary values of ν.
A contribution from the unit operator only occurs in the disconnected part and does not show up in the scattering amplitude. In QCD, the relevant operators of lowest dimension aref f with f = u, d, . . . They have spin zero and are of engineering dimension 3. Chiral symmetry, however, suppresses the contributions from these operators: their Wilson coefficients are proportional to the masses m u , m d , . . . It is convenient to include the mass factor and to work with the operator O f0 = m ff f , which is of dimension 4. Lorentz invariance implies that the spin averaged matrix elements of the operators of lowest dimension can be expressed in terms of the following linearly independent operators, which either have spin 0 or spin 2: Accordingly, the leading terms in the operator product expansion ofT µν (q, X) are given bỹ While the dependence on q resides in the Wilson coefficients, the operators only depend on X.

Leading Wilson coefficients
Lorentz invariance and current conservation imply that the Wilson coefficient of the scalar operator O f0 is proportional to the kinematic tensor K 1 µν specified in equation (6): the contribution from this operator is of the form where the coefficient c f 1 depends on q 2 . The contribution from O G0 is of the same structure.
For operators with spin 2, the situation is not that simple. As shown in Appendix A, Lorentz invariance and current conservation determine the form of their Wilson coefficients only up to two functions of q 2 , which we denote by c 2 (q 2 ) and c 3 (q 2 ). According to equation (A.3), the contribution from O f2 αβ is of the form: and analogously for the contribution from the lowest dimensional gluonic operator of spin 2. This shows that kinematics determines the Wilson coefficients of the lowest dimensional operators in terms of six functions, c f 1 , . . . , c G 3 , that depend on q 2 . The amplitude we are interested in represents the spin average of the one-particle matrix element (Since the initial and final momenta are the same, the matrix element is independent of X). Inserting the expansion (39), we obtain: For scalar operators, the matrix element is a constant, while for spin 2, it depends on the momentum of the particle: According to equation (40), the contributions from the spin zero operators are proportional to the kinematic tensor K 1 µν . For the spin 2 terms, a short calculation is needed to verify that it can be expressed as a linear combination of K 1 µν and K 2 µν : The terms arising from the lowest dimensional gluonic operator of spin 2 are of the same form.
The corresponding asymptotic representations for T 1 (ν, q 2 ) and T 2 (ν, q 2 ) are given by the coefficients of K 1 µν and K 2 µν , respectively: While the leading term in the asymptotic behaviour of T 2 (ν, q 2 ) only depends on q 2 , T 1 (ν, q 2 ) contains a term proportional to ν 2 . The advantage of working with the amplitudeT now becomes visible: T 1 and T 2 contain a common spin 2 contribution. In the combinationT = T 1 + 1 2 T 2 , this term drops out -only the one proportional to the factor ν 2 − 1 4 q 2 remains: As noted above, the angular integration suppresses contributions that are proportional to this factor -this is the reason why in the decomposition (28) only mS contains a divergence.

Difference between proton and neutron
For the mass difference between proton and neutron, only the difference between the Compton amplitudes of the two particles is needed. As far as the asymptotic behaviour is concerned, we thus only need the difference between the spin averaged matrix elements of proton and neutron. In the isospin limit, the neutron matrix elements of the gluonic operators O G0 and O G2 αβ coincide with those of the proton. In reality, since m u differs from m d , the proton and neutron matrix elements of the gluonic operators are slightly different, but in the mass difference between proton and neutron, this generates an effect of second order in isospin breaking and will be neglected. This simplifies matters considerably. Only the matrix elements of non-singlet operators are relevant -operator mixing does not affect these.
Throughout the remainder of this paper, we focus on the difference between the Compton amplitudes of proton and neutron, without explicitly indicating this in the notation: in the following,T and T 2 stand for T p−n and T p−n 2 , respectively.

Perturbation theory
To leading order of the QCD perturbation series, the Wilson coefficients are the same as for free quarks. The explicit expressions are readily obtained by simply replacing the nucleon in the above relations with a free quark of charge Q f = 2 3 or − 1 3 . If the strong interaction is turned off and the e.m. interaction is accounted for only to leading order, the Compton scattering on a quark is elastic. The Sachs form factors are given by In the limit q → ∞ relevant for the OPE, the second term in the denominator becomes negligible compared to the first: for space-like momenta, T f 2 tends to The spin averaged quark matrix elements of the operators O f0 and Q f2 αβ are readily worked out; they yield O f0 = 2m 2 f and O f2 = 4m 2 f . The leading terms in the expansion of the coefficients c f 1 , c f 2 , c f 3 in powers of g 2 can then be read off from the asymptotic relations (46) and (47) [18,35]: In this calculation, the spin of the operators occurring in the OPE does not play an important role. Appendix B contains an alternative derivation of these relations, which is based on the short distance expansion of the quark propagator and explicitly exhibits the spin structure.
The higher order terms in the expansion of the Wilson coefficients have been studied in detail, also for the gluonic operators [15][16][17][18]35] -for a thorough review, we refer to [36]. The qualitative features of the asymptotic structure are intimately related to the fact that the dimension of the spin 2 operators is anomalous. The correction of order g 2 in the perturbative expansion of the quantity Q 4 c f 2 (−Q 2 ) falls logarithmically if Q 2 becomes large. With the renormalization group, the leading logarithms can be summed up to all orders. The contributions from the singlet operators undergo mixing, but as noted above, for the difference between proton and neutron, only the nonsinglet operators are relevant. The matrix element of the term involving c f where Λ QCD is the renormalization group invariant scale of QCD, while d 2 is related to the anomalous dimension of the operator O f2 αβ and depends on the number of flavours: The formula (51) holds provided Q is large, not only compared to Λ QCD , but compared to all of the quark masses. In the intermediate range where Q is large compared to m s , but not large enough to activate the degrees of freedom of the heavy quarks, it should hold approximately, with N f ≈ 3.
Since the perturbation series of c f 3 (−Q 2 ) only starts at order g 2 , the asymptotic behaviour is suppressed by a factor of ln The scalar operatorf f is of anomalous dimension as well, but the same is true of m f and the anomalies cancel: the operator m ff f is renormalization group invariant. This implies that, in the Wilson coefficient c f 1 (−Q 2 ), the correction of order g 2 does not pick up a logarithmic enhancement if Q 2 becomes large and there is nothing to be summed up: Note that the above relations only account for the leading logarithms. The perturbation series of the coefficient c f 1 (−Q 2 ) does contain contributions of order g 2 that are not enhanced by a logarithm -their role in the context of the Cottingham formula will be discussed in section 15.
Inserting the asymptotic expressions for the Wilson coefficients in equations (47) and (48), we obtain: While the coefficient C is determined by the spin averaged matrix elements of a renormalization group invariant operator, C 2 and C 3 are related to the matrix elements of the spin 2 operator O f2 αβ , which do depend on the renormalization convention used.

Moments of the structure functions
Let us now compare the dispersive representation (15) with the asymptotic formulae (55) obtained from perturbation theory. Since the form factors rapidly tend to zero when Q 2 becomes large, the elastic part of the amplitudes does not show up in the asymptotic behaviour. The dispersion integrals approach moments of the structure functions: In our decomposition of the dispersive representation, the asymptotic behaviour (55) boils down to a set of conditions on the subtraction function and on the lowest moments ofF , F 2 and F L : In the literature, the perturbative predictions for the moments have been compared in detail with experiment [36]. The parametrization we will be using for the structure functions is based on the DGLAP equations [37][38][39]. These ensure that the behaviour of the Compton amplitude in the deep inelastic region is consistent with perturbation theory.

Prediction for the constant C
Neglecting isospin breaking effects of second order, the neutron matrix elements of e 2d d agree with the proton matrix elements of e 2ū u and vice versa. The constant C can thus be expressed in terms of proton matrix elements: The matrix element of the operatorūu −dd also determines the leading contribution to the QCD part of the proton-neutron mass difference (see e.g. [40]): This shows that the constant C is related to the value of the proton-neutron mass difference in the absence of the e.m. interaction: Once we have determined the e.m. part, the observed mass difference will provide us with a value of m QCD and hence also with a value of the constant C.
Actually, however, the precise value of C is not crucial in our context. For our purpose, the crude estimate m QCD ≈ −2 MeV is good enough. The quark mass is not yet known very firmly. The FLAG result m u /m d = 0.513 (39) (for N f = 2 + 1 + 1) [41] implies r = 2.16 (42). There is a totally independent determination of the mass ratio Q, based a low energy theorem for the decay η → 3π [42][43][44]. A recent analysis of the data on this basis leads to Q = 22.1(7) [45]. Combining this result with the well-determined ratio m s /m ud = 27.23(10) [41], we obtain m u /m d = 0.450(25) and r = 1.46 (25). As pointed out in [45], the origin of the difference could be identified by calculating the corrections to the low energy theorem on the lattice, but this yet needs to be done. The outcome for the constant C is tiny. With the value r = 1.46, we obtain: The reason why the value turns out to be so small is that C vanishes in the chiral limit. It implies that C is small compared to C 2 and C 3 . Accordingly, it takes very large values of Q 2 for the singularities generated by the operators of spin 0 to finally dominate over those associated with operators of spin 2.

Renormalization
The asymptotic behaviour of the subtraction function in equation (58) implies that the corresponding contribution to m QED is logarithmically divergent. The leading divergence is determined by the coefficient C: A logarithm also occurs in the electromagnetic renormalization of the bare QCD coupling constant and of the bare quark masses (see e.g. [40]): The scale µ of the logarithm is a matter of conventionpicking a value for µ amounts to fixing the ambiguity in the decomposition (1) of the mass difference into contributions arising from the e.m. and strong interactions, respectively.
In the difference between the masses of proton and neutron, the e.m. renormalization of the coupling constant only yields a contribution of second order in isospin breaking -we are neglecting such effects. The renormalization of the quark masses, on the other hand, does not drop out in the difference. In the Lagrangian, the corresponding counter term reads The corresponding shift in the mass of a particle is given by − p|∆L|p /2m. Accordingly, the change in the proton mass generated by the renormalization of m u and m d is given by Again neglecting effects of second order in isospin breaking, the neutron matrix elements can be expressed in terms of those of the proton: Collecting terms and neglecting second order isospin breaking effects, we obtain the following expression for the counter term ∆m Λ = ∆m p − ∆m n : Comparison with the expression (62) for the constant C that determines the asymptotic behaviour of the subtraction function shows that the two quantities are related by where the normalization factor N is specified in equation (31). As it should be, the logarithm in the integral (66) over the subtraction function thus cancels the one in the renormalization (68) of the quark masses: the leading divergences occurring in the expression (34) for mS drop out.

Subleading divergence
As mentioned above, the asymptotic formulae (55) only account for the leading logarithms -they are valid only up to corrections of order g 2 . This applies, in particular also to the Wilson coefficient of the spin 0 operator that is responsible for the logarithmic divergence of the Cottingham formula. The correction of order g 2 gives rise to a theoretical issue, which does not appear to be covered in the literature and which we briefly wish to address. When Q 2 becomes large, the effective strength of the interaction decreases in proportion to 1/ ln Q 2 . Those corrections of order g 2 in the Wilson coefficients or in the counter term ∆m Λ that do not pick up a logarithmic enhancement are asymptotically small compared to the leading terms. This does not ensure, however, that the corresponding contributions to mS remain finite when the cutoff is removed: the corresponding contributions instead grow in proportion to The coefficient of order e 2 g 2 in the renormalization of the quark masses in QCD+QED is known [46,47]: where γ 0 = 2 and γ 1 = 101 12 − 5 18 N f are the well-known coefficients relevant for mass renormalization in QCD. The counter term ∆m Λ considered above is related to the term of order e 2 in equation (76).
On the other hand, the contributions of order g 2 in the Wilson coefficient were considered by Shifman, Vainshtein and Zakharov, more than 40 years ago [48]. Equations (4.15) and (4.18) in this reference indicate that, in the notation used above, the coefficient C picks up the same correction as the counter term, As this amounts to combining results obtained within two different regularization schemes (cutoff in euclidean momentum space, dimensional regularization) it must be taken with a grain of salt, but it does indicate that the divergences of the type ln ln Λ 2 cancel. Irrespective of the regularization used, the renormalization of coupling constants and quark masses must remove the divergences also at the subleading level. Numerically, the perturbative corrections are not important, because, as pointed out above, chiral symmetry suppresses the entire contribution to the mass difference from the region where perturbation theory applies. In that region, the corrections are even smaller than the leading terms -they are in the noise of our calculation and we neglect them. The limit Λ → ∞ in formula (34) can then be done explicitly. The result can be written in the form withμ = µ exp(− 1 2 ).

Input used for the structure functions
For the numerical evaluation of the inelastic contributions, we need a representation for the difference between the structure functions of proton and neutron, and not only for the relatively well explored quantity F 2 , but also for the longitudinal component F L , which is known less well. At low values of Q 2 , we closely follow the analysis of [12] and distinguish three different regions in the centre of mass energy W = √ s (numerical values for W and Q 2 are given in GeV units): (i) For the range W < 1.3, we rely on the parametrizations of the structure functions of MAID and DMT [49][50][51] -we refer to these as MD. Both of them are accessible on the MAID home page [52]. We identify the central values of the structure functions in this region with the mean of the two parametrizations and use the difference as an error estimate (half of the difference would suffice to cover the two). The green error band in Fig. 1, shows the corresponding representation of the structure func-tionF for Q 2 = 1. (ii) In the interval 1.3 < W < 3, we make use of the representation due to Bosted and Christy (BC) [53,54]. It contains a wealth of information, but suffers from several shortcomings that are discussed in detail in section 5.1 of [12]. Part of the problem originates in the fact that the longitudinal cross section is more difficult to measure than the transverse one.
In [53,54] it is assumed that the ratio R = σ L /σ T of the neutron cross sections is the same as for the proton. In the region where the Pomeron dominates, this holds to good accuracy, but we need the difference between the two, where Pomeron exchange drops out. The assumption amounts to using an approximation and introduces a systematic error that is not easy to estimate. In our opinion, the procedure used in [12] to cope with the uncertainties in the region 1.3 < W < 3 is on the conservative side and we adopt it in the present work: we treat the transverse and longitudinal cross sections as independent and assign an uncertainty in σ p−n T and σ p−n L of 8% of σ p T and 8% of σ p L , respectively. In part of phase space, this may well overestimate the uncertainties considerably -a reanalysis of the data in the resonance region would be most welcome. The structure of the brown error band reflects the resonances occurring in this region. (iii) For W > 3, Q 2 < 1, we rely on the parametrization of the proton structure functions due to Alwall and Ingelman (AI) [55]. It represents the amplitude as a sum of a contribution from the Pomeron and one from the Reggeon with the quantum numbers of the f 2 . In the difference between the proton and neu-tron amplitudes, the Pomeron drops out. We assume that the couplings of the f 2 -Reggeon to proton and neutron are approximately SU(3)-symmetric and attach an uncertainty of 30% to the representation for the difference between proton and neutron obtained on this basis. The blue band in Fig. 1 shows the corresponding uncertainty range at Q 2 = 1. For details, we refer to section 5.1 in [12]. (iv) In the region W > 3, Q 2 > 1, we use the solution of the DGLAP equations constructed by Alekhin, Blümlein and Moch, who obained numerical values for the structure functions over a wide range: 1 ≤ Q 2 ≤ 2 · 10 5 and 10 −7 ≤ x ≤ 0.99. The values of F 2 (x, Q 2 ) and F L (x, Q 2 ) are listed for the proton as well as for the neutron on a grid of 60 × 98 points. We thank Johannes Blümlein for providing us with this table, which we refer to with the acronym ABM. The underlying analysis is described in [56][57][58].
In the deep inelastic region, asymptotic freedom leads to very strong constraints, particularly for the structure function F L . The strength of these constraints is clearly visible at leading order of the perturbative expansion, where F L is given by an integral over F 2 . The DGLAP equations extend this relationship to higher orders, by means of the renormalization group. In our framework, the properties of F L play a crucial role in the evaluation of the sum rule for the subtraction function. The theoretical constraints on this quantity are very important for our analysis, particularly also because the raw experimental information for F L is much weaker than the one for F 2 .
The black dots in Fig. 1 show the numbers obtained for W andF from the entries for x, F 2 and F L , at the lowest value of Q 2 listed in the ABM table, Q 2 = 1, and W > 3. While the result agrees very well with AI for W > 5, the two representations do differ at lower values of W . Since the DGLAP equations rely on perturbation theory, we should not be surprised to find deviations at low momenta, i.e. in the region where Q 2 as well as W are small.
where κ is the anomalous magnetic moment of the particle (these relations hold separately for proton and neutron). The dispersion relation for T 2 converts the second one into the Baldin sum rule [59], which expresses the sum α E + β M as an integral over the cross section for photoproduction. For the combinationT = T 1 + 1 2 T 2 we are working with, the low energy theorem involves the difference between the electric and magnetic polarizabilities: It fixes the value of the subtraction functionS(q 2 ) at q 2 = 0 in terms of the polarizabilities: For q 2 = 0, our sum rule for the subtraction function thus represents an analog of the Baldin sum rule: it determines the value of the difference between the electric and magnetic polarizabilities rather than their sum, in terms of the structure functions. While the Baldin sum rule directly follows from the unsubtracted dispersion relation for T 2 , the one for α E − β M relies on Reggeon dominance. The integrals over the structure functions relevant for the evaluation of the subtraction function in the dispersion relation for T 1 at small values of Q 2 were analyzed in detail in section 5 of [12]. As shown there, the prediction for the electric polarizability comes with comparatively small uncertainties: 4 α p−n E = −1.7(4) [12] .
The averages for proton and neutron quoted by the Particle Data Group yield α p−n E = −0.6(1.2) [60]. The fact that experiment agrees with the prediction within errors provides a test of Reggeon dominance.
For a review of the currently available information about the polarizabilities, in particular also of the analysis based on chiral effective theories, we refer to [61,62]. In the framework of χPT, the representation of the virtual Compton scattering amplitude has been worked out to first nonleading order [63]. In this reference, the low energy singularity generated by the ∆(1232) resonance is explicitly accounted for. It will be of considerable interest to compare the result of this analysis for the slope of the subtraction function at Q 2 = 0 with the solution of the sum rule that follows from Reggeon dominance constructed in the present paper.
The Baldin sum rule and the data on photoproduction imply that the sum α E + β M is known more accurately than the individual terms. For this reason, it is 4 As usual, the polarizabilities are given in units of 10 −4 fm 3 .
useful to treat α E ± β M as the two independent quantities. The results quoted for proton and neutron in the compilation of Melendez et al. lead to Combining the prediction (82), which is based on Reggeon dominance, with the result (83) obtained from photoproduction, we obtain a prediction for the magnetic polarizability, which is slightly more accurate than the one given in [12]: There were early attempts at calculating the electric polarizabilities on the lattice [65][66][67], based on turning on a constant external electric field, but they did not reach a level where the results could be compared with the experimental determinations in a meaningful way. The very recent lattice determination of the magnetic polarizabilities, however, which makes use of a constant external magnetic field, does yield a remarkably precise value for β p−n M , in good agreement with our predicton in equation (85).
In connection with the proton-neutron mass difference, the polarizabilities are of interest because they determine the value of the subtraction functionS(q 2 ) at q 2 = 0, according to (81). The prediction (82) for α p−n E and the experimental value (83) of (α E + β M ) p−n implȳ The uncertainty is significantly smaller than the one obtained from the experimental value of (α E − β M ) p−n : On the other hand, combining the lattice result (86) for the magnetic polarizability with the experimental value (83) of (α E + β M ) p−n , we obtain a result for the subtraction function at the origin that is even slightly more precise than the prediction: The fact that, within errors, this result agrees with the prediction (87) amounts to a more stringent test of the Reggeon dominance hypothesis than the one discussed above. It is important to pursue the determination of the polarizabilities; in particular, the pioneering lattice result which provides such a test calls for confirmation.
Next, we discuss the solution of the sum rule (23) for Q 2 < 1, where the parametrizations listed in (i) -(iii) suffice. Fig. 2 displays the contributions arising from the various regions of phase space. The interval of integration in (23) is split into three parts that correspond to the regions where we are using the representations MD, BC and AI, respectively. The values of x where W = 1.3 and W = 3 are denoted by x a and x b , respectively. In the first two parts, the integration over the termF R /x 2 can explicitly be done -we book these contributions together with the term involving the Reggeon residues inS AI . For Q 2 < 1, the solution of the sum rule then takes the form The termS MD (−Q 2 ) includes the most prominent low energy phenomenon, the resonance ∆(1232). Isospin conservation ensures that the couplings of this state to proton and neutron are the same, so that the resonance does not show up at all in the subtraction function relevant for the difference between the two. Indeed, as shown by the green band, the contributions from this region are small.
In the region of the higher resonances, we rely on the BC representation of the structure functions. The brown error band indicates the price to pay with the error estimate specified in section 16: the largest uncertainty in our evaluation of the mass difference stems from there.
The blue band depicts the function S AI . Since the Regge representation AI we are using in this region is restricted to Q 2 < 1, the band stops at Q 2 = 1.
The plot shows that the contributions from MD and AI are negative, while the one from BC is predominantly positive. The net central value ofS(−Q 2 ), which is indicated by the black line, is rather small and negative, but the uncertainty attached to it (gray error band) excludes positive values only in the vicinity of Q 2 = 0.
The cyan-coloured wedge labeled A represents the tangent Q 2S (−Q 2 ) = Q 2S (0) + . . . calculated with the value ofS(0) in equation (89)  netic polarizability plus Baldin sum rule). It confirms that at small values of Q 2 , the subtraction function is negative.
19 Intermediate values of Q 2 The representations of MD and BC are valid also for Q 2 > 1, but for AI, this is not the case. We instead rely on ABM. The formula for the corresponding contribution to the subtraction function is the same as for AI, but the representation forF consists of a numerical table rather than an algebraic parametrization like the one of AI. The leading term in equation (20) stems from the Reggeon with the quantum numbers of the f 2 and α ≃ 0.55. In order to determine the corresponding coefficient b α , we focus on small values of x and approximate the numbers forF obtained from the ABM table at a given value of Q 2 with an approximation of the form The coefficients b α , b ′ α , b ′′ α depend on Q 2 . We determine them by minimizing the sum of the squares of the differences between the parametrization and the ABM values over a suitable interval. At very small values of x, the numerical noise in the entries of the table hides the signal while if x is too large, the approximation used breaks down -we find that 10 −4 < x < x 1 with x 1 = 3 · 10 −2 represents a suitable range. In the grid of x-values used in the ABM table, this range contains points # 15 to 25. We fix the parameter b ′′ α with continuity at point # 24 and treat the coefficients b α , b ′ α as free parameters. For a given value of Q 2 , the minimization then fixes these. In particular, the procedure determines the Reggeon residue, which according to (20) is given by β α = 1 2 Q −2(α+1) b α . Fig. 3 compares the fit (red curves) with the values ofF obtained from the ABM data (black dots), for various values of Q 2 , including the lowest and highest ones listed in the table. In the region where the Pomeron dominates, the values ofF p andF n are nearly the same. It is difficult to reliably determine the difference between the two from the data on inelastic scattering, even if the DGLAP equations provide a strong theoretical constraint. In the ABM table, the problem also manifests itself directly: for Q 2 > 3.5, the results for b α exhibit fluctuations which are generated by the limited numerical accuracy of the entries and are visible in Fig. 3. On the other hand, it is questionable, whether the ABM data can be trusted down to Q 2 = 1, because the DGLAP equations rely on perturbation theory. For these reasons we assign an overall relative error of 30% to the numbers for the difference between the structure functions of proton and neutron obtained from ABM. Fig. 4 compares the Reggeon residue extracted from the ABM table with the values for this quantity obtained from the parametrization of AI. For better visibility, the value of β α is plotted on a logarithmic scale. The figure shows that, at Q 2 = 1, where the two representations meet, the results agree within errors: the two entirely different sources match, both in sign and in size.
Concerning the evaluation of the sum rule for the subtraction function, the only difference compared to the preceding section is that the AI representations for F and b α are replaced by those obtained on the basis of ABM. The black dots in Fig. 5 show the outcome

Vector Meson Dominance
As discussed in detail above, the asymptotic freedom of QCD implies that the subtraction function obeys the asymptotic condition (58):S → C/Q 4 when Q 2 becomes large. The constant C does not represent an unknown, but can be expressed in terms of the mass difference in QCD. Since C is suppressed by chiral symmetry, it is tiny: C ≈ 6 · 10 −4 GeV 2 .
In the subtraction function, the numerical noise mentioned above starts becoming visible at Q 2 ≈ 3.5 and, for Q 2 > 6, it hides the signal completely: there, S vanishes within errors.
In order to interpolate between the values of Q 2 where the ABM table provides significant information and the region where asymptotics sets in, we make use of the Generalized Vector Dominance Model of Sakurai and Schildknecht [69], parametrizing the subtraction function in terms of the contributions from ρ, ω and φ. In the difference between proton and neutron, only the off-diagonal terms survive: The asymptotic condition requires the two terms in the bracket to nearly cancel: This leaves a single parameter free, say c ω . We determine this parameter by fitting the model to the values obtained from MD + BC + ABM in the region 2 < Q 2 < 3.5. This range excludes values of Q 2 below 2, where the validity of the DGLAP equations is questionable as well as the region Q 2 > 3.5, where the fluctuations show up. The minimum occurs at The red band in Fig. 5 shows this fit.
Since the Q 2 -dependence of the VMD parametrization reproduces our results very well, the outcome for mS is not sensitive to the range used in the fit -as long as it does not extend into the region Q 2 > 6, where the numerical fluctuations take over. The dashed red lines indicate the behaviour of the VMD parametrization at low values of Q 2 . Remarkably, although only input for Q 2 > 2 was used, it shows a reasonable behaviour also at low energies. In fact, the central VMD parametrization runs within the error band obtained from the experimental information in the region Q 2 < 1. Evaluating the representation (92) at Q 2 = 0, for instance, and using the relation (81) betweenS(0) and the polarizabilities, we obtain α p−n E − β p−n M = −1.1 (7). This is about four times more accurate than the available experimental information (84) and perfectly consistent with it.
We emphasize, however, that the particular form of the parametrization used to interpolate between low and high values of Q 2 does not play a significant role. A parametrization of the form proposed by Erben et al. [8], is adequate as well, because it does have the proper asymptotic behaviour. Fixing m 0 at the central value used in that reference, treating c 0 as a free parameter and fitting it to the values ofS obtained from MD + BC + ABM in the region 2 < Q 2 < 3.5, the result for the subtraction function can barely be distinguished from the one obtained with the VMD parametrization. observation at µ = M Z . Fitting the numerical results for the moments in the range between Q 2 = 5 · 10 3 and the upper end of the table provided by ABM, we obtain Fig. 6 shows that the asymptotic formulae indeed yield a good approximation all the way down to Q 2 ≈ 100. This property is built in: the ABM analysis is based on the DGLAP equations which in turn rely on perturbation theory. In the region where the effective coupling constant becomes small, the leading terms must dominate. The figure also confirms that M L disappears more rapidly than M 2 by one power of the logarithm, but both moments only fall off very, very slowly. Fig. 7 shows the behaviour of the structure function S at large values of Q 2 , on a logarithmic scale. The red line represents the VMD parametrization (92) of our central result. To make the asymptotic behaviour visible, the vertical axis is stretched with the factor Q 4 . The quantity Q 4S approaches the Wilson coefficient C, which is determined by the proton matrix elements of the spin 0 operator 1 9 (4m u − m d )(ūu −dd) and is indicated by the dashed red line. As discussed in section 15, C picks up a correction of O(g 2 ). The red dots represent the values of the function

Asymptotics
The correction is too small to make a visible difference (at the mass of the Z-boson, which is marked with a star, it increases the value of C by about 1%).
Traditionally, the subtraction function is identified with a multiple of S 1 (−Q 2 ) ≡ T 1 (0, −Q 2 ). The relation between this object and the subtraction function we are working with is readily established by comparing the dispersion relations obeyed byT and T 1 . The quantity to compareS with is the inelastic part of S 1 , The comparison of the two dispersion relations yields In Fig. 7, our result for Q 4 S inel 1 (obtained by subtracting the term Q 4 ∆S from the result for Q 4S ) is shown as a blue line. For large values of Q 2 , the integral over With the asymptotic formulae for the moments, the asymptotic behaviour of S inel 1 thus takes the form This shows that, while the asymptotics ofS is governed by the matrix elements of a scalar operator, S inel 1 picks up additional contributions proportional to the Wilson coefficients C 2 and C 3 , which represent matrix elements of a spin 2 operator.
The qualitative difference in the asymptotic behaviour of Q 4S and Q 4 S inel 1 originates in the fact that (i) the approximate chiral symmetry of QCD suppresses the coefficient C, while C 2 , C 3 are not suppressed -they are larger than C by two to three orders of magnitude; (ii) while the contribution proportional to C is independent of Q 2 , those from C 2 and C 3 fall off logarithmically.
Although, eventually, C dominates Q 4 S inel 1 as well, asymptopia is reached only if Q 2 is so large that the logarithmic suppression of the spin 2 contributions wins over the chiral suppression of those with spin 0 -from Q 2 = 10 2 to Q 2 = 10 3 , the value of Q 4 S inel 1 only shrinks by about 10%.
The counter term ∆m Λ only removes the leading and subleading divergences associated with C. The additional divergence proportional to C 2 does not have anything to do with renormalization and is of purely technical nature: equation (60) shows that the same divergence also shows up in the asymptotic behaviour of T 2 . In the sum of the contributions from S inel 1 and T 2 , the spin 2 divergences cancel [18]. It is difficult, however, to specify the contribution from S inel 1 by itself: the asymptotic formula (100) shows that this contribution diverges unless the non-leading term proportional to C 3 is removed as well as the leading one. Our framework avoids these problems.
22 Numerical evaluation of the mass difference

Form factors, m el
The elastic contribution to the e.m. part of the mass difference is determined by the form factors. In early work, the experimental information about these was adequately described by the dipole formulae (see e.g. appendix A of [12]). They yield 0.63 MeV for the proton and −0.13 MeV for the neutron, so that the elastic contribution to the self-energy difference amounts to m el = 0.76 MeV [5]. In the meantime, the precision to which the form factors are known has increased significantly [71][72][73][74]. Using this information, we obtain The error bar covers the results obtained with the three parametrizations of [71][72][73]. This indicates that, in the difference between the e.m. self-energies of proton and neutron, the departures from the dipole formulae only generate a change of the order of a percent. The uncertainties in the result for the mass difference generated by the elastic part are totally neglible compared to those from the inelastic contributions.

Contribution from the subtraction function
The contribution to mS depends on the scale µ used in the e.m. renormalization of the quark masses. For definiteness, we use µ = µ 2 ≡ 2 GeV. If µ is taken differently, the mass difference changes by 2N C ln(µ/µ 2 ). In the region 0 < Q 2 < 1 our representation of the subtraction function is based on the parametrizations MD, BC and AI (gray band in Fig 5). Inserting this representation in formula (78) we obtain The central value is negative and reduces the elastic contribution by about 5%. The error is twice as large, however, so that small positive contributions from this region are not excluded.
Since the integrand of mS is proportional to Q 2S , small values of Q 2 are suppressed; the fictitious spike occurring there in the parametrization of BC (see Figs. 3-5 in [12]) does not affect the result very strongly, but an improved analysis of the structure functions in the resonance region above the ∆(1232) would allow reducing the quoted uncertainty.
In the region 1 < Q 2 < ∞, we use the VMD parametrization ofS and obtain To account for the correlations between the contributions from the various regions, we determine the net error in mS by evaluating the integral in equation (34) for the upper and lower edges of the error band. This leads to

Contributions from the dispersion integrals
Finally, we evaluate the convergent integrals mF , m F2 in equations (29) and (30). In these integrals, the small x region does not require special care. As mentioned above, the angular integration suppresses the contributions from the deep inelastic region. In fact, a very strong suppression also occurs at low values of Q 2 . Numerically, these integrals are tiny: m F2 = −0.0039(10) MeV .

Result for m QED and m QCD
Collecting the various contributions, the part of the proton-neutron mass difference that is due to the e.m. interaction becomes The observed mass difference then yields The result for m QCD provides a more precise estimate for the leading Wilson coefficient: We have repeated the entire calculation with this input instead of the crude estimate used for this constant. At the quoted accuracy, the results stay put.

Comparison with Lattice calculations
Within QCD, the lattice approach allows a determination of the mass spectrum with steadily increasing precision, not only for the mesons but also for the more difficult case of the baryons. The inclusion of the e.m. interaction gives rise to a serious problem, however, because this interaction is of long range -enclosing the system in a box distorts the results through finite size effects that need to carefully be sorted out. In comparison with the extensive documentation available for lattice determinations of the quark masses within QCD, the literature containing numerical results for m QED is rather scarce. Fig. 8 collects the results we found. Visibly, the likelihood for the results listed to represent statistically independent measurements of the same physical quantity is quite small. Indeed, not all of the errors shown include an estimate for the systematic uncertainties. Also, not all of the listed papers have appeared in print. Some of the results are obtained from a calculation that simulates QCD+QED, others stay within QCD, calculate the part due to the difference between m u and m d and determine the part that comes from the e.m. interaction by comparing the calculated part with the experimental value. It is well-known that the splitting into two parts depends on the convention used, but this is a theoretical problem that does not require numerical simulations. Our numerical result for m QED is dominated by the elastic contribution; the remainder is significantly smaller and negative. The most recent lattice results listed in Fig. 8 are instead larger than the elastic contribution: the remainder is positive and comparable to the elastic term. Clearly, our result is not consistent with that.

Comparison with other evaluations of the Cottingham formula
We are aware of four recent estimates for the protonneutron mass difference based on evaluations of the Cottingham formula: Walker-Loud, Carlson and Miller (WCM) [7,10], Erben, Shanahan, Thomas and Young (ESTY) [8], Thomas, Wang and Young (TWY) [9] and Tomalak [11]. The first three propose models for the subtraction function S 1 (−Q 2 ), using the experimental information concerning the difference between the magnetic polarizabilities of proton and neutron to determine the value of S 1 (0) and making a simple algebraic ansatz for the momentum dependence. A detailed comparison of the models proposed by WCM and ESTY with the results obtained from Reggeon dominance at low values of Q 2 can be found in [12].  [75][76][77][78][79][80][81][82][83], the lower part contains results obtained with the Cottingham formula [5, 7-9, 11, 40], including the outcome of our analysis.
Tomalak [11] also uses the available experimental information about the magnetic polarizabilities, but instead of making an ansatz for the momentum dependence of the subtraction function, he calculates it on the basis of the assumption that -once the contributions from the Reggeons are removed -the amplitudê T 1 = q 2 T 1 +ν 2 T 2 obeys an unsubtracted dispersion relation [84]. Although this assumption resembles Reggeon dominance, we consider it very unlikely that it is correct. For q 2 = 0, for instance, the amplitudeT 1 reduces to ν 2 T 2 . The asymptotic behaviour of this quantity was investigated by Damashek and Gilman [85] and, independently, by Dominguez, Ferro Fontan and Suaya [86]. Their work indicates that f = ν 2 (T 2 − T R 2 ) tends to a nonzero constant when ν becomes large. The assumption used in [11] instead implies that f tends to zero. At any rate, this hypothesis implies a constraint on the imaginary part of T 2 at q 2 = 0, i.e. on the cross section of photoproduction: it leads to a sum rule that requires an integral over the cross section to cancel the Thomson term. We do not know of an argument that would support this assumption.
Incidentally, the assumption used in [11] corresponds to a special case of the universality hypothesis of Brodsky, Llanes-Estrada and Szczepaniak [87][88][89], who do not impose the condition that the differenceT 1 −T R 1 tends to zero for ν → ∞, but postulate that it becomes independent of q 2 . We cannot see any reason for this to be the case in QCD (see also [90][91][92]).

Contributions from the elastic part
Since T 2 obeys an unsubtracted dispersion relation, the corresponding Born term is readily obtained by saturating the dispersion integral with the contributions from the nucleon poles. For T 1 , however, the Born term is not unique -various different expressions are used in the literature. They all obey a subtracted dispersion relation, but differ in the choice of the subtraction function.
Dispersion theory offers a unique solution: since analytic functions are determined by their singularities and their behaviour at infinity, it suffices to impose the condition that the Born term vanishes for ν → ∞. We refer to the resulting expression as the elastic part of the amplitude. It is explicitly given in formula (8) (the unsubtracted dispersion relation used to specify the Born term for T 2 automatically ensures that it disappears if ν becomes large). Accordingly, the elastic part of m QED , which we denote by m el , is an unambiguous notion as well. It is obtained by replacing the amplitudes in (26) by their elastic parts and removing the cutoff -the elastic contributions are convergent. WCM [7] instead represent the elastic part of the mass difference with two terms. 5 The sum of the two, δM el + δM sub el , differs from m el by Numerically, ∆m el is small: using the parametrization of Kelly [72], we obtain ∆m p el = −0.051 MeV, ∆m n el = −0.064 MeV. In the difference between proton and neutron, these numbers even partly cancel.
At the precision at which the nucleon form factors can nowadays be measured, it matters whether the standard expression for m el or the quantity m el + ∆m el is determined. For the decomposition (28) to be valid, it is essential that the nucleon form factors exclusively occur in m el -any other representation of the elastic part must be compensated by a corresponding correction in the term arising from the subtraction function.

Contributions from the subtraction function
As demonstrated in the preceding sections, the inelastic contributions to m QED are totally dominated by the one from the subtraction functionS. The differences in the values quoted for the elastic contributions are small compared to those from inelastic processes. Hence we can compare the various determinations of the mass difference that rely on dispersion theory by comparing the corresponding representations forS. 5 For a critical examination of their line of reasoning, we refer to Appendix E in [12] and to [13,14]. The bands labeled B and C in Fig. 9 show the models for the subtraction function of WCM [7] and ESTY [8], respectively. They are obtained from the representations proposed for S 1 in these references, merely converting numbers for S inel 1 into numbers forS by means of equation (99). The width of the bands exclusively shows the uncertainties arising from the experimental information used for the magnetic polarizabilities -those associated with the freedom in the choice of the model would widen it further. In the Q 2 range shown in the figure, both models are consistent with our analysis, but come with significantly larger errors (as the lower edge of band C runs within our band of uncertainties, it cannot be seen in Fig. 9).
The input used in models B and C for the value of the subtraction function at Q 2 = 0 is the same -it is based on the experimental determination of the polarizabilities of the nucleon. At small values of Q 2 , our uncertainties are smaller because the predictions obtained from Reggeon dominance for the polarizabilities of the neutron [12] are more precise than the experimental values. An improved measurement of the polarizabilities would be most welcome as it would subject Reggeon dominance to an important test. In this connection, we also refer to the new lattice results on the magnetic polarizabilities discussed in section 17.
At large values of Q 2 , the uncertainty band attached to model C is more narrow than the one of B, because the parametrization is improved: asymptotically, model C does reproduce the leading term in the operator product expansion of S inel 1 . As can be seen in Fig. 7, however, the nonleading spin 2 contributions disappear only ex-tremely slowly. In the parametrization of model C, these are neglected.
The net result obtained with model C for β p−n M = −0.5(1.6) is m C QED = 0.95 (25) MeV [8]. The corresponding outcome for the contribution from the subtraction function is obtained by removing the elastic part as well as those from the convergent dispersion integrals. With the entries for the elastic contributions listed in Table  I of [8] and the values given in equation (105) for the tiny terms m F2 and mF , this yields m C S = 0.19 (25) .
The value obtained by integrating the subtraction function of model C only over the low energy region is nearly the same: m C S (Q 2 < 2) = 0.20 (29) MeV. This indicates that in the evaluation of model C in [8], the contribution from Q 2 > 2 is nearly cancelled by the counter term, but we cannot verify this within our own framework. Since the parametrization of S 1 used in model C neglects the non-leading contributions in the asymptotic formula (100), it does not make sense to insert the corresponding representation forS in the expression (78) for mS -the integral diverges. Also, the blue line in Fig. 7 shows that for S inel 1 , asymptopia sets in extremely slowly, because the contributions generated by the short distance singularities of spin 2 fall off only logarithmically.
The numerical results for the subtraction function used by TWY and Tomalak are very similar to model C and they also lead to similar results for the e.m. part of the mass difference: m QED = 1.04(11) MeV (TWY) and m QED = 1.09 (30) MeV (Tomalak). The difference mainly arises from the input used for β p−n M . Note that the value β p−n M = −1.12 (40) used by TWY comes with a remarkably small error and disagrees with the Reggeon dominance prediction (85) by 2.5 σ. This is puzzling, because the determination of β p−n M in TWY is based on the lattice data of Blum et al. [76] -as shown in Fig. 8, these data are perfectly consistent with the range for m QED obtained from Reggeon dominance.
The ambiguities related to the fact that the function S inel 1 approaches asymptotics only very slowly do not arise if the parametrization of model C is used to representS rather than S inel 1 . We refer to this option as model D: the momentum dependence ofS is described by the function specified in equation (95), m 0 is identified with the scale m 2 0 = 0.71 GeV 2 occurring in the dipole representation of the nucleon form factors [10] and the parameter c 0 is fixed with the experimental valueS(0) = −0.2(2.6) GeV −2 given in equation (88), which is based on the determination of the polarizabilities in [64]. The blue shaded region in Fig. 9 shows that the subtraction function obtained with this variant of the models proposed in [7][8][9] agrees perfectly well with our analysis, but comes with a much larger error. Inserting the parametrization of model D in formula (78), we obtain The region Q 2 > 2 does not contribute much to the central value, but is responsible for a substantial fraction of the error: m D S (Q 2 > 2) = −0.02(28) MeV. C and D have the same behaviour at very small and very large values of Q 2 -they only differ in the form of the interpolation used in between. The example shows that -if only the leading terms in the OPE of the subtraction function are accounted for, the outcome is very sensitive to the form of the interpolation: replacing C by D lowers the central value of mS by 0.24 MeV and thus lowers the outcome for the central value of the mass difference to m QED = 0.71 MeV. This is within the uncertainty range attached to our result (106). The sensitivity to the form chosen for the interpolation arises because the subtraction function S 1 reaches asymptotics only very slowly.
Our analysis is not affected by this ambiguity, because we calculate the subtraction function in the region Q 2 < 3.5 on the basis of the experimental information about the structure functionF and rely on the theoretical information about the asymptotics only at higher energies. As pointed out in section 20, the contribution arising from the region Q 2 > 2 is nearly independent of the form of the parametrization used there, provided only that it obeys the theoretical constraints imposed by asymptotic freedom.

Summary and conclusions
1. Dispersion theory determines the amplitude in terms of its physical singularities (poles, cuts), provided the asymptotic behaviour is known. The use of amplitudes that contain kinematic zeros is best avoided, because these make it very difficult to sort out the asymptotic behaviour. We work with the invariant amplitudes introduced by Cottingham which do not contain such deficiencies and which we denote by T 1 , T 2 .
2. In the framework of dispersion theory, the elastic part of T 1 , T 2 is an unambiguous notion, determined by the requirement that it is analytic except for the poles generated by the elastic reaction and disappears when ν → ∞. Accordingly, the elastic contribution to the Cottingham formula is unambiguous.
3. As we do not know the error matrix occurring in the determinations of the form factors, we are not in a position to give a reliable estimate for the uncertainties in m el . We instead rely on the results obtained with the three different parametrizations in [71][72][73], which are covered by A determination of m p el , m n el and m p−n el on the basis of the information about the nucleon form factors available today would reduce the error considerably, but at the precision to which the inelastic contributions can currently be determined, the uncertainty quoted in (112) is too small to affect the error estimate attached to our result for m QED .
4. The leading terms of the operator product expansion of the Compton amplitude involve contributions arising from short distance singularities related to operators of spin 0 as well as spin 2. We make use of the fact that the leading spin 2 contributions to T 1 and T 2 only differ in normalization: in the combina-tionT ≡ T 1 + 1 2 T 2 , they drop out. Replacing the pair {T 1 , T 2 } by {T , T 2 } simplifies the analysis considerably. 5. A further simplification occurs if the dispersion relation forT is not subtracted at ν = 0, but at ν = 1 2 q 2 . This ensures that the contributions from the dispersion integrals over ImT and Im T 2 both contain the factor ν 2 − 1 4 q 2 . The point here is that in the Cottingham formula, only the angular average matters. Since the angular average of ν 2 − 1 4 q 2 vanishes, the contributions from the dispersion integrals are suppressed -numerically, these contributions are tiny. In our decomposition of the amplitude, only the elastic term and the integral over the subtraction function can generate significant contributions to the mass difference.
6. The quarks and gluons Reggeize. The exchange of Reggeons generates moving poles. For large values of ν at fixed q 2 , a Reggeon contributes withT ∝ ν α and T 2 ∝ ν α−2 , where α is the value of the trajectory α(t) at t = 0. Since there are trajectories with α > 0, the dispersion relation forT must be subtracted. The one for T 2 does not require a subtraction. 7. We assume that the asymptotic behaviour ofT is dominated by the contributions from the Reggeons, which we denote byT R . More precisely, we require that T −T R tends to zero when ν → ∞ and refer to this assumption as Reggeon dominance. A nonzero limiting value would represent a fixed pole -we are thus assuming that Reggeization is complete and only moving poles occur. Note that the dispersion relations forT and T 2 imply the presence of contributions that fall off with the power ν −2 . In T 2 , these contributions correspond to a fixed pole at α = 0 -Reggeon dominance is perfectly consistent with fixed poles of this sort. 8. Reggeon dominance implies a sum rule that determines the subtraction functionS in terms of the struc-ture functionF . The explicit expression given in (23) shows that neither the nucleon form factors nor the structure function F 2 enter. A variant of this sum rule was proposed by Elitzur and Harari, long ago [4], on the basis of duality and finite energy sum rules.
9. The value ofS(q 2 ) at q 2 = 0 is related to the polarizabilities of the nucleon. As is well known, the sum of the electric and magnetic polarizabilities is determined by a sum rule involving the cross section for photoproduction. Reggeon dominance implies separate sum rules for the electric and magnetic polarizabilities. The prediction obtained for the difference between the magnetic polarizability of proton and neutron [12] is in agreement with experiment, but this represents only a rather weak test of Reggeon dominance, because the uncertainties in the experimental result are rather large. The errors attached to the recent lattice result of [68] are much smaller -it is encouraging that Reggeon dominance passes this more stringent test as well. More work on the polarizabilities, particularly those of the neutron, would be most welcome. 10. Theory fixes the asymptotic behaviour of the subtraction function: if Q 2 becomes large,S tends to C/Q 4 , where the constant C is given by the proton matrix element of the operator 1 9 (4m u − m d )(ūu −dd). This also holds for S 1 (q 2 ) = T 1 (0, q 2 ), the subtraction function commonly used in dispersive analyses of the Compton amplitude, but the short distance singularities related to operators of spin 2 generate a significant difference in the asymptotic behaviour. Fig. 7 compares the momentum dependence of S 1 andS on a logarithmic scale and shows that, in contrast toS, the asymptotics of S 1 sets in only very, very slowly.
11. An important part of the calculation concerns the determination of the residue of the Reggeon with the quantum numbers of the f 2 , which dominates the asymptotic behaviour of the difference between the amplitudes of proton and neutron. Fig. 4 shows that the result obtained at low values of Q 2 from the Regge representation of [55] matches the outcome of the Regge fit to the numerical ABM table remarkably well.
12. With the values for the subtraction function obtained from the solution of the sum rule, our net result for the e.m. part of the mass difference between proton and neutron reads m p−n QED = 0.58 ± 0.16 MeV .
The conclusions reached in Ref. [5] are thus confirmed: m QED is dominated by the elastic contribution. The uncertainty in the result obtained forty five years ago, m QED = 0.7(3) MeV [5], is reduced by about a factor of two. In the present analysis, the uncertainty is predominantly due to the contributions from the resonance region above the ∆(1232). It could be reduced by an improved experimental determination of the structure functions in that region, particularly for the neutron. 13. It is by no means puzzling that the inelastic contributions are so small: (a) the angular integration suppresses the contributions from the dispersion integrals, (b) in the deep inelastic region, the subtraction function is nearly the same for proton and neutron -in the chiral limit, there is no difference, (c) in the region where Reggeon exchange dominates, the leading term, the Pomeron, is the same, (d) isospin symmetry ensures that the most important resonance, the ∆(1232), contributes equally to proton and neutron and (e) the leading terms of the chiral perturbation series are also the same.
With the experimental value of the mass difference, the above result implies that the part due to the difference between m u and m d is given by 14. The lattice results for these quanitities did not yet reach a level of coherence to be covered by the FLAG report, but the method is steadily being improved and, in the long run, should provide reliable numbers. Fig. 8 indicates that the most recent lattice values are larger than the outcome of the present work. If the value of m p−n QED should turn out to be larger than 1 MeV, we would have to conclude that the Compton amplitude does not fully Reggeize: the amplitude T 1 would then contain a fixed pole that invalidates the Reggeon dominance hypothesis. We would then be left with a puzzle: what is the physical origin of this fixed pole?
16. The evaluations of the Cottingham formula in [7][8][9] lead to values for m QED around 1 MeV. In these references, a simple algebraic ansatz is used to parametrize S inel 1 , the inelastic part of the subtraction function T 1 (0, q 2 ). Fig. 7 shows that, in contrast to these parametrizations, S inel 1 approaches asymptotics only extremely slowly.
The mismatch with the asymptotics disappears if the ansatz is assumed to be valid forS rather than S inel 1 . The central value obtained for m QED then drops by 0.24 MeV and winds up slightly below the elastic contribution, in agreement with what we find. On the other hand, quite apart from the sensitivity to the precise form of the assumptions underlying those models, the uncertainties in the result for m QED are much larger than ours, because the experimental determination of β p−n M , which plays a key role in that approach, is subject to large uncertainties.
The grace files used for the figures can be found in the attached ancillary files.

A Form of the Wilson coefficient for spin 2
Lorentz invariance implies that a tensor C µναβ (q) that is symmetric under µ ↔ ν and α ↔ β and only depends on the four-vector q µ is of the form: C µναβ (q) = a g µν g αβ + b {g µα g νβ + g να g µβ } (A.1) +c g µν q α q β + d g αβ q µ q ν + e {g µα q ν q β + g µβ q ν q α +g να q µ q β + g νβ q µ q α } + f q µ q ν q α q β , where a, b, c, d, e, f can only depend on q 2 . Current conservation, q µ C µναβ (q) = 0, fixes a, b, c in terms of d, e, f : a + q 2 d = 0 , b + q 2 e = 0 , c + 2e + q 2 f = 0 . (A.2) On account of the tracelessness of O αβ , the coefficient d drops out in the sum C µναβ (q)O αβ . To simplify the notation, we replace the coefficients e, f by c 2 , c 3 , with e = 1 2 c 2 , f = c 3 . The contribution from an operator of spin 2 to the OPE then takes the form: 6 C µναβ (q)O αβ = c 3 (q µ q ν − g µν q 2 )O αβ q α q β (A.3) This shows that, while Lorentz invariance and current conservation fix the Wilson coefficients belonging to operators of spin 0 in terms of a single function c 1 (q 2 ), those associated with operators of spin 2 involve two such functions: c 2 (q 2 ) and c 3 (q 2 ). 6 Hill and Paz [35] use a different normalization. In our notation, they work with T HP µν = 2T µν and normalize the spin 2 operator formed with the derivatives of the quark fields differently: O f 2 HP

B Operator product expansion for free quarks
For free quarks, the time-ordered product of two currents can be decomposed as where S f (z) is the quark propagator and N stands for normal ordering. In this expression, the singularities exclusively reside in the propagator -the matrix elements of the normal ordered products are regular at x = y.
The short distance expansion of the propagator starts with The leading singularity is contained in the first line of equation (B.1) and is proportional to (−z 2 +iǫ) −3 -the matrix elements thereof represent the disconnected part of the amplitude. We are interested in the singularities of the connected part, i.e. in the terms that contain one quark propagator and two quark fields. To analyze these, we set x = X + 1 2 z, y = X − 1 2 z and expand in powers of z. The expansion of the connected part starts with T j µ (x)j ν (y) = ǫ µναβ z α π 2 (−z 2 + iǫ) 2 f Q 2 ff γ β γ 5 f (B.3) We have dropped the normal ordering prescription as well as the argument of the quark fields -it is understood that the quark bilinears occurring here are to be normal ordered and evaluated at the point X.
The operator O f αβ contains components with spin 0, 1 and 2:

C Derivation of the sum rule
To calculate the limiting value of the dispersion integral (21), we first use partial fractions and split the integral into two parts: The reason for denoting the first term by −S a is that it is independent of ν and hence amounts to a contribution to the subtraction function.
To perform the limit in the remainder, we decompose the structure function into two parts withF = (F − F R ) + F R . In the differenceF −F R , the terms proportional to x 1−α cancel. We assume that the remainder is sufficiently smooth for x → 0, so that (F − F R )/x 2 is integrable and the integration can be interchanged with the limit. 7 This leads tō with ξ = Q 2 /(2mν). The integral J α represents a hypergeometric function. What remains to be done is to work out the behaviour of this function at small values of ξ. This can be done by making use of the known properties of the hypergeometric functions. Alternatively, one may observe that the contributions from the critical region x ∼ ξ remain the same if the integral is extended to infinity -it can then be done explicitly. On the interval x th < x < ∞, on the other hand, the limit can be interchanged with the integral, which can then be done explicitly as well. The result reads The first term is proportional to ν α . The Reggeon amplitudes specified in equation (18) have the same behaviour when ν becomes large. Indeed, one readily checks that the two expressions agree, so that If the singularity is more complicated, the sum rule does not get lost, but the explicit form must be adapted.
Collecting terms, the formula (22) yieldsS =S a +S b + S c . This agrees with the expression (23) forS quoted in section 5.