Forward doubly-virtual Compton scattering off the nucleon in chiral perturbation theory at NLO: the subtraction function and moments of unpolarized structure functions

The forward doubly-virtual Compton scattering (VVCS) off the nucleon contains a wealth of information on nucleon structure, relevant to the calculation of the two-photon-exchange effects in atomic spectroscopy and electron scattering. We report on a complete next-to-leading order (NLO) calculation of low-energy VVCS in chiral perturbation theory ($\chi$PT). Here we focus on the unpolarized VVCS amplitudes $T_1(\nu, Q^2)$ and $T_2(\nu, Q^2)$, and the corresponding structure functions $F_1(x, Q^2)$ and $F_2(x,Q^2)$. Our results are confronted, where possible, with"data-driven"dispersive evaluations of low-energy structure quantities, such as nucleon polarizabilities. We find significant disagreements with dispersive evaluations at very low momentum-transfer $Q$; for example, in the slope of polarizabilities at zero momentum-transfer. By expanding the results in powers of the inverse nucleon mass, we reproduce the known"heavy-baryon"expressions. This serves as a check of our calculation, as well as demonstrates the differences between the manifestly Lorentz-invariant (B$\chi$PT) and heavy-baryon (HB$\chi$PT) frameworks.


I. INTRODUCTION AND OUTLINE
The forward doubly-virtual Compton scattering (VVCS), Fig. 1, is not a directly observable process. Nevertheless, it is traditionally of high relevance in studies of nucleon and nuclear structure, and of their impact on atomic nuclei. At high energies the VVCS has the apparent connections to deep-inelastic scattering, whereas at low energies it is important for precision atomic spectroscopy, where it serves as input for calculations of the nuclear-structure corrections. Analytical properties of the VVCS amplitude are used to establish useful relations -sum rules -between the static (electromagnetic moments, polarizabilites) and dynamic (photoabsorption cross sections) properties of the nucleon [1][2][3][4][5], see also Refs. [6][7][8][9][10][11] for reviews. In the last decade, with the advent of muonic-atom spectroscopy by the CREMA Collaboration [12][13][14], the interest in nucleon VVCS has resurged in the context of the "proton radius puzzle" (see, e.g., Refs. [15,16] for reviews). The muonic atoms, being more sensitive to nuclear structure than conventional atoms, demand a higher quality of this input in both the Lamb shift [12,13] and, in the near future, the hyperfine structure measurements [17][18][19]. The VVCS enters here in the form of the two-photon exchange (TPE) corrections appearing at O(Z 4 α 5 ), which is the subleading order for the nuclear-structure effects in the Lamb shift (the leading being the charge radius), and leading in the hyperfine structure. In either case, the TPE is the leading theoretical uncertainty and precising this contribution is a challenge for the nuclear and hadron physics community.
In this work we focus on the unpolarized nucleon VVCS, described, for each nucleon (proton or neutron), by two scalar amplitudes T 1,2 (ν, Q 2 ), functions of the photon energy ν and virtuality Q 2 . The discontinuity of these amplitudes is given, respectively, by the two unpolarized structure functions F 1 (x, Q 2 ) and F 2 (x, Q 2 ).
To date, there are two approaches: 1) dispersion relations (DR) and 2) chiral perturbation theory (χPT), used for evaluation of nucleon VVCS, with the goal of quantifying the relevant corrections in muonic hydrogen. It is expected and highly desirable that 3) lattice QCD will join this effort in the near future. In the mean time, however, the DR approach is the most popular one. It employs the well-known dispersion relations expressing the VVCS amplitudes as integrals of the structure functions known empirically from inclusive electron scattering.
Unfortunately, the DRs determine the VVCS in terms of the structure functions only up to a "subtraction function" T 1 (0, Q 2 ). The latter function is not well-constrained em-pirically, which makes this approach prone to model uncertainties. It is worth to mention that there is a new proposal on how the subtraction can further be constrained via the dilepton electroproduction [20]. However, in the foreseeable future, this issue will preclude a systematic improvement of the theoretical uncertainty within the DR approach.
Here we employ the second approach. More specifically, we use an extension of SU (2) χPT [21,22] to the single-baryon sector [23][24][25], referred to as the baryon χPT (BχPT), augmented by inclusion of the explicit ∆(1232)-isobar in the δ-counting scheme [26]. In this framework we compute the inelastic (non-Born) part of the VVCS amplitudes to next-toleading order (NLO). A first version of this calculation was briefly considered in Ref. [27].
Here we provide a few important improvements, in particular, the inclusion of the Coulombquadrupole (C2) N → ∆ transition, and a more comprehensive comparison of our results with the DR approach. The impact of this calculation on the muonic-hydrogen Lamb shift, extending our previous evaluation [28] to higher orders, will be discussed elsewhere.
The paper is organized as follows. In Sec. II, we recall the general formulae for VVCS and its relation to structure functions, form factors and polarizabilities. In Sec. III, we discuss the main ingredients of our NLO calculation. In Sec. IV, we examine results for the proton and neutron scalar polarizabilities, and some of the other moments of structure functions. In the concluding section (Sec. V), we summarize and give a brief outlook for the near-future work. In App. A, we discuss the structure functions, in particular, the πN , π∆ and ∆ production channels relevant to our calculation. In App. B, we give analytical expressions for the πN -loop and ∆-exchange contributions to the static values and slopes of the polarizabilities and moments of structure functions. The complete expressions, also for the π∆-loop contributions, can be found in the Supplemented material. Figure 1 schematically shows the VVCS amplitude, which for an unpolarized target (of any spin) can be decomposed into two independent Lorentz-covariant and gauge-invariant tensor structures [9]:

