Nucleon in a periodic magnetic field: Finite-volume aspects

The paper presents an extension and a refinement of our previous work on the extraction of the doubly virtual forward Compton scattering amplitude on the lattice by using the background field technique, Phys. Rev. D 95, 031502 (2017) (arXiv:1610.05545). The zero frequency limit for the periodic background field is discussed, in which the well-known result is reproduced. Further, an upper limit for the magnitude of the external field is established for which the perturbative treatment is still possible. Finally, the framework is set for the evaluation of the finite-volume corrections allowing for the analysis of upcoming lattice results.


I. INTRODUCTION
Using the background field technique on the lattice for the extraction of various hadronic observables has proven to be extremely efficient. As examples we mention the measurement of magnetic moments, polarizabilities and axial-vector matrix elements of baryons and light nuclei in constant background fields [2][3][4][5][6][7]. Moreover, non-uniform background fields have been used for the calculation of the hadronic vacuum polarization tensor, hadronic formfactors and the nucleon structure functions, as the fields, which are periodic in space, allow one to measure current matrix elements at a given non-zero three-momentum transfer [8][9][10].
Different scenarios for implementing periodic background fields in lattice QCD calculations are considered in Ref. [11].
In Ref. [1] we described a framework, based on the background field method, which enables one to extract the doubly-virtual forward Compton scattering amplitude from lattice QCD calculations (note that later, in Ref. [10], a very similar formula was given without a derivation, see Eq. (12) in that paper). Low-energy Compton scattering plays an indispensable role in probing the electromagnetic structure of hadrons (for recent work, see Ref. [12]).
For example, it enters the expression for the proton-neutron electromagnetic mass difference [13], as well as the expression for the Lamb shift of the muonic hydrogen, which is used to extract the value of the proton radius (see, e.g. Ref. [14]). In this paper, we in particular focus on the relevant spin-independent invariant amplitudes, denoted as T 1 and T 2 , respectively. The experimental data on the structure functions completely determine the amplitude T 2 through dispersion relations. The fixed-q 2 dispersion relation for the amplitude T 1 , however, requires a subtraction. Thus, the subtraction function S 1 (q 2 ) ≡ T 1 (0, q 2 ) remains the only input in the calculations, which is not fixed by experimental data. Its elastic part is essentially given by the Born terms, but the inelastic piece is known only at the real photon point q 2 = 0 (the low-energy theorem, see, e.g., Ref. [15]).
In the past, there have been attempts to model the subtraction function, using phenomenological parameterizations [16][17][18]. However, this type of approach inherently contains a systematic error, which is very hard to control. Further, the subtraction function can be extracted from the Compton scattering amplitude, calculated in the low-energy EFT of QCD [19][20][21]. However, here the difficult question about the convergence of the chiral expansion arises, namely, up to which value of q 2 the results of the chiral expansion can be trusted. Recently, the authors of Ref. [15] have been able to determine the subtraction function by invoking the so-called Reggeon dominance hypothesis, considered first in Ref. [22]. In particular, it is assumed that the forward Compton scattering amplitude does not contain any fixed pole. In Regge theory, such a pole generates an energy-independent contribution to the amplitude (such as, e.g., local two-photon couplings in scalar QED). If the fixed poles are present, the subtraction function, in general, deviates from the predicted one. In some cases, e.g., the q 2 -independent fixed pole [23], the behavior of the S 1 (q 2 ) can be also predicted (see Ref. [24]), and is different from the one calculated in the absence of the fixed pole.
Hence, lattice QCD provides a model-independent approach to the verification of the Reggeon dominance hypothesis. The question whether there is a fixed pole in the Compton scattering is of conceptual interest and is still open.
In Ref. [1], we considered the case of a nucleon placed in a static periodic magnetic field B = (0, 0, B 3 ) with B 3 = −B cos(ωx) and ω = (0, ω, 0). It has been shown that, for ω ≠ 0, the measurement of the spin-averaged energy shift of the nucleon in this field allows one to extract the value of the subtraction function S 1 (−ω 2 ) at non-zero values of q 2 = −ω 2 . The relation between these two quantities takes the form: Here, δE denotes the spin-averaged energy shift, and m stands for the nucleon mass.
The result, obtained in Ref. [1], still leaves room for improvement. In particular, one should find the answer to the following questions: i) In the limit ω → 0, we arrive at the case of a constant magnetic field. This case is studied very well, both analytically and numerically. However, our expressions become singular in this limit, contradicting the expectations. One needs to understand how this limit can be approached smoothly.
ii) Our approach relies on a perturbative expansion of the energy shift δE in the external field strength B. What is the radius of the convergence of this expansion? Note that, for example, in the zero-frequency limit ω → 0, the radius is equal to zero in case of a charged particle, since the Landau levels are formed for any value of B. In Ref. [1], using heuristic arguments, for a given value of ω, we gave a very rough estimate of the maximal value of B, for which the perturbative expansion should still work. These arguments should be refined in order to obtain a more reliable result.
iii) Our expressions were obtained in the infinite-volume limit. The issue of the finitevolume corrections in the presence of the external fields is a rather subtle one, since gauge-invariant non-local objects (Wilson lines) can be formed in a finite volume. For this reason, it is mandatory to re-formulate the problem in a finite volume from the beginning and to give a consistent interpretation of the finite-volume result it terms of the subtraction function, defined in the infinite volume.
The aim of the present paper is to answer the questions given above. The plan of the paper is as follows. In section II we give a collection of basic definitions and discuss two different implementations of the external field on the lattice. In sections III and IV, we give two alternative derivations of the energy shift formula in the periodic external field, based on the matching to the non-relativistic EFT, as well as the direct derivation within the relativistic framework. Both settings are complementary to each other. For example, the first derivation is more intuitive and uses ordinary quantum-mechanical Rayleigh-Schrödinger perturbation theory for the energy levels. In particular, the zero-frequency limit ω → 0 as well as the issues related to the convergence of the perturbative expansion can be considered more easily in this formulation. By contrast, the relativistic formulation allows one to investigate the exponentially suppressed finite-volume corrections in a direct manner. For completeness, in the appendix we give yet another derivation of the expression for the energy shift, considering the behavior of the nucleon two-point function at large time separations.

