Finite volume corrections to forward Compton scattering off the nucleon

We calculate the spin-averaged amplitude for doubly virtual forward Compton scattering off nucleons in the framework of manifestly Lorentz invariant baryon chiral perturbation theory at complete one-loop order $O(p^4)$. The calculations are carried out both in the infinite and in a finite volume. The obtained results allow for a detailed estimation the finite-volume corrections to the amplitude which can be extracted on the lattice using the background field technique.


I. INTRODUCTION
Recent years have observed a rapidly increasing interest in calculations of nucleon structure functions on the lattice. Different algorithms, which enable one to extract these quantities from lattice measurements, have been proposed. For example, in Ref. [1] a method for a direct calculation of the quark and gluon distribution functions on Euclidean lattices by Lorentz-boosting of the nucleons was suggested. Recently, a lattice calculation of the Euclidean four-point function, describing the virtual Compton amplitude, and its relation to the leptoproduction cross section has been considered [2]. In the present paper we shall concentrate on an alternative proposal which is based on the use of the background field technique (or, the Feynman-Hellmann method) for measuring the forward doubly virtual Compton scattering amplitude off nucleons, see Refs. [3][4][5][6][7][8][9]. This amplitude is directly related to the moments of the structure functions. For a review of the present status of lattice studies of the structure functions, see, e.g., Ref. [10].
In Refs. [3][4][5][6][7][8][9] a comprehensive theoretical assessment of the feasibility of the extraction of the Compton amplitude has been carried out. Here, one has to note that a similar technique has been already successfully used for the extraction of the magnetic moments and polarizabilities of certain hadrons [11][12][13]. The study of Compton scattering, however, implies another level of sophistication. Namely, whereas the static characteristics of the nucleon can be measured in constant background magnetic and electric fields, the dependence of the forward Compton amplitude on the photon virtuality, q 2 = −Q 2 , cannot be studied similarly. Therefore, one has to use periodic background fields (in space), which enable one to consider non-zero values of the photon three-momentum, while the time-component of the photon momentum q stays zero. Several subtle issues had to be addressed in this context, for example, a consistent realization of the periodic background field on a finite lattice [14], or the zero-frequency limit [9]. It must also be mentioned that, according to Ref. [15], the interpretation of the lattice measurements, which are done in a finite volume, might be ambiguous both for constant and periodic fields. More precisely, the quantity that is obtained as a result of such a measurement (for instance, the polarizability) could be different from what one has previously identified as a finite-volume counterpart of the polarizability. This point of view has been countered in Ref. [9], where it has been argued that the finite-volume lattice results allow for a unique interpretation in terms of well-defined physical quantities (at least when the photon virtuality is not zero). Since this issue is important in the context of the problem considered here, we shall briefly address it later.
It should be stressed that the measurement of the forward Compton amplitude on the lattice is a useful endeavor by itself, even beyond its relation to the nucleon structure functions. Indeed, let us point out that the forward Compton amplitude represents an important building block in many long-standing fundamental problems that have recently come under a renewed scrutiny. In particular, the knowledge of this amplitude is needed for the evaluation of the Lamb shift in muonic hydrogen [16], as well as the proton-neutron mass difference. The study of the latter problem has a decades-long history [17,18], but still continues to attract quite some interest that is reflected in a string of recent publications [19][20][21][22][23][24]. To a large part, this upsurge of interest can be related to the fact that the present lattice arXiv:2010.10917v1 [hep-lat] 21 Oct 2020 studies are in a position to address the calculation of the purely electromagnetic proton-neutron mass shift in QCD plus QED and hence the results of phenomenological determinations can be directly confronted with lattice data. Further, in Refs. [23,24], under the assumption that the high-energy behavior of the Compton amplitude is fully determined by Reggeon exchange (the so-called Reggeon dominance hypothesis), a sum rule has been derived that involves this amplitude in a particular kinematics (a variant of this sum rule has been known in the literature already for fifty years [25]). Notably, the latter enables one to express the Compton amplitude through the experimentally measured electroproduction cross sections. Calculations on the basis of the above sum rule have been performed recently [24], where the uncertainties emerging from the use of all presently available experimental input have been thoroughly analyzed. A direct evaluation of the Compton amplitude on the lattice would allow one to compare the outcomes of these two different theoretical calculations. Should it happen that the results are very different, this could be attributed to the failure of the Reggeon dominance hypothesis, i.e., to the existence of so-called fixed poles in the Compton amplitude. At present, we are not aware of any mechanism within QCD that would lead to such poles. Hence, their discovery would challenge our understanding of the asymptotic behavior of QCD and stimulate a quest for new mechanisms, which are responsible for this behavior.
One of the most important questions, which so far has not been addressed in the context of the extraction of the Compton amplitude from lattice data, is the issue of the finite-volume corrections to the physical quantities of interest. It is very important to estimate, prior to performing any calculations on the lattice, how large lattices should be used to suppress the unwanted finite-volume artifacts. Note that, even though on general grounds these artifacts are exponentially small, due to possible large prefactors they might be still substantial for the presently used lattice sizes. The systematic study of this problem that is carried out in what follows within the framework of Baryon Chiral Perturbation Theory (BChPT) at order O(p 4 ) is intended to fill this gap.
The layout of the paper is the following. In Sect. II we collect all definitions and input, which will be needed later. This concerns both purely infinite-volume calculations as well as the finite-volume setting used on the lattice for the extraction of the Compton amplitude. Further, the calculation of the Compton amplitude is carried out in the infinite as well as in a finite volume. Namely, Sect. III contains the full expression of the infinite-volume Compton amplitude at O(p 4 ) in BChPT. Also, a comparison to the results available in the literature is carried out. The expression of the finite-volume amplitude at the same order is given in Sect. IV. In Sect. IV B, the results of the numerical estimations of the finite-volume artifacts are discussed. Sect. V contains our conclusions.