II. VVCS FORMALISM
where p and q are the four-momenta of the target particle and the photon, respectively; M N is the target (here, nucleon) mass. The scalar amplitudes T i are functions of the photon energy ν = p · q/M N and virtuality Q 2 = −q 2 .
The optical theorem relates the absorptive parts of the VVCS amplitudes to the structure functions, or equivalently, the inclusive electroproduction cross sections: with the fine-structure constant α = e 2 /4π, and the Bjorken variable x = Q 2 /2M N ν. The two response functions σ T and σ L are cross sections of total photoabsorption of transversely (T) and longitudinally (L) polarized photons. The flux of virtual photons is conventionally defined up to the flux factor K(ν, Q 2 ). The experimental observables do not depend on it, only the definitions of the response functions σ T and σ L do. Throughout this work we adopt Gilman's flux factor (for other common choices, cf. Ref. [6]): where q is the photon three-momentum in the lab frame. The VVCS amplitudes satisfy the following dispersion relations derived from the above statement of the optical theorem, combined with general principles of analyticity and crossing symmetry (cf., for example, Refs. [6,9,10] for details): We refer to M 1 (Q 2 ) as the generalized Baldin sum rule, see Sec. IV A. More generally, we have the following relation (for an integer n): arising from electromagnetic gauge invariance. One way to derive it is to introduce the longitudinal amplitude and to show that lim Q 2 →0 T L (ν, Q 2 ) = 0. Incidentally, the same is true for asymptotically large Q 2 , because of the Callan-Gross relation: In what follows, we consider some of these moments, obtaining them from the χPT results for the VVCS amplitudes, and compare them with the results of empirical parametrizations of the nucleon structure functions. The dispersion relations (4) are used by us to cross-check the results, using the tree-level photoabsorption cross sections discussed in App. A.

III. CALCULATION OF THE VVCS AMPLITUDE AT NLO
Our goal here is to obtain the χPT predictions of the non-Born parts of the nucleon VVCS amplitudes T 1,2 . The present NLO calculation is still within the "predictive powers" of χPT for Compton scattering (CS) amplitudes, i.e., the results are given in terms of wellknown parameters (see Table I) obtained from non-Compton processes. In this sense, it is complementary to the existing calculations of real CS (RCS) [32,33], and the virtual CS (VCS) [34]. All of these studies, including the present one, are done in the same framework, using the same set of parameters.