A. Basic definitions
The matrix element of the electromagnetic current between one-nucleon states is given by Here, j µ (x) is the electromagnetic current, and q = p ′ − p, and p (p ′ ) and s (s ′ ) are the fourmomenta and spin projections of the initial (final) nucleon, respectively. Further, F 1 and F 2 denote the Dirac and Pauli form factors. The Sachs form factors are defined by The Dirac spinors are normalized asū(p, s ′ )u(p, s) = 2mδ s ′ s .
The Compton tensor is defined as: where q is the photon momentum. Taking into account Lorentz invariance, current as well as parity conservation, one arrives at the well-known decomposition of the matrix element in Eq. (4) in terms of Tarrach's amplitudes [25,26]. For our purposes, it is sufficient to consider the process in the forward direction p ′ = p and perform spin-averaging in Eq. (4): The tensor T µν (p, q) is related to the aforementioned invariant amplitudes T 1 , T 2 through the decomposition (see, e.g., Ref. [15]): where the kinematic structures K µν 1 , K µν 2 read Here, ν ≡ p ⋅ q m.
According to the asymptotic behavior of the structure functions at large values of the parameter ν, the dispersion relation for the amplitude T 1 (ν, q 2 ) requires one subtraction. It is usually performed at ν = 0, and hence the subtraction function S 1 is defined as The function S 1 (q 2 ) can be formally split in two parts: The elastic term S el 1 (q 2 ) is associated with the one-nucleon exchange in the s-and u-channels. The inelastic piece S inel 1 (q 2 ) is a regular function of q 2 . We use the same definition of the elastic part as in Ref. [1]: Little information is available on the inelastic part of the subtraction function S inel 1 (q 2 ). According to the low-energy theorem, its value at q 2 = 0 is given by Here, β M denotes the magnetic polarizability of the nucleon, α ≃ 1 137 is the fine structure constant and κ = F 2 (0) is the anomalous magnetic moment of the nucleon. At large values of q 2 , the asymptotic behavior of the subtraction function is fixed by the operator product expansion (see, e.g., Ref. [17,27]). Otherwise, it is unknown in the intermediate kinematic region 0 < −q 2 ≲ 2 GeV 2 , which is amenable to lattice simulations.

B. External field configuration
In Ref. [1], it was proposed to place the nucleon in the time-independent periodic magnetic field where B denotes the strength of the field and the frequency ω takes nonzero values. The components of the gauge field A µ (x) are chosen as follows: The magnetic flux is quantized in a finite box of size L: As discussed in Ref. [11], this quantization condition can be implemented on the lattice in two different ways. In the first scenario, the frequency ω is constrained and no constraint is imposed on the magnetic field strength B. In the second scenario, the situation is reversed.
Hence, we have: Only the first scenario was considered in Ref. [1]. In the present paper, we will exploit both quantization possibilities and demonstrate that the obtained results are quite different.
Obviously, the limit ω → 0 can be directly performed in the second setting only, where ω is a free parameter, unrelated to the box size L.

A. Method
The framework, which is based on the use of the non-relativistic EFT, consists of two steps. At the first stage, one matches the parameters of the non-relativistic Lagrangian to the expression of the relativistic two-point function of the nucleon in an external field. At the next step, one uses the resulting non-relativistic Hamiltonian to carry out the calculation of the spectrum. The advantage of the method is its transparency: the calculations of the spectrum are done by using ordinary perturbation theory in quantum mechanics. The setting is, however, not well suited for the calculation of the finite-volume corrections, which are proportional to exp(−M π L) where M π denotes the pion mass. Within this approach, these corrections should be included in the couplings of the non-relativistic Lagrangian through the matching procedure.
Let us consider the two-point function of the nucleon, placed in the external field A µ (x).
The path integral representation in Minkowski space reads: where the integration over all possible gluon, quark and antiquark field configurations is performed. Further, j µ (x) is the electromagnetic current, built of the quark fields, and Ψ(x) denotes the composite nucleon field operator in QCD. Expanding the right-hand-side of Eq. (17) up-to-and-including O(A 2 ), one obtains where the subscript "0" refers to the quantities evaluated in QCD without any external field, and we have used the fact that ⟨0 j µ (x) 0⟩ 0 = 0. Note that the expansion in Eq. (18) is written down for connected matrix elements (the subscript "conn" is omitted everywhere for brevity).
In the above quantities, the nucleons are in general off the mass shell. Performing the Fourier transform in Eq. (18), amputating the external nucleon legs, and putting the external nucleons on the mass shell, we see that the nucleon electromagnetic vertex ⟨p ′ , s ′ j µ (0) p, s⟩ emerges at order A. At order A 2 , the Compton tensor, defined in Eq. (4), is obtained from the matrix element ⟨0 T Ψ(x)Ψ(y)j µ (z 1 )j ν (z 2 ) 0⟩ 0 . Note that, in general, the described procedure is equivalent to replacing the nucleon fields Ψ(x) andΨ(y) by the out-and ingoing nucleon states, ⟨p ′ , s ′ and p, s⟩ respectively. This fixes the overall normalization of the quantity we are considering below.