II. BASIC DEFINITIONS AND NOTATIONS
A. Doubly-virtual Compton scattering in forward direction in the infinite volume In this paper we follow the notations of Ref. [20]. In order to render the paper self-contained, below we collect all formulae that will be used in the infinite-volume calculations. The Compton scattering amplitude is defined as: T µν (p , s , q |p, s, q) = i 2 d 4 x e iq·x p , s |T j µ (x)j ν (0)|p, s , where (p , s ) and (p, s) are the four-momenta and spin projections of incoming and outgoing nucleon, respectively, and q and q are the momenta of the (virtual) photons in the initial and final state, respectively. Further, j µ denotes the electromagnetic current. The state vectors of the nucleon are normalized as: We define the unpolarized forward scattering amplitude as an average over the nucleon spins: Using Lorentz-invariance, current conservation and parity, this amplitude can be expressed through two invariant amplitudes: where ν = p · q/m and m is the nucleon mass.
At this place we mention that the choice of the invariant amplitudes is not unique. In the literature, another choice is often made, using the setT 1 ,T 2 withT 1 = q 2 T 1 + ν 2 T 2 andT 2 = −ν 2 T 2 . This alternative choice, however, produces kinematic singularities, which complicate the discussion of the asymptotic behavior of these amplitudes. As a result, the issue with the fixed poles may become obscure. For further details on this subject we refer the reader to Refs. [20,24] and also to [26,27]. Further, in Refs. [23,24], another set of invariant amplitudesT = T 1 + 1 2 T 2 , T 2 has been introduced instead of T 1 , T 2 . The advantage of using this set consists in the fact that the leading asymptotic behavior ofT at large values of Q 2 is governed by spin-0 operators, whereas the contribution from the spin-2 operators in the operator product expansion cancels in this particular linear combination. A thorough discussion of this question is given in Ref. [24]. Here, we only mention that the setT , T 2 is obviously free from the kinematic singularities, as well as the set T 1 , T 2 .
Further, the invariant amplitudes can be split into the elastic and inelastic (or, equivalently, into the Born and non-Born) parts. Again note that such a splitting is not uniquely defined and the definition of Ref. [20] differs from the ones used in Refs. [28,29] 1 . Here we would like to mention only that the definition, which will be used in the following, unambiguously follows from the requirement that the elastic amplitude vanishes in the limit ν → ∞, for fixed q 2 , and thus obeys an unsubtracted dispersion relation in the variable ν. Under this requirement, the elastic part is given by: where G E , G M denote the electric and magnetic (Sachs) form factors of the nucleon.
The inelastic invariant amplitudes are defined as T inel The amplitudes T inel i obey dispersion relations in the variable ν: Here, one has already taken into account the fact that, according to Regge theory, the dispersion relations for T inel 1 , T inel 2 require one subtraction and no subtractions, respectively. The lower integration limit is equal to ν th = (W 2 th − m 2 − q 2 )/(2m), with W th = m + M π , where M π is the pion mass. The quantities V 1 , V 2 denote the absorptive parts of T inel 1 , T inel 2 . They can be expressed through the experimentally observed total (transverse, longitudinal) electroproduction cross sections σ T (ν, q 2 ), σ L (ν, q 2 ).
The choice of the subtraction point ν 0 is arbitrary. In the literature, the choice ν 0 = 0 is often used. The quantity S inel 1 (q 2 ) = T inel 1 (0, q 2 ) is usually referred to as the subtraction function. Analogously, one can define the full subtraction function that includes the elastic part as well: . At q 2 = 0 the inelastic part of the subtraction function is given by where κ, β M denote the anomalous magnetic moment and the magnetic polarizability of the nucleon, respectively, and α 1/137 is the electromagnetic fine-structure constant.
Recently, a different subtraction function was introduced in Refs. [23,24] . The subtraction point has been chosen at ν 0 = iQ/2, where Q = −q 2 . The new subtraction function is expressed through the amplitudeT : At Q 2 = 0 one has:S The two subtraction functions are closely related to each other. Namely, the difference S inel 1 (q 2 ) −S(q 2 ) is given through a convergent integral over the experimentally measured electroproduction cross sections. Hence, it suffices to calculate one of these subtraction functions. Since the choice ν 0 = 0, in contrast to ν 0 = iQ/2, can be implemented on the lattice in a straightforward manner [8,9], we stick to this choice.
B. Extraction of the subtraction function on the lattice Below, we shall collect all formulae which are needed for the extraction of the subtraction function on the lattice with the use of the background field method. More details are contained in the original papers [8,9]. Here, we consider the nucleon placed in a periodic magnetic field on the lattice with a spatial size L (the temporal size of the lattice is assumed to be much larger and is effectively set to infinity). The configuration of the magnetic field is chosen as: where e denotes the proton charge. Because of the periodic boundary conditions, the available values of ω are quantized: 2 The energy levels of a nucleon in the magnetic field depend on the projection of the nucleon spin along the z-axis. In Ref. [9] it has been shown that the spin-averaged level shift in the magnetic field with a given configuration is given by: Note that here T 11 L (p, q) denotes the 11-component of the full Compton scattering amplitude in a finite volume (in other words, T 11 L (p, q) includes both inelastic and elastic parts). Further, q 2 = −ω 2 . Hence, placing a nucleon in the periodic magnetic field enables one to extract the amplitude at non-zero (albeit discrete) values of q 2 < 0. The other variable is ν = p · q/m = 0 in the given kinematics. Thus, in order to obtain a non-zero values of ν, one has to put the nucleon in a moving frame.
Note also that due to the lack of Lorentz invariance on a finite hypercubic lattice, the decomposition of this amplitude into invariant amplitudes in a form given in Eq. (4) does in general not hold. However, all quantities in Eq. (12) are well-defined in a finite volume. For example, in perturbation theory, T 11 (p, q) is given by a sum of all diagrams at a given order, where all momentum integrations are replaced by finite-volume momentum sums. In the infinite-volume limit one has: The finite-volume corrections in the above formula are suppressed by a factor exp(−M π L), multiplied by powers of L. As already mentioned in the introduction, despite the exponential factor, the corrections can still be sizable for the present-day lattices. Last but not least, let us stress once more that, for all values of L, Eq. (12) enables one to extract a perfectly well-defined quantity T 11 L (p, q), which in the infinite volume limit yields the quantity S 1 (q 2 ) that we are after. This demonstrates explicitly that in this setup one could avoid any ambiguous interpretation of the results as mentioned in Ref. [15].
We conclude this section by briefly specifying the scope and aims of the present paper. It is clear that the extraction of the Compton amplitude on the lattice can be carried out only in a restricted kinematical domain. For instance, if the variable ν lies above the pion production threshold ν th = (2m) −1 (W 2 th − m 2 − q 2 ), then the extracted matrix element does not possess an infinite-volume limit. This can be seen immediately since, in the infinite-volume limit, the corresponding amplitude is complex, whereas the amplitude extracted from an Euclidean lattice is always real. Hence, in order to arrive at the infinite-volume amplitude, one has either to take into account the proper Lellouch-Lüscher factor in analogy with Ref. [31], see also Refs. [32][33][34][35][36][37][38] (a first step in this direction has been made in Ref. [39]), or to use an approach that resembles the optical potential method of Ref. [40], see also Refs. [41,42]. All this is very complicated and not even needed to achieve the goals we have stated in the beginning. Indeed, given the subtraction function, which is obtained from the Compton amplitude at ν = 0, one may restore the whole Compton amplitude by using dispersion relations. The whole uncertainty related to the fixed poles then resides in the subtraction function, and the rest is uniquely determined by analyticity, unitarity and the experimental input.