A. Remarks on power counting
We shall employ BχPT, which is the manifestly-covariant extension of χPT to the singlebaryon sector in its most straightforward implementation, where the nucleon is included as in Ref. [23]. The power-counting concerns raised in Ref. [23] have been overcome by renormalizing away the "power-counting violation" using the low-energy constants (LECs) available at that order. This has been shown explicitly within the "extended on-mass-shell renormalization scheme" (EOMS) [25], but is not limited to it. The inclusion of the explicit ∆(1232) here will follow the "δ-counting" framework of Ref. [26] (see also Refs. [35,36] for concise overviews).
To explain the power counting in more detail, let us recall that chiral effective-field theory is based on a perturbative expansion in powers of pion momentum p and mass m π over the scale of spontaneous chiral symmetry breaking Λ χ ∼ 4πf π , with f π 92 MeV the pion decay constant. Each operator in the effective Lagrangian, or a graph in the loopwise expansion of the S-matrix, can have a specific order of p assigned to it.
To give a relevant example consider the following operator: with δβ the coupling constant, N the Dirac field of the nucleon, and F 2 the square of the electromagnetic field strength tensor, . This is an operator of O(p 4 ). Two of the p's come from the photon momenta which are supposed to be small, and the other two powers arise because the two-photon coupling to the nucleon must carry a factor of α (the charge e counts as p, since we want the derivative of the pion field to count as p even after including the minimal coupling to the photon). This operator enters the effective Lagrangian with an LEC, which we denote δβ. It gives a contribution to the CS amplitude in the form of 1 and leads to a shift in the magnetic dipole polarizability as: β M 1 → β M 1 + δβ. Now, two remarks are in order.
i) Naturalness. The magnitude of the LEC is not arbitrary. It goes as δβ = (α/Λ 3 χ )c, with the dimensionless constant c being of the order of 1, or more precisely: This condition ensures that the contribution of this operator is indeed of O(p 4 ), as inferred by the power counting.
ii) Predictive powers. This LEC enters very prominently in the polarizabilities and CS at tree level, which means its value is best fixed by the empirical information on these quantities. If this is so, the O(p 4 ) result is not "predictive", as it could only be used to fit the χPT expression to experiment or lattice QCD calculations. On the other hand, contributions of orders lower than p 4 are predictive, as they only contain LECs fixed from elsewhere.
As already mentioned, the "predictive" contributions to CS and polarizabilities have been identified and computed for the case of RCS [32], VCS [34], and VVCS [27]. Our present calculation is quite analogous to those works and hence we refer to them for most of the technical details, such as the expressions for the relevant terms of the effective Lagrangian. It is crucial to first study these predictive contributions. We note, however, that here we choose to also include the p 4 LEC that shifts the magnetic polarizability. In doing so, we fit the value of δβ so as to reproduce the Baldin sum rule values: 1 Throughout this paper we use the conventions summarized in the beginning of Ref. [9].
taking the values of α E1 obtained at O(p 4 /∆) as BχPT predictions. This choice reflects the fact that the most prominent scalar moments considered here, the second moment of F 1 (x, Q 2 ) and the first moment of F 2 (x, Q 2 ), both have the Baldin sum rule as their static limit. The values of the magnetic polarizabilites that result from this fit are where the error bar does not include the theoretical uncertainty. One has to admit that this procedure results in somewhat smaller values of β M 1 than, for instance, those obtained in the recent HB and covariant chiral analyses: 3.2(0.5) × 10 −4 fm 3 [39,40] for the proton and 3.65(1.25) × 10 −4 fm 3 [41] for the neutron. We will, however, use this simplified procedure since the only affected quantity studied by us is the proton subtraction function . . , and the discrepancy for β M 1p is tolerable.
We also include the Coulomb-quadrupole (C2) N → ∆ transition, described by the g C term in the following non-minimal γ * N ∆ coupling [42,43] (note that in these references the overall sign of g C is inconsistent between the Lagrangian and Feynman rules): with M + = M N + M ∆ and the dual of the electromagnetic field strength tensorF µν = 1 2 µνρλ F ρλ . The electric, magnetic and Coulomb couplings (g E , g M and g C ) are known from the analysis of pion photoproduction P 33 multipoles [43]. The corresponding numerical values, as well as those of other physical constants used in this work, are given in Table I. The Coulomb coupling is subleading compared with the electric and magnetic couplings, and it was not included in the previous calculations. However, the relatively large magnitude of g C hints at its potential numerical importance, which we examine in this work.
The counting of the ∆(1232) effects is done within "δ-counting" [26], where the Deltanucleon mass difference, For the non-Born VVCS amplitudes and polarizabilities the predictive orders are O(p 3 ) and O(p 4 /∆). The O(p 3 ) contribution comes from the pion-nucleon (πN ) loops. We refer to it here as the LO contribution. The ∆-exchange graph is described by the γ * N ∆ interaction in Eq. (18). For the magnetic coupling, one assumes a dipole behavior to mimic the form expected from vector-meson dominance: with the dipole mass Λ 2 = 0.71 GeV 2 . The inclusion of this Q 2 dependence is motivated by observing the importance of this form factor for the correct description of the electroproduction data [42]. Since electroproduction off the nucleon is directly connected with the nucleon polarizabilities via sum rules, it is understandable that a better description of the electroproduction data will lead to a better description of the Q 2 behavior of the polarizabilities. The effect of the dipole form factor in g M is illustrated in Figs. 4, 6, and 7.
A feature of the δ-counting is that the characteristic momentum p distinguishes two regimes: the low-energy (p m π ) and resonance (p ∆) regimes. The above counting is limited to the low-energy regime. Since we are interested in the low-energy expansion of the VVCS amplitudes (i.e., the expansion in powers of small ν with Q 2 finite), we do not consider the regime where one-Delta-reducible graphs are enhanced (resonance regime). However, going to higher Q 2 one does need to count the Delta propagators similar to the nucleon propagators, which, in turn, calls for inclusion of π∆ loops with two and three Delta propagators, which have been omitted here. They are only included implicitly by adjusting the isospin coefficients of the one-nucleon-reducible π∆-loop graphs to restore current conservation, as explained in the next section. Apart from that, π∆ loops have a rather mild dependence on momenta and the missing loops are unlikely to affect the Q 2dependence of the moments of structure functions significantly, even for Q 2 comparable to ∆ 2 .  [44] at the order they first appear.
The πN ∆ coupling constant h A is fit to the experimental Delta width and the γ * N ∆ coupling constants g M , g E and g C are taken from the pion photoproduction study of Ref. [42]. The free parameters δβ p,n are fitted to the Baldin sum rule for the proton and neutron [37,38], respectively.

