Charged pion electric polarizability from lattice QCD

We present a calculation of the charged pion electric polarizability using the background field method. To extract the mass-shift induced by the electric field for the accelerated charged particle we fit the lattice QCD correlators using correlators derived from an effective model. The methodology outlined in this study (boundary conditions, fitting procedure, etc.) is designed to ensure that the results are invariant under gauge transformations of the background field. We apply the method to four $N_f=2$ dynamical ensembles to extract $\alpha_{\pi^\pm}$ at pion mass of $315$ MeV.


I. INTRODUCTION
Electromagnetic polarizabilities are important properties that shed light on the internal structure of hadrons. The charged quarks inside a hadron respond to applied electromagnetic fields, revealing the deformation (or 'stiffness') of the composite system as a fraction of its volume for uniform fields, and charge and current distributions for varying fields. There is an active community in nuclear physics partaking in this endeavor. Experimentally, polarizabilities are primarily studied by low-energy Compton scattering. On the theoretical side, a variety of methods have been employed to describe the physics involved, from dispersion relations [1], to chiral perturbation theory (ChPT) [2][3][4] or chiral effective field theory (EFT) [5],to lattice QCD. Refs. [2,5] also contain reviews of the experimental status.
In this work, we focus on the electric polarizability of charged pions in the background field method, which presents its own challenge because the entire system accelerates in the presence of electric field. This motion is unrelated to polarizability and must be isolated from the deformation due to quark and gluon dynamics inside the hadron. Standard plateau technique of extracting energy from the large-time behavior of the two-point correlator fails. Our goal is to find a robust methodology to take the motion into account so that the polarizability can be reliably extracted. In this study we employ Dirichlet boundary conditions and address a number of related theoretical issues.
In Section II, we discuss the effective model we use to capture the behavior of the charged pion. In Section III, we describe our fitting procedure incorporating the effective model and present results on the polarizability from lattice QCD simulations. Conclusions are given in Section IV.

A. Background field method
To extract the polarizability, we compute the pion mass shift induced by an external electric field. For technical convenience we use imaginary electric fields and the mass shift is then The imaginary background electric field is introduced by multiplying the SU (3) gauge links by the U (1) links Note that here the charge q is the charge of the relevant quark field: for down quarks q d = −e/3 and for up quarks q u = 2e/3, with e the absolute value for the electron charge. There are an infinite number of U (1) gauges that produce the same electric field; one such gauge is the t-dependent vector potential on the x-links, where t 0 is the origin of the potential. Since the calculation will be performed in a gauge-invariant manner, the value of t 0 can be shifted via a gauge transformation; we set t 0 = 0. The potential produces a constant electric field in the x direction except at the edges of the lattice.
Note that the interaction with the background field is physical only if the electric field is real. In this case the electric field leads to a real-valued exponential factor in Eq. (2). Using an imaginary electric field is convenient since in this case the factor in Eq. (2) is a complex factor of modulus one and the implementation for the discretizied Dirac operator is basically unchanged. Using imaginary fields is justified because the theory, with Dirichlet boundary conditions, is analytic for small values of E and we can use analytically continuation to make the fields imaginary. The two formulations are equivalent as long as analytic continuation is valid and the sign is corrected in the interpretation of Eq. (1) (see discussion of the issue in Refs. [13,31]). We note in passing that there is no such ambiguity for introducing a constant magnetic field on the lattice where only spatial coordinates are involved. A real magnetic field leads to a U (1) phase factor and the magnetic polarizability term in the mass shift has the negative sign: ∆m = −βB 2 /2.