C. Effective Lagrangian
In this paper the forward Compton scattering amplitude will be calculated in BChPT at order p 4 , both in the infinite and in a finite volume. Below we specify the effective Lagrangian with the pions and nucleons, which is needed for such a calculation. The leading-order effective Lagrangian of pions interacting with external sources has the form [43]: where χ = 2B(s + ip), D µ U = ∂ µ U − ir µ U + iU l µ and the 2 × 2 matrix U represents the pion field. The parameter B is related to the quark condensate, F is the pion decay constant in the two-flavor chiral limit, and s, p, l µ = v µ − a µ and r µ = v µ + a µ are external sources. The notation · · · denotes the trace in flavor space.
The full order four effective Lagrangian of nucleons interacting with pions and external fields is given in Ref. [44]. Below we specify only those terms, which contribute in our calculations: 3 where In the above expressions, Ψ denotes the nucleon field, u = √ U contains pion fields, v (s) µ is the part of the vector current that is proportional to the unit matrix in the flavor space, It is assumed that the above Lagrangian is used to generate Feynman diagrams, which will be evaluated by using the EOMS renormalization scheme. Acting in this manner, there is no need to explicitly display the counterterms of the effective Lagrangian that remove the power-counting breaking contributions from the loop diagrams. Hence, the numerical values of the finite parts of the LECs correspond to the EOMS scheme.
To obtain the expressions that correspond to diagrams with external photons, we need to substitute s = M, p = a µ = 0, v µ = −eτ 3 A µ /2 and v (s) µ = −eA µ /2, where A µ is the electromagnetic field. We work in the isospin limit m u = m d =m, and M 2 = 2Bm is the pion mass at leading order.