B. Renormalization
The calculation of the πN -and π∆-loop graphs is analogous to Ref. [32], with the obvious extension to the case of a finite photon virtuality. The renormalization is also done in the exact same way; namely, subtracting the loop contribution to the Born term of the the VVCS amplitude. The π∆-loop graphs still contain divergences after this subtraction. These divergences are of higher orders, O(p 5 /∆ 2 ) and O(p 4 ), and will be cancelled by the corresponding higher-order contact terms. In practice, they are removed by taking the MS values of the divergent quantities.
As mentioned above, π∆-loop graphs where photons couple minimally to the Delta contain more than one Delta propagator and therefore should be suppressed by extra powers of p/∆. However, their lower-order contributions are important for electromagnetic gauge invariance and therefore for the renormalization procedure. In particular, the lower-order contributions of chiral loops should not affect the result of the low-energy theorem [45,46], and this condition is automatically satisfied by a subset of graphs if it obeys gauge invariance. The π∆-loop graphs, cf. Ref. [27, Fig. 2], form such a subset for the case of the neutral Delta. In reality, the Delta comes in four charge states (isospin 3/2), thus, a gauge-invariant set will in addition have the higher-order graphs where the photon couples minimally to the Delta. To make the subset of π∆-loop graphs gauge invariant without the higher-order graphs, we used the same procedure as in Ref. [32]. The one-particle-irreducible (1PI) graphs in Ref. [27,Fig. 2] are computed with the correct isospin factors, i.e., summing over all charge states of the Delta, whereas the isospin factors for the one-particle-reducible (1PR) loop graphs are chosen such that their ratio to the isospin factors of 1PI graphs is the same as in the neutral Delta case. This procedure automatically ensures exact gauge invariance and thus effectively includes the relevant contributions of the one-loop graphs with minimal coupling of photons to the Delta. If the latter graphs are included explicitly, the isospin factors of 1PR graphs can be restored to actual values. The corresponding change in the values of the polarizabilities, however, is of higher order than what we consider in this work.