B. Relativistic effective correlator
In the presence of the background field the two-point correlator for charged particles is not expected to be a simple sum of exponentials, as is the case for the neutral particles. A work around is to fit the hadron correlators against the functional form expected for a charged particle propagator in the presence of an electric field [32]. As it turns out the proposed correlation function computed in the infinite-volume is not a good fit for our correlators since the finite volume effects are important. We proposed a method that incorporates these effects [33]. In this section we outline the method used to compute the fit correlator. We show that our correlator converges in the infinite volume limit to the expected correlator and that the finite volume effects are important.
To extract polarizabilities, we fit the lattice QCD correlator for pions to the correlator for a relativistic scalar particle in two-dimensions. As we discuss later in the paper we use Dirichlet boundary conditions in the direction of the electric field. For transverse directions we use periodic boundary conditions and we expect that the low energy pion states have zero momentum in the transverse directions. As such calculating the particle propagator in 3 + 1 dimensions reduces to a two-dimensional problem. We start with the continuum action in Euclidean space for a massive scalar particle where φ is a complex scalar field. A background electromagnetic field can be introduced via a vector potential A µ through minimal coupling to the charge, Note that here the charge corresponds to the pion charge q = e. Our effective model is a two-dimensional, discretized version of the action given by [34] with where n denotes the sites on the lattice, a is the lattice spacing, andμ accounts for the nearest neighbors of every site on the lattice, i.e. ±x and ±t. The two-point correlation function is the inverse of the K-matrix: This effective correlator can be compared to its infinitevolume, continuum version which can be derived using Schwinger's proper time method [11,35] The comparison is shown in Fig. 1 for the effective mass function where the continuum and infinite volume limits are taken. The effective model correlator is computed using Dirichlet boundary conditions. We see that our effective model converges to the correct infinite-volume, continuum version, with or without background electric field. The effective mass in the presence of the field is not constant even at large times. In fact, it grows as a function of time, probably due to the acceleration of the particle in the field. As we advertised earlier, a different fitting procedure must be developed for the new functional form. We also note large differences between the finite-volume and infinite-volume forms for values of mL ≈ 5 encountered in typical lattice QCD simulations (in this study we focus on the pion and m is the pion mass). Since we are relying on tiny mass shifts to extract polarizabilities, such differences can have a large impact on the results. We will develop in this work a fitting method incorporating the full content of the effective correlator and test it out on a number of lattice QCD ensembles.

C. Boundary conditions
In a finite volume the background field method requires careful handling of the external field at the boundaries.  (9). For plots on the left we study the continuum limit: we make the lattice spacing finer and finer at fixed box size until convergence. For plots on the right, we focus on the finite volume effects. For the lower plots we use a value of the electric field of a 2 qE = 0.005. For all plots the effective correlator is computed using Dirichlet boundary conditions in both directions. The lattice size in the time direction is Lt = 10m −1 and the source for the propagator is at t0 = m −1 .
Even though the electric field is constant, the electrostatic potential A µ has discontinuities at the edges of the box. For real electric fields, these discontinuities cannot be avoided. One workaround is to use imaginary fields with quantized values such that the ratio qE/(2πL x L t ) is an integer. This ensures uniform electric flux through all the xt plaquettes and while the electrostatic potential is still discontinuous. One problem with this approach is that there is no obvious analytical continuation that connects this setup to a situation involving physical electric fields, so it is unclear whether the results computed in this setup are directly relevant for polarizability. Another problem is that in this setup different equivalent gauge choices for the external fields lead to results that are not connected by a gauge transformation. For example, in a gauge similar to the one used in Ref. [24], of the form corresponding to Eq. (3), the Polyakov loops in the spatial directions depend on the choice of t 0 and, since the Polyakov loops are invariant under gauge transformations, different t 0 choices could lead to different results. This might not be a problem if the two gauges produce the same results, but it has to be verified, a posteriori. We choose to work with Dirichlet boundary conditions in both x and t directions. In this setup gauge invariance in the external field is guaranteed as we will discuss later. We can use either real or imaginary electric fields, and non-quantized values for the field strength, as small as needed for the quadratic terms due to polarizability to dominate.
On the other hand, using Dirichlet boundaries introduces a number of issues that must be dealt with. First issue is that the lowest energy state corresponds to a moving hadron. This is the main reason for the discrepancy present in the upper-right panel of Fig. 7 between the effective mass extracted from the finite-volume correlator and the infinite-volume one. For neutral hadrons, when extracting the energy from exponential decay of the correlators, we need an estimate of the hadron momentum to correct for the effect on the mass shift [7]. When fitting the lattice QCD correlator to the one from the effective model, this ajustment is done automatically and the value extracted is directly the rest mass of the particle.
Another important effect is due to the pion interactions with the Dirichlet walls. Since hadrons are extended objects, they interact with the walls differently than the point-like particles described by the effective model. Away from the walls the hadrons move freely and the centerof-mass wavefunction corresponds to a free particle wave, but this does not hold true near the walls. The net effect is that the free wavefunction for the hadrons away from the walls has nodes that are not exactly aligned with the walls. To account for this difference we use a different distance between the walls in the effective model. The difference is related to the reflection scattering length introduced to capture the effect of this interaction for extended objects [36,37]. This effect will be studied using the pion correlator in the absence of an external field.