A. The workflow
We calculate doubly-virtual Compton scattering on the proton and on the neutron separately up-to-and-including O(p 4 ). There are tree and one-loop diagrams contributing to this process at the given accuracy. Standard powercounting rules of low-energy chiral effective field theory apply to these diagrams [46,47]. More specifically, we are assigning chiral order −2 to pion propagators, nucleon propagators count as of order −1, the interaction vertices originating from the effective Lagrangian of the order N count also as of order N , and the integrations over loop momenta are assigned order 4. While the power counting rules are directly applicable to the tree diagrams, the loop diagrams of our manifestly Lorentz-invariant formalism contain pieces that violate the counting rules. However, powercounting violating terms are polynomials in the quark masses and external momenta and thus can be systematically absorbed in the redefinition of the parameters of the effective Lagrangian. In our calculations, we use dimensional regularization supplemented with the EOMS scheme [48,49]. In this scheme, the polynomials in quark masses and external momenta which break the power counting up to a given chiral order, are dropped from each diagram. This naturally guarantees that the renormalized one-loop diagrams satisfy power counting 4 . Note that, in difference to the results of the heavy baryon formalism [52,53], our expressions for loop diagrams contain an infinite number of higher-order terms.
In the calculations of the diagrams we used the programs FeynCalc [54,55] and X-package [56]. We have verified that the sum of all tree and one-loop diagrams satisfy current conservation. This guarantees that the whole unpolarized amplitude is parameterized in terms of two invariant functions T 1 and T 2 . However, this does not apply to the individual diagrams. In order to extract the contributions of separate diagrams to T 1 and T 2 , we single out the contributions to the coefficients of the structures q µ q ν and −q 2 p µ p ν /m 2 , since these appear exclusively in K µν 1 and K µν 2 , respectively, see Eq. (4). The individual contributions, which are listed below, should be interpreted in this sense. Finally, using the LSZ scheme, we add these contributions and multiply the result with the residue of the nucleon propagator at the pole corresponding to one-nucleon state. Acting in this manner, one gets the full expressions of T 1 (ν, q 2 ) and T 2 (ν, q 2 ) we are interested in.
Next, we need an algorithm for the separation of the elastic contributions from T 1 and T 2 . For instance, let us first extract the s-channel pole. To this end, we multiply the full expressions for T 1 (ν, q 2 ), T 2 (ν, q 2 ) by 2mν + q 2 and then substitute ν = −q 2 /(2m). Apparently, as a result of this procedure, one obtains the residue of the s-channel pole. In order to determine the residue of the u-channel pole, we multiply the amplitudes by 2mν − q 2 . Finally, adding both pole terms together, we arrive at the elastic amplitudes T el 1 (ν, q 2 ) and T el 2 (ν, q 2 ). These vanish in the limit when ν → ∞ and q 2 stays fixed. This exactly coincides with our definition of the elastic amplitudes.

B. The amplitude in the infinite volume
The invariant amplitudes T 1 , T 2 up-to-and-including O(p 4 ) are given as sums over the contributions of the diagrams shown in Figs. 1-4: The Here, we remind the reader that, under the individual contributions to the invariant amplitudes T 1 , T 2 , we understand the scalar factors that multiply the structures q µ q ν and −q 2 p µ p ν /m 2 , respectively. Further, Z N is the residue of the nucleon propagator at the pole: which counts as Z N = 1 + O(p 2 ). For this reason, one has to take this factor into account only together with the tree-level diagrams at O(p) and O(p 2 ). Note also that, to this accuracy, one may replace m, F by their physical values g A , m, F π . The loop functions which enter the above expression are tabulated in Appendix A. The derivative in the function B 0 (denoted by the prime) is taken with respect to the first argument.
In the Appendices B and C we list the individual contributions to T 1 (ν, q 2 ) and T 2 (ν, q 2 ) for ν = 0. The full expressions for a generic ν are much more complicated and we do not display them here explicitly. These can be extracted from the Mathematica notebook, which is posted on the archive as an ancillary file. Note also that the expressions, given in these appendices, should be understood as 2 × 2 matrices in the isospin space, folded by the isospin wave functions of a proton or a neutron (not shown separately). For example, the factor 1 + τ 3 is equal to 2 and 0 for the proton and the neutron, respectively.

C. Numerical input
In order to make numerical predictions, one has to fix the values of the low-energy constants (LECs) that enter the amplitudes. It should be noted that, albeit these LECs do not depend on the quark masses by definition, such a dependence sneaks in when these are determined through the fit of the amplitudes calculated in a given chiral order to the experimental data. In certain cases, such remnant quark mass dependence can be numerically significant and, in addition, lead to the different results when different schemes, say, the Infrared Regularization (IR) or the EOMS scheme, are used. This fact should be kept in mind once input from different fits is used in the amplitudes.
The LECs that appear in the amplitudes fall into different groups. We shall use g A = 1.2672 and F π = 92.3 MeV throughout the paper, and M π , m will be identified with the charged pion and the proton masses, respectively. The order p 2 LECs c i are studied in most detail and rather precise values for these are available in the literature. Moreover, different fits (see, e.g., Refs. [57][58][59]) yield results which are compatible with each other at an accuracy that is sufficient for our purposes. The recent and very precise determination of c 1,2,3,4 from πN input has been performed in Ref. [60] using the matching of the chiral representation to the solution of the Roy-Steiner equations: In the following calculations, we shall use these values. Owing to the fact that the quoted uncertainties are so small, one may neglect their impact on the total uncertainty and safely stick to the central values.
Turning to the LECs c 6,7 , we note that these appear in the amplitudes in combination with the O(p 4 ) LECs: where M 2 = M 2 π at this order. These are exactly the contributions which appear in the anomalous magnetic moments of the proton and the neutron, κ p = 1.793 , κ n = −1.913, and can be fitted to the latter. A consistent extraction of these couplings has been performed in Ref. [61], which gives the results for the IR and EOMS schemes separately (no errors are attached to their results). Here, we quote the result for the EOMS scheme only, as this is used in our calculation:c Next, the O(p 3 ) LECs d 6,7 and O(p 4 ) LECs e 54,74 enter the expression of the electric and magnetic radii of the proton and the neutron, and can be fitted by using experimental data on the nucleon electromagnetic form factors. This was done in Ref. [61], where the values of these LECs are given, again separately for different schemes and with no uncertainties attached: and Note that the values for d 6 , d 7 are consistent with the earlier determination in Ref. [62].
Let us now turn to the last group of the O(p 4 ) LECs, which are related to the nucleon polarizabilities. At order p 3 , the polarizabilities are predictions free of LECs [63]. The O(p 4 ) LECs that contribute to the amplitudes at q 2 = 0 can be related to the polarizabilities, and we use (experimental or lattice) input for the latter. This allows to determine four linearly independent combinations of the O(p 4 ) LECs which, according to Ref. [45], can be defined as: e ± x = 2e 90 + e 94 + e 117 ± e 92 , e ± y = 2e 89 + e 93 + e 118 ± e 91 .
In Ref. [45] the results of three different fits for e + x , e + y are presented. In the following, however, we shall not use these results.
In order to carry out the comparison with the results from the literature, we perform a numerical evaluation of the subtraction functions S inel 1 (q 2 ) andS(q 2 ), which can be written as: In other words, in order to minimize the uncertainty, we aim at a description of the q 2 -dependence of the subtraction functions only, their values at q 2 = 0 are considered as input. Stated differently, up-to-and-including order p 4 , the q 2 -dependence of the subtraction functions S 1 (q 2 ) andS(q 2 ) (but not their normalization at q 2 = 0) is determined by the LECs, which are rather well known from the fit to the data on the low-energy πN scattering and nucleon electromagnetic form factors. Thus, the q 2 -dependence can be determined very accurately from BChPT.
The experimental values for the electric and magnetic polarizabilities are summarized in the recent paper by Melendez et al., [64] (see e.g. Refs. [65][66][67][68] for some earlier work): For the difference proton-neutron one gets: All quantities are given in units of 10 −4 fm 3 .
In Ref. [20], using Reggeon dominance, the isovector electric and magnetic polarizabilities have been predicted with an accuracy that supersedes the experimental precision. For instance, the value for the electric polarizability, extracted from the recent Review of Particle Physics [69], is given by α p−n E = −0.6(1.2). On the other hand, using Reggeon dominance and the experimental value for α p−n E + β p−n M from Ref. [64], which was determined by using the Baldin sum rule, one gets: Finally, recently a very accurate lattice calculation of the magnetic polarizability has become available [70]: One can combine this with the experimental result for α p−n E + β p−n M from Eq. (27) in order to get a more accurate estimate: To summarize, we have now three sets of polarizabilities, referred to as model A [Eq. (27)], model B [Eq. (28)] and model C [Eq. (30)]. These correspond to the purely experimental input, the Reggeon dominance hypothesis and the combination of the lattice results with experimental data. Below, we shall evaluate the difference of the subtraction functions for the proton and the neutron 5 using the input from the three distinct models:

D. The subtraction function
The main goal of the present paper is to evaluate the finite-volume corrections to the Compton amplitude. However, having the expression of the infinite-volume amplitude at hand, one may compare it to the known results from the literature. For instance, here we shall discuss the comparison to the subtraction functions S inel 1 (q 2 ) andS(q 2 ) for proton minus neutron, obtained from the experimental input by using the Reggeon dominance hypothesis in Refs. [20] and [24], respectively. Note that the experimental input, used in these papers, leads to a very large uncertainty at small values of Q 2 , which comes mainly from the resonance region above the ∆-resonance. On the other hand, the results of Chiral Perturbation Theory become generally unreliable at higher values of Q 2 . Thus, combining both calculations, one can get a coherent picture of the Q 2 -dependence of the subtraction function in a wide interval and check the consistency of the Reggeon dominance hypothesis: Should it turn out that there is an apparent mismatch between the low-Q 2 and high-Q 2 regions, this might cast doubt on the above hypothesis.
In order to reduce the error, which stems from the poor knowledge of the higher order LECs, in our calculations we have attempted to evaluate the Q 2 -dependence of the subtraction functions by subtracting their values at the origin. Eventually, the latter quantity is expressed in terms of the electric and magnetic polarizabilities, see Eq. (25). The polarizabilites can be fixed from the different inputs, leading to what we term the models A,B and C, see Eq. (31). The corresponding results in the interval 0 < Q 2 < 0.3 GeV 2 are displayed in Fig. 5 for the function S inel 1 (q 2 ) and in Fig. 6 for the functionS(q 2 ). Note that the uncertainty, which is displayed here, comes entirely from the poor knowledge of the polarizabilities. In other words, we assume that the uncertainties coming from other LECs are much smaller and do not contribute significantly to the error (it is clear that even attaching a reasonable uncertainty to other LEcs and adding uncertainties in quadrature, the changes will barely be visible due to the huge uncertainty in the polarizabilities). Note also that, in this scheme, the subtraction functions at O(p 3 ) are no more considered as a parameter-free prediction, but contain polarizabilities as input.
From Figs. 5 and 6 one may conclude that the results obtained in BChPT and by using the Reggeon dominance are reasonably consistent with each other within the error bars. In case of the subtraction functionS(q 2 ) this can be seen more clearly, since the results of the calculations, based on the Reggeon dominance, are available for all values of Q 2 down to Q 2 = 0, see Ref. [24] (in case of S inel 1 (q 2 ), the calculations extend down to Q 2 = 0.5 GeV 2 , see Ref. [20]). Note, however, that the uncertainty in different models, considered in the present paper, lead to very large error bars in these plots. What is more important in our opinion, is that the amplitudes calculated in the effective field theory show a smooth behavior in the vicinity of Q 2 = 0, stated differently, no rapid variations are observed. Note also that the convergence is quite poor, and the picture changes significantly when going from O(p 3 ) to O(p 4 ). Still, we do not observe any apparent disagreement to the Reggeon dominance. Also, it should be noted that both the O(p 3 ) and O(p 4 ) contributions are part of the complete one-loop amplitude, so a true test of convergence could only be achieved by going to the two-loop level. This, however, is beyond the scope of this paper.
We finish this section by briefly mentioning related calculations in the literature. Our result for the polarizabilities fully agrees with that of Ref. [71], where the calculations were done in the relativistic BChPT at O(p 3 ). Next, our subtraction function S inel 1 (q 2 ) coincides with the one from Ref. [29] at order p 3 . Moreover, in the recent paper [72], these calculations were extended up-to-and-including order p 4 . However, at O(p 4 ), the contribution from the ∆resonance comes into play and this renders the direct comparison more complicated (note, however, that for proton minus neutron the leading contribution of the ∆ drops out, as already noted in Ref. [73]). For this reason, in this paper we restricted ourselves to O(p 3 ) and verified that the sum of the polarizabilities α E + β M is algebraically reproduced in our calculations, both for the proton and the neutron. Further, expanding our result in inverse powers of m, one should reproduce the non-analytic pieces of the Heavy Baryon ChPT (HBChPT). The result of Ref. [74], however, cannot be obtained in this way, whereas the same quantity, evaluated at O(p 4 ) within HBChPT in Ref. [75], can be readily reproduced after the expansion. Finally, note the calculation of the imaginary part of the subtraction function, carried out in Ref. [76] in the framework of HBChPT (the dispersion relations are used to recover the real part). The result in that paper is given in a form that is not very suitable for a straightforward comparison.
Another group of the papers deals with the calculation of the subtraction function by using dispersion relations and experimental input (see, e.g., Ref. [77]), or modeling it, taking into account the constraints at Q 2 = 0 and at Q 2 tending to infinity [19,21,22]. This work has been discussed in great detail in Refs. [20,24], to which the interested reader is referred to.

IV. FINITE VOLUME CORRECTIONS
A. Analytic expression for the finite-volume amplitude The diagrams that contribute to the Compton amplitude are the same in the infinite and in a finite volume. The only difference consists in replacing the three-dimensional integrals with the sums over discrete lattice momenta (we take that the effects related to a finite size of a lattice in the temporal direction are already taken into account during the measurement of the energy levels). Assuming periodic boundary conditions, in the loop integrals one has to replace: Some of the loop integrals in the infinite volume diverge. In a finite volume, this divergence can be dealt with using dimensional regularization in the same manner as in the infinite volume. The counterterms that remove the divergences are the same in both cases. In the following, we shall use this fact and write down the finite-volume sums in a dimensionally regularized fashion, without specifying how this is done. These sums will be further split into the infinite-volume parts and the corrections. The regularization is relevant only for the first parts, where the standard prescription can be applied. The finite-volume corrections are ultraviolet-finite and can be calculated in four dimensions.
A closely related question concerns the presence of the power-counting breaking terms in the covariant BChPT. It must be stressed that there are no such terms in the finite-volume part of the amplitude, as all such terms are suppressed by a factor exp(−mL) containing the nucleon mass. On the contrary, in the infinite-volume part these terms are present and can be dealt with, e.g., by using the EOMS prescription, as described above. The reason for this is, of course, that the breaking of the power-counting rules is a high-energy phenomenon that emerges for loop momenta of the order of the nucleon mass. On the other hand, the large-L behavior is governed by the momenta of the order of the pion mass. This region does not contribute to the breaking of the power counting.
In order to carry out the calculations, one needs the expression of the nucleon Z-factor in a finite volume. This quantity is defined in the rest-frame of the nucleon, where the nucleon propagator is given by: In the above expression, the integration is carried out in a finite box. Further, where Σ L (p) denotes the self-energy of the nucleon in a finite volume, and The finite-volume mass of the nucleon is implicitly given through the solution of the equation that contains the scalar functions A L (p 0 ), B L (p 0 ): This equation can be solved by iteration, expressing m L order by order through the infinite-volume parameters. Further, the Z-factor in a finite volume is given by the residue of the nucleon propagator at the pole p 0 = m L . It can be also expressed in terms of the functions A L (p 0 ), B L (p 0 ) and the derivatives thereof with respect to the variable p 0 : The explicit expression for this quantity is given by: The derivatives of the loop integrals have been reduced again to loop integrals utilizing the algebraic identities that are specified in Refs. [45,78].
Finally, note that Lorentz symmetry is broken in a cubic box and one can no more use the decomposition of the Compton amplitude into two scalar amplitudes. This is, in fact, not needed, because the energy shift of a nucleon in the periodic field is directly given by the 11-component of the Compton tensor [8,9]. Putting together all contributions, one can write: The individual contributions are given by: O(p): protonT O(p 2 ): neutronT O(p 2 ): protonT  (2c 6 + c 7 )(8p α p β + 9p α q β + 4q α q β )B αβ (1,1) (M, m; p + q) Note that in the above expressions, we have replaced m L , emerging from the kinematics, by the infinite-volume nucleon mass, m. Up to the chiral order we are working, this is a perfectly valid procedure. Finally, the expressions for the various finite-volume sums, which enter the above formulae, are listed in appendix D.

B. Numerical results
In this section, we evaluate the finite-volume corrections to the 11-component of the Compton tensor, which enters the expression of the energy shift in the periodic background field. The quantity, which will be calculated here, is given by: We calculate this quantity for the physical value of the pion mass and for several different values of q 2 , separately for the proton and the neutron. Having explicit expressions for the amplitude, it is straightforward to carry out the calculations for unphysical quark masses as well, if needed. It should be stressed that we are mainly interested in the order-of-magnitude estimate of the correction, which is needed to answer the following question: How large the box size L should be so that one can safely neglect the finite-volume artifacts?
The results at O(p 3 ) can be obtained directly from Eqs.
This choice covers model A, as well as model C.
The finite-volume corrections to the Compton amplitude are shown in Figs. 7 and 8 for the following values of the variable Q 2 : These figures contain our main result, answering the question about the feasibility of the extraction of the subtraction function on the lattice. It is seen that, both for the proton and the neutron, the finite-volume corrections are encouragingly small for M π L ≥ 4 (the fact that they are slightly smaller for the proton than for the neutron stems from the presence of the large second-order pole term proportional to the combination 2c 6 + c 7 , in the infinite-volume proton amplitude). Further, it is very comforting to see that the convergence of the result at fourth order is reasonable. Moreover, the uncertainty caused by the poor knowledge of the O(p 4 ) LECs is indeed moderate in the final results (for example, in case of the proton, it is hardly visible by the bare eye). Taking this fact into account, one might wonder whether the uncertainty in the LECs at lower orders might play a more significant role. Following our expectation, a 20-30 % error in a final result generously covers the effect coming from these LECs. Another easy way to estimate the uncertainty of calculations (not limited necessarily to the poorly determined LECs) is to compare the results at O(p 3 ) and O(p 4 ). Taking into account the fact that the present study was primarily intended to serve as a rough estimate of the size of the exponentially suppressed corrections to the amplitude, we did not try to investigate this question further.
Finally, note that the relative correction stays almost constant from Q 2 0 to Q 2 2M 2 π and, possibly, even for higher values of Q 2 . In other words, the finite-volume artifacts do not hinder an accurate extraction of the amplitude at large Q 2 . Here we remind a reader that an accurate measurement of the inelastic part on the lattice becomes more difficult as Q 2 → 0, because the elastic contribution dominates in this limit [8,9]. Hence, the finite-volume artifacts do not further restrict an interval in Q 2 , where an accurate calculation of the subtraction function is possible.