B. Matching
The first few terms of the effective Lagrangian, which describe the interaction of the nucleon with an external electromagnetic field, are given by: where the ellipses denote higher order terms with derivatives, Here, ψ(x) denotes the two-component nucleon field, µ is the magnetic moment and α E , β M are the electric and magnetic polarizabilities, respectively. The coupling constant g takes the values g = 0, +1 for the neutron and proton, respectively, and E, B are the electric and magnetic fields (in order to simplify the notations, we include the factor e in the definition of the vector-potential A µ ). The NRQED Lagrangian at order m −4 is given, e.g., in Ref. [28].
If the properties of a nucleon at rest are calculated, it suffices to write down only the first few terms in the Lagrangian. However, we are studying the nucleon in a periodic field, with ω corresponding to the magnitude of the momentum transfer from the field to the nucleon.
Consequently, we have to retain all terms in the derivative expansion of the Lagrangian given by Eq. (19). In this case, the exact relativistic dispersion relation for the energy of the free nucleon with the three-momentum p, w(p) = m 2 + p 2 , is satisfied.
An important remark is in order. In the infinite volume, one is allowed to use partial integration in the Lagrangian. The same is true in a finite volume, if all fields (including the external electromagnetic field) are subject to periodic boundary conditions, since the surface terms vanish in this case. The above remark will be relevant, if the realization of the external field is carried out according to the scenario b) from section II B: the Lagrangians, which differ only by surface terms and lead to the same amplitudes in the matching condition, might yield a different spectrum in the finite volume. Bearing this in mind, we must for instance ensure that all terms of the Lagrangian that we write down are explicitly gaugeinvariant (not only up to the surface terms), otherwise, one is not guaranteed that the resulting finite-volume spectrum is gauge-invariant. Note also that we did not pay special attention to this issue in our previous paper [1], where it was anyway not relevant, since only the scenario a) was considered.
Taking the above issue into account, below we write down the explicit non-relativistic Lagrangian, describing the nucleon in an external field up-to-and-including O(A 2 ) in this field. Note also that we change the normalization of the fermion field by a factor (2W ) 1 2 with W = √ m 2 − D 2 as compared to Eq. (19), in order to ensure that the free one-nucleon states obey the relativistic normalization condition (see, e.g., Refs. [29,30]). The Lagrangian takes the following form: Here, and where Γ E B and Π EE EB BB denote the pertinent combinations of the effective couplings with the invariant tensors like g µν or ε µναβ and Pauli matrices for the spin (in the following, for brevity, we shall refer to Γ E B and Π EE EB BB merely as to the effective couplings). The Latin indices run from 1 to 3 (only space derivatives), whereas the Greek indices run from 0 to 3. The derivatives in the square brackets act only on the function within the brackets Also, as a convention, the values m, n, k = 0 correspond to no derivatives in Eqs. (23,24).
Expanding explicitly in powers of the external field A and using partial integration and the equations of motion, one may rewrite the above Lagrangian in a simpler form, already displayed in Ref. [1]. whereL whereas Here, the effective couplings Γ and Π are the linear combinations of the couplings appearing in Eqs. (22,23,24). This form of the effective Lagrangian is better suited for carrying out the matching to the relativistic amplitudes. For example, the couplings Γ should be matched to the current matrix element in Eq. (2). Calculating the same vertex function in the effective field theory with the LagrangianL 1 , we get ∞ m,n=0 This means that, expanding the nucleon form factor in a Taylor where T µν (p ′ , s ′ ; p, s; q) is the Compton tensor defined in Eq.
The second iteration of the LagrangianL 1 gives another term, M 2 , with where the tensor U µν (p ′ , s ′ ; p, s; q) is given by the sum of the nucleon pole terms: Note that, in order to derive the above expression, the matching at O(A) has been used.
An important remark is in order. The aim of the matching is to determine the couplings Γ and Π, which encode the physics at short distances. It can be carried out in the infinite volume, where no specific care about the quantization of the magnetic flux should be taken.
The latter will be, however, important in the calculation of the energy shift.