III. TECHNICAL DETAILS AND RESULTS
To fit the effective model we need to construct the hadron wavefunction profile between the Dirichlet walls. This requires us to calculate lattice QCD correlators as a function of both x and t, the directions where we use Dirichlet boundary conditions. In the y and z directions we use periodic boundary conditions and we project the correlator to zero momentum in these directions, whereÔ =dγ 5 u is the interpolator for π + in the twopoint function, is (x s , y s , z s , t s ) is the source for the quark correlators. Note that due to isospin symmetry between the u and d quarks, π + and π − (Ô =ūγ 5 d) have identical polarizability.
For the effective model we use the propagator derived in Eq. (8) Above we emphasized that the correlators are functions of the external field A µ and the distance between the walls: L for the effective correlator and L for QCD calculation.
We have generated the correlators with the position of quark sources fixed at x s /a = N x /2 + 1 and t s /a = 6. We use a coordinate system where the Dirichlet walls are positioned at x = 0 and x = L in the lattice QCD box. The effective model is aligned such that the source is at the same position and the walls are at x = (L −L)/2 and x = (L +L)/2. In effect the setup keeps the source position the same and changes the distance to the xdirection walls by the same amount. For both QCD and the effective model the time boundaries are t = 0, and t/a = (N t + 1). For QCD correlator, we take advantage of translation invariance of the system in y and z directions, and use sources at 64 locations for each configuration to obtain more statistics. The sources are spaced regularly at 8 positions on a y-z grid and we use 8 different time translations. Note that since the gauge fields are generated with (anti-)periodic boundary conditions in time, time-translation is a symmetry of the ensemble. We apply the Dirichlet boundary conditions in time after the translation. A similar symmetry exists also in the x-direction but we did not take advantage of it in this study.
To find the appropriate distance between the walls for the model, we fit the lattice correlator to the effective model using different sizes for the model. We fit the correlators away from the walls to insure that the difference in the interactions with the walls do not affect the fit. The fit ranges, in both x-and t-directions are reported in Table I.

A. Fitting procedure
To extract the mass from the lattice correlator we determine the mass in the effective model that minimizes the distance between the QCD correlator G and the effective one G 0 . The distance between the two correlators is defined via a χ 2 function. The design of this distance function was guided by two main concerns: the fit procedure should produces the same results for all possible equivalent external field gauges and we want to accommodate the fact that the correlator is a complex valued function.
Assume that we perform on each configuration the measurements y i that are complex valued. The effective model value that is designed to match y i is assumed to be f i . We define the χ 2 -function (13) and the covariance matrix The covariance matrix is hermitian and the χ 2 -function is real and positive-definite. The minimization process will vary the parameters of the effective model to minimize We are now interested in the action of the gauge transformations on the distance functions through the observables y i . Assuming that the transformation acts linearly on the observables such that y = T y (that is y i = T ij y j ), with T the same matrix for all configurations. Then y = T y and if we have f = T f , then δ = T δ. Similarly for the covariance matrix we have C = T CT † and the distance function is then invariant, that is (χ 2 ) = (δ ) † (C ) −1 δ = χ 2 . This condition is satisfied if the observables y i are any subset of the pointto-point propagator function G(x, t; A µ , L) since under a gauge transformation, (A, φ) → (A , φ ) with we have G(x, t; A µ , L) = e iq[Λ(x,t)−Λ(xs,ts)] G(x, t; A µ , L) . (16) Above we assumed that we use imaginary electric field. The same relation holds for the effective theory propagator G 0 so the distance between the point-to-point propagators is invariant under gauge transformations. Given the discussion above it is then tempting to fit the lattice propagator at a subset of the (x, t) positions, in the range included in Table I. The problem is the presence of high frequency modes in the propagators close to the ultraviolet cutoff that are distorted due to lattice discretization. These modes can have an important contribution to the χ 2 -function. The smallest eigenvalues of the covariance matrix C will have the biggest contribution to the χ 2 . This means that the lowest eigenvalues will drive the fitter and force the minimizer to a point in the parameter space that matches the distorted high frequency modes in the lattice QCD data instead of the physical long-range two-point function we are interested in. This is indeed the case for our calculation. In Fig. 2 we show the modes corresponding to the lowest and the highest eigenvalue in the covariance matrix, when we take the observables to be all space-time points in the range indicated in Table I for ensemble EN3. We see that the most important modes for the fitter, the ones that correspond to the lowest eigenvalue of the correlation matrix, are the ones that fluctuate rapidly along the spatial directions, the distorted high-frequency modes. To remedy this problem, we need to filter out the these high frequency modes from our observables.
For the zero-field propagators, we can just sum over the spatial position to filter out the high-frequency modes. This is akin to zero momentum projection used when fitting propagators with periodic boundary conditions. For our case this is only used to reduce the influence of the high-frequency modes, not to project to a particular momentum state. Our summation only extends over the interior points, away from the Dirichlet walls. However, this does not work for non-zero field propagators since gauge transformations act non-linearly on the observable defined via this summation. A different form is required to preserve gauge-invariance. In this work we use the following observables: where  and similarly for f (t), and the distance function becomes gauge invariant. The energy shift induced by the background field is very small, much smaller than the stochastic error on the extracted hadron masses. To extract the shift reliably we need to take into account the correlations between zero-field and non-zero field correlators. We perform simultaneous fits of the zero-field and non-zero field correlators. The residue vector δ = {δ 0 , δ E } includes the zero-field residue δ 0 and the residue in the presence of the field δ E defined using the fit form The enlarged covariance matrix can be written as The enlarged χ 2 remains gauge invariant. Above we made explicit the dependence on the mass parameter. The fit parameters are A, am, ∆A, and ∆am. The mass shift ∆am is used to compute the polarizability.

B. Effective range
As discussed earlier, the effective sizeL, the distance between the Dirichlet walls in the effective model, is slightly different from the distance between the walls in the QCD ensemble due to the interaction between the hadrons and the walls. To determineL we use zerofield lattice QCD correlators. The ensembles used in this study are summarized in Table I. For details about the generation of the ensembles see [38].
For each ensemble, we have done multiple fits to map out the dependence of the extracted mass onL. We choose the x and t ranges as wide as possible while maintaining a reasonable χ 2 , and look for an effective sizeL that reproduces the mass obtained from periodic boundary conditions on the same ensembles [39]. For the three larger ensembles the finite-volume corrections for the pion mass are small and we used am π = 0.1936 (2). For the smallest ensemble, EN1, we used am π = 0.1986(22) [7]. The error-bars onL were determined by jackknife. The dependence of the extracted mass onL for our ensembles is shown in Fig. 7. We will discuss later the sensitivity of the mass-shift extraction on the effective size. TheL values determined from our ensembles and the fit ranges used are added to Table I and will be used in all of our subsequent fits. Note that theL values were rounded up to the nearest half lattice-spacing to make the effective model calculations easier.
To see how well the effective model describes the lattice QCD correlators, we show in Fig. 3 some snapshots of the spatial profile (wavefunction) in the absence of electric field. Indeed, the Dirichlet walls force a standing wave in the x direction. For the effective model there are no excited states in contrast to QCD where other excitations are present. However, in the pion channel the excited states lie much higher and the ground state pion becomes dominant as we move away from the source. We can see from the figure that this is indeed the case: the two correlators match very well as time progresses. This suggests that our effective models correctly captures the relevant dynamics. We note here that the correlators match well only after we adjust the effective distanceL appropriatedly. The same match is observed in the time decay of the correlation at fixed x locations, as shown in Fig. 4.
In the presence of the background electric field the correlators become complex. We show in Fig. 5 and Fig. 6 the real and imaginary parts separately of the charged particle correlator. For the effective correlator we use the same parameters, am and A, as in the zerofield case above. The agreement is very good since the corrections, ∆am and ∆A, associated with the external field (a 2 q d E = 10 −2 in this case) are minute. Interestingly, the wavefunction develops nodes between the walls in the presence of background electric field. It is tempting to think that this is due to the acceleration in the presence of the field, but in fact this is due to the gauge choice for the external field.

C. Polarizability extraction
To determine the pion polarizbility, we computed the lattice QCD correlators G(x, t) on the four ensembles in Table I for charged pions, and carried out correlated fits FIG. 7. Dependence of am on the distance between the Dirichlet walls for zero-field case. This is used to determine the effective lattice sizeL in Table I. as described above. We have used the same fit region that we had used when determiningL of the effective model in Table I. Furthermore, for ensembles EN2, EN3, and EN4 which have the same temporal extension, we have fixed the temporal fit range to be the same for all of them.
To make sure that we are extracting the quadratic response to the electric field, we compute the mass shift for a sequence of increasing values for the electric field strength. In Fig. 8 we plot the mass shift as a function of field strength for the smallest ensemble (EN1) and the largest (EN4). We see that the value of the electric field used in this study is indeed in the regime where the mass shift scales quadratically with the electric field strength. For our study we used the smallest value of the electric field where the mass-shift is not statistically compatible with zero on the smallest ensemble (EN1): a 2 q d E = 0.01.
The fit results are summarized in Table II. The mass shifts are indeed small when compared to the particle mass. The conversion from lattice units to physical units from Eq. 1 is given by where a = 0.1245 fm, e 2 = 1/137, and a 2 q d E = 0.01. The dominant sources of systematic errors are the choice of the fit range in the spatial and temporal direction, and the value ofL using in the effective model. To determine their contribution we analyze the sensitivity of the polarizability to expected changes in these parameters. For the fitting range, we reduced both the x and t-ranges by one unit of lattice spacing. ForL we vary it within the error bands reported in Table I. The results are included in Table III. We note here that for the three largest ensembles used in our study the finite-volume corrections estimate from ChPT is less than 10% [40]. In Fig. 9 we plot the results for polarizability as a function of the lattice extent in the field direction. Our results are compared with the ChPT value of [41]  For the largest volume (EN4) the polarizability turns, surprisingly, negative. We do not have an explanation for this behavior but we note that there are other lattice calculations of pion polarizability for comparable quark masses that are at variance with ChPT expectations. One is on magnetic polarizability that is expected to satisfy β π + = −α π + , a calculation for m π = 296 MeV finds β π + = 0.35(11) × 10 −4 fm 3 . Another calculation of the "connected" neutral pion electric polarizability-where the disconnected diagrams are not included-finds that the polarizability turns negative for m π < 350 MeV [7,8] although the ChPT expectations are that this value should be positive [32].

IV. CONCLUSION
Computing electric polarizability for charged hadrons has been a challenge for lattice QCD due to the acceleration experienced by charged particles in the background  electric field. Using the lightest charged hadron, the pion, as an example, we presented a method to extract the electric polarizability. The method relies on an effective field model that is designed to capture the behavior of the charged pion in a box with the same geometry as that used in lattice simulations. This has proven crucial since the infinite-volume version of the model fails to account for the significant finite-volume effects. This is important to isolate the energy shift due only to the deformation of the hadron from the effects due to its motion in the electric field. To fit the lattice QCD correlators we need to compute its x and t dependence separately to match the effective model. We construct a χ 2 -function that utilizes information in both the real and imaginary parts of the correlator simultaneously and is invariant under gauge transformations of the background field.
Since we use Dirichlet boundary conditions, we are not limited to just quantized values under periodic boundary conditions. We can dial the electric field continuously and choose values that give better signals.
We extract the pion polarizability on a set of ensembles with different distance between the Dirichlet walls. For the pion mass m π = 315 MeV used in this study, the results on the smallest three ensembles are compatible with a small positive value for α, a value in agreement with ChPT predictions. However, on the largest ensemble the polarizability is negative, a puzzling result. We have checked the sensitivity of the polarizability on the fitting range and the effective distance between the walls employed in the effective model. The effect of varying these parameters is smaller than the stochastic error and cannot explain the shift to negative values for the polarizability. We note that there are other lattice calculations of electromagnetic polarizabilities at similar pion mass that in tension with ChPT predictions.
In this paper we focused more on the technical aspects of the extraction method. Looking to the future, in order to compare with experiments, we need to extrapolate to physical quark masses. We expect that as we approach the physical mass ChPT predictions become more reliable and should offer a strong check for our results. In particular we plan to investigate the puzzling result we got on the largest volume in this study. ChPT predicts a 1/m π rise as the pion mass is lowered (see Eq. 23) so the signal might be stronger at smaller pion masses. Another direction we plan to investigate is the effect of lattice discretization. As an additional check of the fitting procedure, we plan to apply the same method to neutral mesons on the same set of ensembles where simple exponential fits were previously employed [7,8]. Another systematic effect, namely charging the sea quarks, is not addressed here, though such effects are expected to be small [9,10]. In the long run, an extension to baryons would be desirable with an eye towards the proton electric polarizability which is more precisely measured than charged pions but is less well studied on the lattice.

ACKNOWLEDGMENTS
This work was supported in part by DOE Grant No. DE-FG02-95ER40907. AA gratefully acknowledges the hospitality of the Physics Department at the University of Maryland where part of this work was carried out. The computations were performed on the GWU Colonial One computer cluster and the GWU IMPACT collaboration clusters.