V. SUMMARY
i) Using Baryon Chiral Perturbation Theory and the EOMS renormalization scheme of Refs. [48,49], which guarantees that the renormalized expressions satisfy the standard power counting, we have evaluated the doubly virtual spin-averaged Compton scattering amplitude off nucleons up to O(p 4 ) in the infinite volume. We have further calculated finite-volume corrections to the so-called subtraction function at the same chiral order.
ii) The inelastic parts of the infinite-volume subtraction functions S inel 1 (q 2 ) andS(q 2 ) show a rather regular behavior at small values of q 2 , that is of the order of the pion mass squared. None of the loop diagrams up-to-and-including order p 4 leads to a rapid variation of the calculated subtraction functions at small scales. Note that in these calculations we have fixed the numerical values of some of the O(p 4 ) LECs through the proton and neutron magnetic polarizabilities which, at present, are not known to high accuracy (especially, the one for the neutron).
iii) The main result of this work is the calculation of the finite-volume corrections in the Compton scattering amplitude. These calculations are interesting, first and foremost, in view of the perspective of a measurement of the subtraction function on the lattice using periodic external fields. Note also that, since the cubic lattice does not preserve Lorentz invariance, the definition of the subtraction function in a finite volume is ambiguous. On the other hand, what is extracted from the nucleon mass shift, measured on the lattice, is a particular component of the second-rank Compton tensor T 11 L , which is a well-defined quantity and which tends to S 1 (q 2 ) (modulo an overall kinematic factor) in the infinite-volume limit. The numerical results quoted in this work refer to this quantity. iv) At this stage, we do not know, on the one hand, how the other subtraction functionS(q 2 ) can be measured on the lattice. On the other hand, the two subtraction functions are related to each other: Their difference is a convergent integral containing experimentally measured electroproduction cross sections.
v) Our results show that the exponentially vanishing finite-volume corrections to the quantity T 11 L amount up to 2-3% percent or less at M π L 4 both for the proton and the neutron. This means that one can extract the infinite-volume subtraction function with a good accuracy already using reasonably large lattices. We also note that the convergence of our results is rather good, and the poor knowledge of the O(p 4 ) LECs does not pose a real obstacle as the resulting uncertainty is very small.
vi) As pointed out in Refs. [8,9], the large elastic part renders the extraction of the subtraction function at low  The one-loop integrals appearing in our calculations are defined as follows: One factor in the denominator: Two factors in the denominator: = (g αν g βµ + g αµ g βν + g αβ g µµ )B (2,1) 0000 (q 2 ; M, M ) + q µ q ν q α q β B (2,1) 1111 (q 2 ; M, M ) + (q β q ν g αµ + q µ q ν g αβ + q µ q β g αν + q α q β g µν + q µ q α g βν + q ν q α g µβ )B (2,1) Three factors in the denominator: O(p) contributions: O(p 2 ) contributions: Below one finds the one-loop contributions to the amplitudes T 1 (0, q 2 ) and T 2 (0, q 2 ): O(p 3 ) contributions: Contributions of tree diagrams to the elastic parts of T 1 (0, q 2 ) and T 2 (0, q 2 ): Here, V denotes an integral in the finite volume, which really is a sum. The calculation of these sums by using the Poisson formula is considered in appendix E.
where n µ = (0, n) is a unit spacelike vector, whose components take integer values. Further, K ν (z) denote the modified Bessel functions of the second kind.
The finite-volume sums, which are displayed in appendix D, contain an infinite-volume piece and finite-volume correction. The ultraviolet divergences are contained only in the former, while the latter is ultraviolet-convergent and vanishes exponentially for large values of L. In order to ease the notation, we list only the finite-volume corrections. The following notation is used: The full list of the finite-volume sums entering the amplitude at O(p 4 ) is given below. Note that in the expressions, which contain only nucleon propagators, the finite-volume corrections are extremely small (proportional to the factor exp(−mL)) and can therefore be neglected. We shall indicate these quantities by writing ≈ 0 at the end. Also, note that the structure of the integrands, which appear in the infinite and in a finite volume, is in general different. This is related to the fact that Lorentz-invariance is used in the infinite volume to reduce tensor integrals to scalar ones. Some factors in the denominator get canceled during this procedure. One cannot apply the same trick in a finite volume.