C. Perturbation theory for the energy levels
In the previous section, an effort was made to match the relativistic and non-relativistic theories at order A 2 . In this section, we shall be rewarded for this effort, using the resulting non-relativistic Lagrangian for the calculation of the energy spectrum of the nucleon in an external field. Also, up to this moment, we have not specified the external field. Here we assume that A µ is the static field described in section II B and the scenario a) is chosen.
Stationary energy levels exist in such background field configurations.
Consider the canonical Hamiltonian H, which is obtained from the non-relativistic Lagrangian L. In order to arrive at the non-relativistic normalization of states, used in quantum mechanics, it is convenient to rescale back the nucleon field, entering in this Hamiltonian, as ψ → (2w) −1 2 ψ. Further, we define the quantum-mechanical Hamiltonian H, which is given by the matrix element of H between the free one-nucleon states. H is a differential operator, which acts on the nucleon wave function: where Note that H is a 2 × 2 matrix in spin space.
The wave function of the nucleon in the external field obeys the Schrödinger equation: where the ψ n,s (x) denote stationary solutions in a finite volume, satisfying periodic boundary conditions. The eigenfunctions ψ n,s (x) and the eigenvalues w(k n ) of the unperturbed Hamiltonian H 0 satisfy the equation The unperturbed spectrum has the form and the normalized solutions are given by Here, we have introduced a double-bracket notation that is different from the relativistic case.
Namely, k n , s⟫ denotes the state vector in the non-relativistic theory, which corresponds to the unperturbed solution.
Next, we apply perturbation theory in order to calculate the shift of the ground state in the external field. By doing this, we implicitly assume that the structure of the spectrum is not changed by the background field which is sufficiently small. Below, we shall put this condition under scrutiny.
Let us start from the ground state energy shift at order A. The unperturbed spectrum is degenerate (the same energy for both spin projections), so the perturbation theory for the degenerate levels should be applied. As is well known (see, e.g., Ref. [31]), the first-order energy shift is the solution of the secular equation: where V s ′ s = ⟪0, s ′ H 1 0, s⟫.
The matrix element of the operator H 1 is given by Here,Ã µ (q) denotes the Fourier transform of the field A µ (x, 0) which, for the field configuration described in Eq. (13), gives The sum in Eq. (42) has precisely the same form as in the matching condition at O(A), Eq. (29). Accordingly, the matrix element of the operator H 1 takes the form: Setting p = p ′ = 0 in Eq. (45) and taking into account that ω ≠ 0, it is seen that the matrix elements V ss ′ vanish: Thus, there is no first-order correction to the energy shift: As expected, this result follows from the three-momentum conservation at the vertex of the three-point function. We again stress that it holds only for ω ≠ 0. Note also that, since the off-diagonal matrix elements vanish as well, the correct wave functions at this order are still given by Eq. (40).
The second-order contribution to the energy shift can be found again from the secular equation, which differs from Eq. (41) by the replacement The first term emerges from the second iteration of H 1 and another one is the matrix element of H 2 . The spin-averaged energy correction at O(B 2 ) consists of two pieces: where The first term is evaluated by using Eqs. (44,45). Taking into account the fact that w(0) = m, we obtain: wherep = (m, 0). Performing the summation over k n , we get: where the quantity F (ω) reads For the second piece we need to evaluate the matrix element of the operator H 2 : The integration leads to the expression where the integral I µ 1 ...µ l ,µν reads If q = 0, this integral has a non-zero value for µ 1 = ⋅ ⋅ ⋅ = µ l = 2, µ = ν = 1 and for even l, Inserting this expression into Eq. (56), one gets Further, using the matching condition at O(A 2 ) in Eq. (34), and setting p = 0, we see that the expression for the second energy correction δE ′′ s takes a compact form: whereq = (0, ω).
It remains to put all pieces together. Noting that the quantity U 11 (0, s; 0, s;q), defined by Eq. (33), is exactly equal to − 1 2 [F (ω) + F (−ω)], we finally arrive at the expression of the spin-averaged energy shift of the ground state, derived first in Ref. [1]: D. The zero-frequency limit The equation (61) does not posses a smooth zero-frequency limit. This is seen from the fact that, e.g., the quantity T 1 (0, −ω 2 ) includes the elastic contribution that diverges in this limit as 1 ω 2 . On the other hand, the result for the energy shift in the constant field is well known: it is finite and is proportional to the pole-subtracted part of the forward Compton amplitude. In this section, we shall discuss this apparent contradiction.
Let us start with the energy shift at O(A). As the frequency of the magnetic field tends to zero, the field approaches a constant value and the first-order correction to the energy shift does not vanish anymore. It is immediately seen that approaching smoothly the limit ω → 0 is not possible in the scenario a), where ω is quantized, according to Eq. (15) so that ω is either zero from the beginning or not. If ω ≠ 0, then the energy shift of the ground state, caused by the perturbation Hamiltonian in Eq. (45), is strictly zero. Consequently, one has to turn to scenario b). Here, one can immediately visualize the problem with the non-vanishing surface terms, which were mentioned above. For example, the matrix element of the current, entering Eq. (45), has the following representation: where and In the matching procedure, which is carried out in the infinite volume, we always have where q = (0, q, 0). Assuming now p ′ = p = 0 in this expression and then letting ω → 0 leads to the vanishing matrix element for any non-zero ω. On the other hand, using partial integration, one obtains ε kim q iÃk (q) = −iB m (q), where B is the magnetic field, and the matrix element of the Hamiltonian takes the form: This expression does not vanish anymore and, in the limit ω → 0, yields the well-known result for the first-order energy level splitting in the constant magnetic field. The reason for this inequivalence is immediately seen: in case of an external field, which does not obey periodic boundary conditions, the 3-momentum conservation is not guaranteed, and the equality q = ω does not hold anymore. At threshold, the vector q vanishes, but the vector potential contains the factor 1 ω and the result depends on the way the limit is performed.
It can be also seen that the above ambiguity disappears, if an explicitly gauge-invariant Lagrangian, defined in Eqs. (21,22,23) and (24), is used from the beginning. The equation (66), which leads to the correct result in the zero-frequency limit, is directly obtained by using the gauge-invariant Lagrangian without performing the partial integration.
The situation is pretty much the same in case of the second-order energy shift. Below, for simplicity, we shall consider the case of the neutron only. Further, we shall stick to the particular field configuration, described in section II B. The nucleon pole contribution to the energy shift (the analogue to Eq. (51)) in case of the external field that does not obey periodic boundary conditions, is given by where It is seen that, in case of the a periodic field, this factor does not reduce to the Kronecker delta-symbol, corresponding to the conservation of the total three-momentum. In the limit ω → 0 the above expression simplifies considerably, and we have where k ⊥ n denotes the components of the vector k n , perpendicular to the vector ω. The expression for the matrix elements, entering Eq. (68), can be read off from Eqs. (62) and (63). First of all, because the three-momentum, perpendicular to the direction of ω is conserved, only the terms that contain a 3 and a 4 , can potentially contribute. Further, comparing Eq. (65) and Eq. (66), it is clear that using the gauge-invariant Lagrangian (23) in our case boils down to the following heuristic prescription: write down the vertices in terms of two linearly independent vectors p ′ + p and q = p ′ − p, and replace everywhere q through ω as if the three-momentum was conserved. Now, one can ensure that the contributions from both terms, containing either a 3 or a 4 , vanish in the limit ω → 0. Indeed, as seen from Eq. (62), the term with a 3 contains q, which is eventually replaced by ω and the limit ω → 0 is performed afterwards. Further, using Eq. (63), one sees that a 4 is proportional to w(p ′ ) − w(p) or, equivalently, to q(p ′ + p). This expression also vanishes, when q is replaced by ω and the limit ω → 0 is performed (we remind the reader that the Fourier transform A 1 (q) stays finite for a non-zero q and ω → 0). Hence, the entire pole term does not contribute to the energy shift in the limit ω → 0, and the latter is given solely by the contact contribution: Only a single coupling Π 33 BB,ss contributes in this limit, since the derivative terms give vanishing contributions. Comparing the first term in the expansion of Eq. (24) to Eq. (19) and taking into account the different normalization of the nucleon field in these two Lagrangians, one immediately sees that Π 33 BB,ss is given by the magnetic polarizability and, thus, the standard formula for the spin-averaged energy shift δE = −2πβ M B 2 is reproduced in the limit ω → 0 (note that δE ′′ s , given by the above expression, does not depend on the spin orientation).
To summarize this part, we note that, in order to perform a smooth zero-frequency limit, one has to use the realization b) of the external field on the lattice, in which the frequency ω is not quantized. Using this realization for a finite ω, however, is not very convenient.
Apart from the subtleties, arising in the treatment of the surface terms, the final expression for the energy shift is rather complicated and simplifies only in the limit ω → 0. For this reason, in the following we stick to the scenario a).

E. Landau levels
Here we consider, how the Landau levels emerge from the periodic potential in the zero frequency limit ω → 0 . Further, we give an estimate for the maximum value of the field strength B, for which our method still works (note that a crude estimate was provided already in our first paper [1]). In order to simplify the discussion, we merely discard the whole string of non-minimal couplings of the (charged) nucleon to the external field, since they only give corrections to the Landau levels (in the zero frequency limit).
We look for a stationary solution of the Dirac equation (see, e.g., [32]): where F µν = ∂ µ A ν − ∂ ν A µ denotes the electromagnetic field strength tensor and g = +1, 0 for the proton and the neutron, respectively. Writing the wave function ψ as we obtain: with p = −i∇. Further, it is convenient to consider the non-relativistic limit, in which the mass m is the largest term on the l.h.s of (76). Assuming that eB ≪ m 2 , the function G can be easily expressed from Eq. (76) and hence one gets an equation for F (x): For the field configuration given in Eq. (12), the solution F (x) can be searched by using the ansatz where p 1 , p 3 are the conserved components of the three-momentum, and f is also an eigenvector of the σ 3 matrix: Using now the periodic field configuration from Eqs. (12) and (13), Eq. (78) takes the form where we have set p 1 = p 3 = 0 to focus on the ground state.
Let us first consider the case of the proton with g = 1. Introducing a new variable z = ωx 2 2, this equation can be brought to the standard form of the Whittaker-Hill equation: where a = 4 According to Floquet's theorem (see, e.g., Ref. [33]), Eq. (82) has solutions with the pseudoperiodic property where ν denotes the characteristic exponent. As is well known, solutions are bounded only for certain values of a (the parameters p, q are fixed), which form the band structure. Only in this case, ν has a vanishing imaginary part.
To see how the Landau levels emerge in the limit ω → 0, it is useful to go back to Eq. (81).
It is clear that, in the limit ω → 0, the cosine in the last term on the left-hand-side can be replaced by unity. Then, Eq. (82) takes the form of the Mathieu equation: where z ′ = 2z and A = (a − 2q) 4, Q = −p 4. In Fig. 1 we show the stability chart for this equation. In particular, for the colored regions in the A−Q plane, the solutions are bounded.
The band structure is clearly seen for Q ≠ 0. As Q → ∞, each band smoothly transforms into a Landau level. We note that the stability chart for Eq. (82) will be slightly different, but the picture is very similar, in particular, in the large Q region. This qualitative result can be explicitly verified by using the asymptotic expansion of the eigenvalues A n (n = 0, 1, . . . ) for large Q (see, e.g., [33]): or, The same conclusion can be drawn in a finite volume. In this case, the solutions in a band (the Bloch wave functions), do not obey, in general, periodic boundary conditions.
The requirement of periodicity, f (z + π) = f (z), picks out one level from the band. Such levels form the so-called characteristic curves A n (Q), which eventually approach the Landau levels as Q → ∞. To this end, we consider another limiting case, when B → 0 while the frequency ω is fixed.
Since q = O(B) and p = O(B 2 ), from Eq. (82) we again arrive at the Mathieu equation: The structure of the finite-volume spectrum is shown in Fig. 2. The expansion of the ground state level a 0 in small q reads [33] a 0 = − A similar formula can be written by using exact power series for a 0 in small p and q, see Ref. [34]. In particular, the first terms in these expansions coincide. One obtains for the proton and neutron, respectively: with κ p = 1.79 and κ n = −1.91. In case of the neutron, the numerical value of the radius of convergence reads q crit ≈ 1.47 (see also Ref. [35]). For an estimate we use the same value for the proton. This is justified since in the limit B → 0, Eq. (82) can be well approximated by the Mathieu equation. Accordingly, one obtains eB < 0.26ω 2 (proton) and eB < 0.38ω 2 (neutron).
The upper bound on the magnetic field strength in Eq. (90) can be improved by noting that our perturbative result, Eq. (61), might not be applicable at q = q crit . To estimate higher order corrections, we can again resort to Eq. (89). In Fig. 3 we plot the function a 0 (q) as well as the first term in Eq. (89). As is seen, the higher order terms become large at q = q crit and hence can not be neglected anymore. Accordingly, it is plausible to choose a certain value q = q max , for which they, e.g., amount to a 10% of the leading piece. This gives q max ≈ 1.05 (see also Fig. 3). Using Eq. (90) we get an improved bound on the magnitude B: eB < 0.19ω 2 and eB < 0.27ω 2 for the proton and the neutron, respectively.
It interesting to note that Eq. (89) allows us to verify the main result, Eq. (61), in the approximation where the proton is treated as a point-like particle but with a non-zero anomalous magnetic moment. This is equivalent to setting F 1 = 1, F 2 = κ and β M = 0 in Eqs. (10,11). The subtraction function takes the value It is seen from Eq. (89) that there is no spin-dependent contribution to a 0 . Using Eq. (83), we get directly the spin-averaged energy shift: where we have used the relation E 2 − m 2 ≈ 2mδE. This is precisely the expression which is obtained from the main formula in Eq. (61), when we substitute the subtraction function given in Eq. (91).
The above discussion is equally applicable for the neutron, for which g = 0. In particular, the differential equation for F (x), Eq. (78), simplifies: It can be brought into the form of the Mathieu equation: where As expected, no Landau levels emerge in the zero frequency limit ω → 0. Further, the main formula in Eq. (61) can be verified in a similar manner by setting F 1 = 0, F 2 = κ and β M = 0 in Eqs. (10,11). The spin-averaged energy shift reads

IV. PROPAGATOR IN AN EXTERNAL FIELD
The derivation based on the non-relativistic framework, which was given in the previous section, is not well suited for the study of the (exponentially suppressed) finite-volume effects. For example, the final result, displayed in Eq. (61), contains the infinite-volume Compton amplitude on the right-hand side, i.e., the finite-volume effects are neglected there.
In the non-relativistic framework, these effects may emerge from different sources. with the external field in a finite volume uniquely determines all these couplings. For a detailed discussion of these issues we refer the reader to Refs. [36][37][38][39][40][41].
From the above discussion it is clear that, in order to evaluate the finite-volume effects, it is better to work directly with ChPT in a finite volume, abandoning the non-relativistic framework, which has proven very convenient for discussing the zero-frequency limit. The exponentially suppressed finite-volume effects, which are not taken into account in ChPT, go as exp(−Λ H L) instead of exp(−M π L) (here, Λ H denotes a typical hadronic scale of order of one GeV), and thus can be neglected.
One important remark is in order. A procedure, which is used in Refs. [36][37][38][39][40][41] for the extraction of the polarizabilities in a finite volume, boils down to the derivation of the finitevolume one-particle effective action and to the identification of the different terms in this action. In this way, one again encounters the problem with operators containing Wilson lines that makes e.g., the authors of Ref. [38] to conclude that "At finite volume, there is no longer a discernible relation between polarizabilities and the Compton tensor." In our framework we shall choose a different path, directly relating the amplitude for forward Compton scattering in a finite volume to the second-order energy shift of the nucleon in the external magnetic field. We are not asking ourselves, what the subtraction function is in a finite volume -this question anyway does not have an unique answer and, making an inconvenient choice, one can easily obscure the relation between the infinite-and finitevolume quantities. Rather, we can uniquely identify the quantity that is extracted from the nucleon energy shift in a finite volume and which reduces to the subtraction function in the limit L → ∞ (as we shall see below, this is a certain component of the spin-averaged Compton tensor in a particular kinematics). This fully suffices to define a finite-volume counterpart of the subtraction function S 1 (q 2 ) and to calculate finite-volume corrections in an unambiguous way.
In this section, using the framework of the effective field theory in a finite volume, we shall derive the expression for the nucleon energy shift in an external field (a finite-volume analog of Eq. (61)). Note also that we shall never specify the Lagrangian of this theory -it is only used to catalyze the proof and produce the diagrammatic expansion of all amplitudes in terms of hadronic propagators.
We start from the nucleon two-point function in the external field in Minkowski space and define:D where Ψ(x) denotes the four-component spinor field, describing the nucleon. Note that the Dirac indices are not shown explicitly. Since the external field does not depend on time, we haveD(x, y) =D(x 0 − y 0 ; x, y). Further, it is convenient to define the Fourier transform in the fourth componentD as well as in vector components, One can invert this expression, giving The free propagator takes the form: where m is the physical nucleon mass. The diagrammatic representation of the propagator in the external electromagnetic field is schematically shown in Fig. 4. We have where the self-energy part Σ(p, k; E), which is a matrix in the space of Dirac indices, can be expanded in powers of the magnitude of the external field: Here, Σ 0 (p, k; E) = L 3 δ pk Σ 0 (p; E) is the sum of all one-particle irreducible diagrams in the absence of the external field, and the functions Σ 1 , Σ 2 will be determined below. Note also that we use the MOM renormalization scheme, i.e., Σ 0 (p; m 2 + p 2 ) = 0. The sum in Eq. (102) can be written in a compact form: where the amplitude T (p, k; E) satisfies the relation which is similar to the Lippmann-Schwinger equation.
In order to find the energy shift of the nucleon ground state, we have to determine the pole position in the propagator D(0, 0; E). For this purpose, we single out the term with l = 0 in the sum (corresponding to the unperturbed ground state) and rewrite the amplitude T (0, 0; E) as follows: where the quantity T ′ (p, k; E) satisfies the equation Note that now the sum runs over all l ≠ 0. We get Here, I denotes the unit 4 × 4 matrix. Inserting this expression into Eq. (104) for the propagator D(0, 0; E), one obtains Obviously, D(0, 0; E) and T (0, 0; E) have the same pole structure.
In order to simplify the matrix equation (106), we can use the octahedral symmetry of the cubic lattice. For zero momenta, the symmetry requires that: where g denotes an arbitrary element of the octahedral group and R αβ (g) is the matrix of the linear representation of the octahedral group which is obtained by restricting the (1 2, 0) + (0, 1 2) representation of the Lorentz group to its octahedral subgroup (here the Greek letters denote Dirac indices). The requirement of invariance restricts T αβ (0, 0; E) to the form: where T (1) and T (2) are scalar functions. Accordingly, we see that Similar relations can be established for the amplitude T ′ (0, 0; E). For example,T ′ is defined through the equationū Further, it is convenient to write the free propagator D (0) αβ (p; E) in the form Multiplying the Eq. (106) byū(0, s) from the left and by u(0, s) from the right, we get: This equation is not a matrix equation anymore. The pole position is given by As expected the free pole at m − E = 0 has disappeared. This result is similar to the "master equation" in case of hadronic atoms [42].
Next, let us proceed with the calculation of the amplitude T ′ (0, 0; E). In perturbation theory, up-to-and-including O(B 2 ) this quantity reads The quantity Σ 1 (p, k; E) can be expressed through the three-point vertex function. Indeed, consider the linear coupling to the external field which is described by the Lagrangian Here, we have used Eq. (13) which states that the vector potential has only one nonzero Further, using translational invariance, one can write This is a definition of the vertex function Γ(p; k). Substituting this expression into Eq. (121) and integrating over u 0 and k 0 , we get where Γ(p, k; E) is obtained from Γ(p, k) by substituting p 0 = k 0 = E (we remind the reader that the field A 1 is static). Accordingly, the Fourier transform ofD 1 (x, y) takes the simple form: Comparing this result with the expansion of the propagator in Eq. (102) at O(B), it is seen Next, we evaluate the quantity Σ 2 (p, k; E) which consists of all one-particle irreducible diagrams with amputated nucleon legs, with two external fields attached. Let Υ(p, k, l) denote the sum of all such diagrams in momentum space. Here, p and k denote the momenta of the outgoing and ingoing nucleon, respectively, and the momenta of two external "photons" are equal to l + (k − p) 2 and l − (k − p) 2, respectively. Further, denoting Υ(p, k, l; E) = Υ(p, k, l) it is easy to check that Σ 2 is given by We now have all ingredients for the evaluation of the energy shift. Here, following the discussion in the previous section, we consider the scenario a) for the external field, when the frequency ω is quantized and the three-momentum conservation holds. In this case, the term linear in B gives: The second-order term at threshold takes the form On the other hand, it is straightforward to verify that the expression on the right-hand side of this equation is proportional to the "11" component of the forward Compton scattering tensor. Note also that, since ω ≠ 0, the sum over the intermediate nucleon states does not contain the term with l = 0 and hence, there is no difference between T 2 and T ′ 2 at threshold. Thus where p µ = (m, 0), q µ = (0, ω) and the nucleon wave function renormalization constant (in the absence of the external field) is given by It is now straightforward to determine the spin-averaged energy shift at order B 2 from Eq. (118). Taking into account the fact that the linear term in B vanishes, and expanding the quantityT ′ in Eq. (118) in Taylor series in E − m, it immediately follows that the spin-averaged energy shift is given by from which we finally obtain This equation is the finite-volume version of Eq. (61) and contains the "11" component of the Compton tensor, evaluated in a finite volume. Up to the corrections, proportional to exp(−M π L), this quantity is given by the subtraction function S 1 (q 2 ) via T 11 (p, q) = −ω 2 S 1 (q 2 ), see Eqs. (6) and (7). Hence, the equation (133) provides a framework for the systematic calculation of such corrections.

CALCULATIONS
In Ref. [1] we presented a brief discussion of the lattice parameters, which could be used in the numerical extraction of the subtraction function S 1 (q 2 ). With the use of the new, more stringent constraints on the value of the magnetic field we are now able to refine this analysis. For a moment, we neglect the finite-volume corrections altogether. As in Ref. [1], the elastic and inelastic parts of the amplitude are parameterized as: where κ and β M denote the anomalous magnetic moment and the magnetic polarizability of the nucleon (we use the same numerical values in the estimates as given in Ref. [1]). The electric and magnetic form factors of the proton and the neutron are given by: with the same dipole form factor as above.
One of the estimates, which goes through exactly in the same way as in Ref. [1], is related to our ability to separate the physically interesting inelastic part from total amplitude. As seen, the elastic part is singular at threshold and falls off very fast at higher ω 2 . Thus, the separation will be difficult for very small values of ω 2 . One may require, for instance, that at the minimum value of ω 2 , the inelastic contribution amounts up to a 10% of the elastic contribution. In this manner, we get ω 2 min = 0.086 GeV 2 for the proton and ω 2 min = 0.045 GeV 2 for the neutron (a slight difference to the numbers given in Ref. [1] is caused by the fact that here we take into account the exact momentum dependence of all amplitudes). If one requires instead that elastic and inelastic parts are equal, one gets ω 2 min = 0.40 GeV 2 for the proton and ω 2 min = 0.26 GeV 2 for the neutron. In any case, the lower cutoff on the available frequencies is rather comfortable and does not put significant restrictions on the parameters of the lattices which can be used in the calculations.
The conditions, which involve the magnitude of the magnetic field, are more restrictive.
On one side, the magnetic field should be strong enough, in order to measure the effect at all. On the other hand, it must be weak enough, so that the perturbation theory still applies. In section III F we have made a more stringent estimate It is now clear that the window for the available values of eB exists if and only if δE inel can be taken sufficiently small or, in other words, if the uncertainty in the determination of δE inel does not exceed certain value. Generously allowing this uncertainty to be δE inel = 0.05m, as done in Ref. [1], is no more an option -the inequalities in Eq. (137) can not be satisfied. A better accuracy in the determination of the energy shift δE inel = 0.01m would lead to the lower cutoff on the available frequencies ω 2 min = 0.90 GeV 2 for the proton and ω 2 min = 0.41 GeV 2 for the neutron. If one wants to increase the range of available frequencies and reach lower values of ω 2 , one has to improve on the accuracy further. The range of the magnitudes for the magnetic field can be determined from the above equations, if the accuracy is given.

VI. CONCLUSIONS AND OUTLOOK
i) We have presented three alternative derivations (the third one is contained in the appendix) of the formula for the energy shift of the nucleon, placed in a periodic external field. Namely, the non-relativistic effective Lagrangian was used for this purpose, as well as the relativistic framework. All these alternative settings are advantageous for discussing different issues arising in the treatment of the problem. The aim of the whole exercise is to extract the forward Compton scattering amplitude in a certain kinematics, the so-called subtraction function S 1 (q 2 ), from lattice simulations. In its turn, measuring this function would enable one to gain important information about the properties of QCD at low energy.
ii) The result or Ref. [1] has been refined and extended in various aspects. For example, in this paper we discuss in detail the zero-frequency limit (constant magnetic field) of the expression for the energy shift. This limiting case is studied in the literature in detail.
Here, it is shown that, to have a smooth transition to this limit, the external field on the lattice should be implemented in a specific way, corresponding to the scenario b).
There exists no zero-frequency limit for scenario a).
iii) Another important issue is the limit of validity of the perturbative treatment of the external magnetic field (the convergence radius of the perturbative expansion in B).
In Ref. [1], using heuristic arguments, we gave a rough estimate of the maximal field strength, for which the perturbative approach is still applicable. In the present paper we improved the argument and give a new, much more stringent estimate, which is based on the properties of the solutions of Mathieu's equation. It should be pointed out that the non-relativistic EFT approach provides the most convenient framework for the discussion of the above two problems. For the sake of completeness, we present yet another derivation of the main formula for the energy shift given in Eq. (133), which is based on the study of the two-point function of the nucleon field in the external field at large time separation and, hence, has a closer resemblance to the methods used on the lattice. It is assumed that the external field is implemented according to scenario a), i.e., the frequency is quantized. In this derivation, we directly expand the nucleon two-point function in the external field Here, the integration over d 3 xd 3 y projects onto the states with zero initial and final threemomenta.
Note that, strictly speaking, for a rigorous derivation one should perform the Wick rotation into the Euclidean space and pick up the leading terms in the two-point function at large (Euclidean) times. For simplicity, however, we stay in the Minkowski space and identify the leading exponentials there -in our case, the identification is easy and no ambiguities arise.
The completeness condition, which we shall be using, takes the form 1 L 3 ks k, s⟩⟨k, s 2ω(k) where ellipses stand for the excited states contributions; they will be neglected altogether by taking the limit x 0 − y 0 → ∞ (the time extent of the lattice is assumed to be infinite).
Let us start with the matrix element C (0) that describes the propagation of the nucleon in the absence of the external field. Using the translation invariance, it takes the form Note that the second term in the T -product, containing θ(y 0 − x 0 ), picks up the antiparticle pole for large x 0 − y 0 instead of the particle pole, and thus can be neglected.

(A.13)
After the summation over the three-momentum the above expression simplifies and we get: The integral over the periodic electromagnetic potential vanishes, and so