C. Uncertainty estimate
To estimate the uncertainties of our NLO predictions, we define the running expansion parameterδ such that the N 2 LO is expected to be of relative sizeδ 2 [42]. To estimate the uncertainty of a polarizability P (Q 2 ) due to the neglected higher-order terms in the chiral expansion, we separate that polarizability into the static piece P (0) and the Q 2 -dependent remainder P (Q 2 ) − P (0). The uncertainty of P (Q 2 ) is obtained by adding the estimates for these two parts in quadrature: The uncertainties in the values of the parameters have a much smaller impact compared to the truncation uncertainty and are therefore neglected.

IV. RESULTS AND DISCUSSION
We now consider the numerical results for some of the moments of the nucleon structure functions which appear in the expansion Eq. (9). We shall also consider the proton subtraction function T 1 (0, Q 2 ). The complete NLO values will be decomposed into three individual contributions: the πN loops, the ∆ exchange, and the π∆ loops. In practice, we extract all results from the calculated non-Born VVCS amplitudes. For a cross-check, we used the photoabsorption cross sections described in App. A.
A. M The electric and magnetic dipole polarizabilities, α E1 (Q 2 ) and β M 1 (Q 2 ), encode information about the dipole response of the nucleon to an electromagnetic field. For finite momentum transfers, the sum of dipole polarizabilities is given by the generalized Baldin sum rule: where ν 0 is the lowest inelastic threshold, in this case the one-pion production threshold ν 0 = m π + (m 2 π + Q 2 )/2M N , and x 0 = Q 2 /2M N ν 0 . The electric and magnetic dipole polarizabilities of the nucleon enter the nucleon-structure contributions to the Lamb shift of muonic hydrogen and other muonic atoms [28,[52][53][54], and thus are of major interest for an accurate extraction of the nuclear charge radii.
Our BχPT predictions for α E1 + β M 1 are shown in Fig. 3 {upper panel}, both for the proton and the neutron, up to photon virtualities of 0.3 GeV 2 . Our main result is given by the blue solid lines and the blue error bands, where we used the p 4 LEC δβ to fit the static polarizabilities to the empirical Baldin sum rule values (green and purple dots) given in Eq. (16), see discussion in Sec. III A. Inclusion of δβ, cf. Table I, merely leads to a constant shift, as can be seen by comparing to the pure O(p 4 /∆) predictions (blue longdashed lines), which include the πN -loop, the ∆-exchange and the π∆-loop contributions. In order to illustrate the effect of the Delta in these predictions, we also plot the LO πN -loop contributions separately (red solid lines). We compare our results for the Q 2 evolution with the LO HBχPT predictions [47] (purple dashed lines) and the MAID model predictions [48,49] (black dotted lines). The latter are based on the generalized Baldin sum rule (22) [47]. The black dotted line is the MAID model prediction [48][49][50]; for the proton we use the updated estimate from Ref. [6] that includes the π, η, ππ channels. At Q 2 = 0 GeV 2 , we show the Baldin sum rule value for the proton (green dot) [37] and neutron (purple dot) [38]. For the proton, the Q 2 = 0. evaluated with (π + η + ππ) photoproduction cross sections [6]. The data points are also evaluations of the (generalized) Baldin sum rule [37,51,55]. One can see that the BχPT predictions seem to systematically overestimate the MAID model in the Q 2 range shown here. One has to note that the MAID model, on the other hand, slightly underestimates the empirical Baldin sum rule evaluations. The LO HB results seem to agree with the empirical values at the real-photon point [55] both for the proton and the neutron. However, they do not fall off with increasing Q 2 in contrast to the BχPT predictions. This asymptotic behavior is the reason for the large proton-polarizability effect on the muonic-hydrogen Lamb shift found within HBχPT [47,56,57], much larger than the phenomenological value. As shown in Ref. [28,31], this issue is solved within the relativistic formulation, which gives a result closer to calculations based on the dispersive approach.
The static dipole polarizabilities α E1 and β M 1 have been studied within both the HB and the BχPT. While HBχPT gives results remarkably close to the experimental determinations already at LO [58], the contribution of the ∆(1232) is harder to accommodate in this framework [59]. In contrast to that, LO BχPT [60,61] yields smaller values for the sum of dipole polarizabilities, in disagreement with the empirically extracted values based on evaluations of the Baldin sum rule with modern photoabsorption data [37,55,62]. However, the NLO contributions from ∆ exchange and π∆ loops improve the situation [32,33]. In the case of the proton, they bring the BχPT result in agreement with the experimental extraction, while for the neutron the total result is slightly bigger. The ∆(1232) contributions are, therefore, naturally accommodated in BχPT, and not in HBχPT (where they can be reconciled with the empirical values only by means of the p 4 LEC δβ, see, e.g., Refs. [40,63] for the recent calculations and review).
The BχPT contributions from πN loops, ∆ exchange and π∆ loops to the static polarizabilities are, in that order and in the usual units of 10 −4 fm 3 :  Fig. 4 {upper panel}. For the proton, the dominant contribution in the studied Q 2 range is that of the ∆ exchange, while for the neutron the πN -loop and ∆-exchange contributions are of roughly the same size. The importance of the Delta is related to the fact that the nucleon-to-Delta transition is dominantly of the magnetic dipole type, therefore it gives a huge contribution to β M 1 .
In addition to the static values, we can now investigate the slopes of the polarizabilities at the real-photon point. Decomposing the results as before into the three contributions, we observe that BχPT predicts large contributions to the slopes both from πN loops and ∆ exchange. The Q 2 dependence generated by π∆ loops, on the other hand, is negligible, as can be clearly seen from The dipole form factor in the magnetic coupling g M generates the Q 2 fall-off of the dipole polarizabilities, cf. Fig. 4, which is also observed in parametrizations of experimental cross sections [64]. Due to cancellations between the πN -loop and the ∆-exchange contributions, the dipole also crucially affects the overall sign of the slope, as can be seen in Fig. 4. Note that due to these cancellations we estimate the relative error of the slope byδ instead ofδ 2 . Evaluating the Baldin sum rule radius, we obtain r (α+β)p = 0.29(9) fm and r (α+β)n = 0.52 (16) fm, where we estimated the relative error to beδ. Here, we again used our result including the δβ contribution, i.e., we fixed the static polarizabilities to the Baldin sum-rule values in Eq. (16), while the slope is still a prediction of BχPT. The result for the proton is in tension with the sum-rule evaluations [51,64,65] which use empirical parametrizations of the structure function F 1 (x, Q 2 ), e.g. [64]: From Fig. 3 one can see that the MAID empirical parametrization also leads to a steeper slope than BχPT. This calls for a careful revision of the low-momentum behavior of the empirical parametrizations in the near future.
B. α L (Q 2 ) -the longitudinal polarizability The low-energy expansion of the longitudinal VVCS amplitude goes as with α L called the longitudinal polarizability. Note that, in terms of the moments α L = M 1 (0). The generalized longitudinal polarizability is given by, with Our BχPT prediction for α L (Q 2 ) is shown in Fig. 3 {lower panel}, where we compare our results, with and without the Delta contributions, with the MAID model predictions [6,[48][49][50] and the HB limit of the πN -loop contribution. One can see that the Delta plays a negligible role in the low-Q 2 evolution of α L , which in the BχPT approach is dominated by πN loops. Our results run very close to the MAID curves, with small discrepancies in the intermediate Q 2 region. At higher virtualities, these discrepancies decrease. The HB approach, on the other hand, seems to systematically overestimate the value of α L in the considered Q 2 range. This relatively big mismatch can be traced back to the slow convergence of the 1/M N expansion, as one can see from the analytic expression for the πN -loop contribution to α L (Q 2 = 0) given in App. B. For the static values of α L , we obtain the following contributions from πN loops, ∆ exchange and π∆ loops, in units of 10 −4 fm 5 : The corresponding individual contributions to the Q 2 dependence of α L (Q 2 ) are demonstrated in Fig. 4 {lower panel}. One again notices that ∆ exchange and π∆ loops give negligible contributions in this Q 2 range. The smallness of the ∆-exchange contribution is explained by the fact that the magnetic coupling g M does not contribute to α L .
At Q 2 = 0, the first moment of the structure function F 2 (x, Q 2 ), reproduces the Baldin sum rule: M 1 (0), cf. Eq. (23). However, at finite Q 2 this moment is independent of M  Fig. 3 {upper panel}, one can indeed see that the two moments noticeably diverge as one departs from the static limit. Note that the contribution of the p 4 operator in Eq. (14) simultaneously shifts M at Q 2 = 0, we find the following contributions from πN loops, ∆ exchange and π∆ loops, in units of 10 −4 fm 5 : Interestingly, the slope does not show such a drastic cancellation between the πN -loop and the ∆-exchange contribution as one encounters in the generalized Baldin sum rule. Correspondingly, the shape of the M   The red line represents the LO BχPT result, while the purple dashed line is the LO HB limit [47].
The black dotted line is the MAID model prediction [48][49][50]. At Q 2 = 0 GeV 2 , we show Baldin sum rule evaluations for the proton (green dot) [37] and neutron (purple dot) [38]. Lower panel: Fourth-order generalized Baldin sum rule for the proton (left) and neutron (right). The NLO BχPT prediction of this work is shown by the blue solid line and the blue band; the legend for the remaining curves is as in the upper panel. At Q 2 = 0 GeV 2 , we show fourth-order Baldin sum rule evaluations for the proton (dark green pentagon) [37] and neutron (magenta pentagon) [66].  Let us now consider the fourth moment of the structure function F 1 (x, Q 2 ): In the real-photon limit, this moment is related to a linear combination of dispersive and quadrupole polarizabilities [67,68], resulting in the fourth-order Baldin sum rule (see Ref. [9] for review): M Here we obtain the following NLO results for the proton and neutron (showing also the separate contributions from πN loops, ∆ exchange and π∆ loops), in units of 10 For the slopes at Q 2 = 0, we find, in units of 10 −4 fm 7 : dM (4) 1n (Q 2 ) dQ 2 The corresponding individual contributions to the Q 2 dependence of M 1 (Q 2 ) are demonstrated in Fig. 6 {lower panel}.
In Fig. 5 {lower panel}, we show our BχPT predictions compared to the MAID model predictions [48][49][50], the HB limit of the πN -loop contribution [47], and empirical evaluations of the fourth-order Baldin sum rule [37,66] for proton and neutron, respectively. Our NLO BχPT predictions are in good agreement with MAID, while the HB results fail to describe the decrease with growing Q 2 .

the proton subtraction function
The knowledge of the proton subtraction function T 1 (0, Q 2 ) is needed to evaluate the leading contribution of the nucleon structure to the (muonic-)hydrogen Lamb shift, see Refs. [9,15,16] for reviews. At very low momenta the non-Born part of the subtraction function is given by the magnetic dipole polarizability, T 1 (0, Q 2 )/Q 2 = 4πβ M 1 + O(Q 2 ). Since the Lamb-shift integrals are weighted towards low Q 2 , the low-momentum features of T 1 (0, Q 2 ) have a more pronounced effect. In particular, the uncertainty in the (empirical) extraction of β M 1 contributes the bulk of the theoretical uncertainty in Ref. [69]. At the same time, the slope of this function at Q 2 = 0 could potentially be important, and the different models and mechanisms could lead to rather different values of that slope. This is illustrated in Fig. 7 {left panel}, which shows the low-Q 2 behavior of T 1 (0, Q 2 )/4πQ 2 . One can see that both β M 1 and the slope change significantly when one adds the ∆ contributions (∆-exchange and the π∆ loops contribution), cf. Fig. 7 {right panel}. The resulting curve in the left panel is compared with the HBχPT evaluation of Ref. [69], showing an appreciable disagreement in the slope, with the Q 2 dependence of the two curves being noticeably different. Our NLO prediction of the slope at Q 2 = 0 is given by (in units of 10 −4 fm 5 ): 1 8π Extracting the slope of the subtraction function experimentally should in principle be possible through dilepton electroproduction as proposed in Ref. [20]. It remains to be seen whether such a measurement is feasible in the near future.

V. SUMMARY AND CONCLUSIONS
We have completed the NLO calculation of the unpolarized VVCS amplitudes in SU(2) BχPT, wiht explicit ∆(1232). We have calculated the non-Born amplitudes, which at this order come out as a parameter-free prediction of BχPT. These amplitudes are used to examine several notable combinations of the (generalized) polarizabilities that are expressed through the moments of the nucleon structure functions, i.e.: M   2 (Q 2 ) -the first moment of the structure function F 2 (x, Q 2 ). The dispersion relations between the VVCS amplitudes and the tree-level photoabsorption cross sections served as a cross-check of these calculations.
These results can be compared with the dispersive evaluations using the empirical parametrization of the nucleon structure functions. The biggest discrepancy is observed for the low-Q behavior of the generalized Baldin sum rule, calling for a future revision of the low-momentum behavior of the empirical parametrization of the structure function F 1 (x, Q 2 ).
Concerning the ∆(1232) contribution, we have seen that it plays an important role in transverse quantities, whereas in the longitudinal quantities, such as the longitudinal polarizability α L , its role is negligible. The Coulomb γ * N → ∆ transition coupling g C has generally a small effect in the unpolarized moments considered here.
We have also obtained an NLO prediction for the proton subtraction function T 1 (0, Q 2 ), which cannot be deduced from dispersion relations. This is an important step towards a systematic improvement of the LO χPT evaluation [28] of the proton polarizability contribution to the muonic-hydrogen Lamb shift. The ∆-production cross sections are related to the tree-level ∆-exchange shown in Fig. 2. The threshold for production of the ∆(1232)-resonance is at lab-frame photon energies of: Therefore, the ∆-production cross sections contain to the following Dirac's δ-function: δ(ν − ν ∆ ).
The explicit form of these cross sections is given by:   It is important to note that the above cross sections only describe the ∆-pole contributions to the tree-level ∆ exchange. In general, the VVCS amplitudes described by the ∆-exchange diagram in Fig. 2 can be split as follows [71]: (0, Q 2 ) is the usual subtraction function: with ω ± = (M 2 ∆ − M 2 ± Q 2 )/2M ∆ . T ∆-pole i are the ∆-pole contributions which feature a pole at the ∆(1232)-production threshold, and thus, are proportional to: where s and u are the usual Mandelstam variables. T ∆-exch.
i are the (∆-)non-pole terms in which the pole has cancelled out [71]: T ∆-exch.
To describe the non-pole terms in Eqs. (A6a) and (A6b) within the standard dispersive framework, Eq. (4), we define the auxiliary structure functions: