First extraction of the scalar proton dynamical polarizabilities from real Compton scattering data

We present the first attempt to extract the scalar dipole dynamical polarizabilities from proton real Compton scattering data below pion-production threshold. The theoretical framework combines dispersion relations technique, low-energy expansion and multipole decomposition of the scattering amplitudes. The results are obtained with statistical tools that have never been applied so far to Compton scattering data and are crucial to overcome problems inherent to the analysis of the available data set.


I. INTRODUCTION
Real Compton scattering (RCS) is one of the fundamental processes to access information on the internal structure of the nucleon. The RCS amplitude can be separated into a Born contribution, describing the scattering off a pointlike nucleon with anomalous magnetic moment, and a structure-dependent part referred as non-Born term. The non-Born contribution is parametrized by polarizabilities, which describe the response of the nucleon's internal degrees of freedom to an external electromagnetic field. In the low-energy expansion of the non-Born amplitudes, the leading-order effects are given by static polarizabilities, that are defined in the limit of zero frequency of the photon field and therefore measure the response to a static external electromagnetic field. The leading-order spin-independent polarizabilities are the scalar dipole electric and magnetic polarizabilities, α E1 and β M1 , respectively, while four spin-dependent polarizabilities appear at the next order and involve the nucleon-spin degrees of freedom. They have been the subject of intense research both experimentally and theoretically [1][2][3][4][5][6]. The currently accepted values for the proton scalar polarizabilities are α E1 = (11.2±0.4)·10 −4 fm 3 and β M1 = (2.5 ∓ 0.4)·10 −4 fm 3 [7], while the first extraction of the individual four spin polarizabilities has been obtained only recently from double-polarized Compton scattering [8].
As it is well known from many branches of physics, polarizabilities become energy dependent due to internal relaxation mechanisms, resonances and particle production thresholds in a physical system [9][10][11]. This energy dependence is subsumed in the definition of dynamical polarizabilities, that parametrize the response of the internal degrees of freedom of a composite object to an external, real-photon field of arbitrary energy. The enriched information encoded in the dynamical nucleon polarizabilities has been pointed out in different theoretical calculations, using dispersion relations or effective field theories [12][13][14][15][16]. In this work, we attempt for the first time to extract information on the scalar dipole dynamical polarizabilities (DDPs) from a fit to all available unpolarized RCS data below pion-production threshold. To this aim, we apply a statistical analysis based on the parametric bootstrap technique (see, for instance, [17] and references therein). Such a method has never been exploited to analyze RCS data and it is crucial to deal with problems inherent to both the low sensitivity of the RCS cross section to the energy dependence of the dynamical polarizabilities and to the poor accuracy of the available data sets.
A feasibility study for the extraction of two spin dynamical polarizabilities from unpolarized Compton scattering data has been presented previously in Refs. [18,19], following a different strategy from the present work. These fits did not turn out to be conclusive, mainly because of the scarce accuracy of the data set at disposal and the smaller sensitivity of the unpolarized cross section to the spin polarizabilities rather than the scalar polarizabilities. More recently, a partial-wave analysis of the unpolarized Compton scattering data has been discussed also in Ref. [20], pointing out the need to improve the accuracy of the experimental data set to pin down the values of the static scalar dipole polarizabilities with more precision.
The paper is organized as follows: In Sec. II, we introduce the theoretical framework to analyze the unpolarized Compton scattering cross section. In Sec. III we describe and motivate our fitting procedure based on the bootstrap technique, in comparison with the standard chi-squared minimization technique. Section IV contains our results for the fit of the scalar DDPs, and section V summarizes our conclusions.

II. THEORETICAL FRAMEWORK
The theoretical framework for the analysis relies on the multipole expansion of the scattering amplitude, the low-energy expansion (LEX) of the scalar DDPs and dispersion relations (DRs) for the calculation of the higherorder multipole amplitudes and the energy dependence of the scalar DDPs near the pion-production threshold. The definition of the dynamical polarizabilities rests on the multipole expansion of the RCS amplitude in the center-of-mass (cm) frame [1,[21][22][23]. The multipole amplitudes f l± T T ′ with T, T ′ = E, M , correspond to transitions T l → T ′ l ′ and the superscript indicates the angular momentum l of the initial photon and the total angular momentum j = l ± 1/2. In particular, for the scalar DDPs one has [13]: wheref indicates the non-Born contribution to the multipoles. In the calculation of the cross section, we compute the full Born contribution, given by pole diagrams involving a single nucleon exchanged in s-or u-channels and γN N vertices taken in the on-shell regime [1]. For the non-Born part, we use the multipole expansion. As also observed in Ref. [13], the multipole expansion of the non-Born contribution has a very fast convergence below pion-production threshold. In our analysis, we take into account the non-Born amplitudesf l± T T ′ up to l = 3: the scalar DDPs are fitted to the data, and the remaining contributions are calculated through subtracted DRs.
DRs have been proven to be a powerful tool to analyze RCS data [3,8,[24][25][26], as they allow to minimize the model dependence using as input available experimental information from other processes. The dispersion calculation is performed in terms of six independent invariant amplitudes A i , i = 1, . . . , 6, which can be recast in terms of the multipole amplitudes f l± T T ′ as explained in App. A of Ref. [13]. In the subtracted formalism, the six invariant amplitudes are obtained from subtracted dispersion integrals in both the s and t channels, and subtraction constants that are directly related to the six leading-order static polarizabilities. The subtracted integrals are saturated by πN , ππN and heavier-meson intermediate states in the s channel, and ππ intermediate states in the t channel [27]. In the present analysis, the input for the pion-photoproduction amplitudes has been updated to the most recent version of MAID [28], while we refer to [27,29] for the calculation of the contributions beyond πN in the s-channel and for the t-channel contribution. Four of the subtraction constants are fixed to the values of the static leading-order spin polarizabilities extracted in Ref. [8]. The two remaining constants are given in terms of the static dipole scalar polarizabilities as specified in the following.
We are interested in the energy dependence of the scalar DDPs below pion-production threshold, where they are real functions. By performing a LEX of the scalar DDPs, one recovers the limiting values of the static dipole polarizabilities at zero energy. Higher-order terms contain dispersive or retardation effects which can be parametrized in terms of higher-order polarizabilities [1,16,30].
The static polarizabilities are best defined via the effective non-relativistic Hamiltonian in the Breit frame, and hence the next-to-leading-order coefficients of the LEX of Eq. (1) are not given only in terms of higher-order static polarizabilities related to retardation effects of the E1 and M 1 radiation. The direct link is spoiled by recoil corrections in the cm frame. We have calculated the relevant recoil corrections up to O(ω 5 ), obtaining the following expressions for the LEX of the scalar DDPs: In Eqs. (2) and (3), the terms with even power of ω contain both retardation effects of the dipole radiations and recoil terms. The terms with odd powers of ω are recoil contributions, which, in addition to the contributions from the static polarizabilities of lower orders, can include terms with the static scalar and spin-dependent polarizabilities of higher multipolarity.
In particular, the first dispersive contributions enter at O(ω 2 ) and correspond to the static polarizabilities α E1,ν and β M1,ν . The recoil terms at O(ω 3 ) are given in terms of the static dipole scalar and spin polarizabilities, the fourth-order dipole scalar polarizabilities α E1,ν and β M1,ν and the quadrupole scalar polarizabilities α E2 and β M2 . In particular, α E1,ν , β M1,ν and the quadrupole polarizabilities enter as recoil terms with the same suppression factor in 1/M . The static spin dipole polarizabilities enter with a coefficient in 1/M 2 and the static scalar dipole polarizabilities enter with a factor in 1/M 3 .
At O(ω 4 ), the recoil terms contain different combinations of the same polarizabilities entering at O(ω 3 ), weighed with an additional power in 1/M . Furthermore, they involve the higher-order spin polarizabilities defined in Ref. [30], with a coefficient in 1/M , and the dispersive coefficients α L 4 and β L 4 corresponding to sixth-order scalar polarizabilities, which have never been defined in literature. Following Ref. [1], we can write them as combinations of the second-order derivatives of the non-Born contribution to the Lorentz invariant amplitudes A i , i.e.
At O(ω 5 ), one finds a recoil contribution in 1/M given by combinations of these 18 constants, which correspond to a combination of the dispersive effects of the sixth-order scalar polarizabilities entering at O(ω 4 ) and new scalar sixth-order polarizabilities, which have never been discussed so-far in literature. These terms are collectively indicated with the α L 5 and β L 5 coefficients in Eqs.
In Fig. 1, we show the predictions for the scalar DDPs from the full DR calculation, without LEX, in comparison with the results obtained from Eqs. (6), using the predictions from DRs for all the static polarizabilities entering the α L0 E1 (ω) and β L0 M1 (ω) contributions and the results from the fit to the DR calculation for the residual functions f α (ω) and f β (ω). We note that the parametrization in Eq. (6) is able to reproduce very well the full energy dependence of the scalar DDPs in the energy range considered in the present fit, giving us confidence that it can be conveniently adopted for our fitting procedure of the scalar DDPs to the Compton scattering data.

III. FITTING STRATEGY
In this section, we outline the fitting strategy for the coefficients of the LEX in Eq. (6). As we are interested into the genuine dispersive effects of the E1 and M 1 radiation in the LEX of the scalar DDPs, we fixed the recoil contributions in α L0 E1 (ω) and β L0 M1 (ω) from the static leading-order spin polarizabilities to the experimental values of Ref. [8], and the recoil terms from higher-order spin polarizabilities as well as from quadrupole scalar polarizabilities to the values predicted from subtracted DRs [30].
As input for the subtracted dispersion integrals we used the updated MAID solution, which is employed also for the calculation of the multipole amplitudes which are not fitted. We used two different data sets, i.e., all the available experimental data for the unpolarized cross sections below pion-production threshold, denoted as "FULL" data set 1 (for a total of 150 data points) [26,[31][32][33][34][36][37][38][39][40][41][42], and the data set given by the TAPS experiment alone (55 data points) [26], which is, by far, the most comprehensive available subset. "Improved" data sets have been defined in recent fits of RCS observables, by discarding some data points from different experiments [4,20,43]. Here we do not apply any selection to the data, and we postpone a more detailed discussion of the statistical consistency of the different data subsets to a future work [44], entirely devoted to the extraction of the static scalar dipole polarizabilities from data. As a first attempt, we tried to fit α E1 , β M1 , α E1,ν , β M1,ν and the four coefficients parametrizing the residual functions f α (ω) and f β (ω), for a total of eight parameters. This choice has the advantage of fitting the full energy dependence of the scalar DDPs with a minimum of model dependence. We then applied the gradient method in the χ 2 minimization procedure of MINUIT [45]. Unfortunately, this method did not show convergence since the positive-definiteness condition of the covariance matrix could not be achieved. This is because of too strong correlations between the fitted parameters, resulting from the very low sensitivity of the available experimental data to the higher order dispersive coefficients (note that the fitting parameters enter to all orders in the LEX of the scalar DDPs, both as genuine dispersive effects and as recoil contributions). To circumvent this problem, we used a combination of the simplex [46] (that is a purely geometric minimization algorithm) and bootstrap (that is a Monte Carlo technique) methods. Each bootstrap "measurement" is assumed to be Gaussian distributed around a given experimental data point with a standard deviation given by its statistical error. All bootstrapped points of a given subset are then shifted by the same quantity proportional to the published systematic error, assumed to be uniformly distributed. If we define as cycle a number of bootstrapped points equal to the total number of points in the considered experimental data set, the bootstrap sampling can be finally described as: where S stands for the differential cross section, with the superscript "exp" indicating the experimental values for the mean value S exp i and the statistical error σ exp i . In Eq. (9), the index i runs over the data points in the whole set, the index k labels each subset, and the index j indicates the bootstrap cycle. Furthermore, ξ k,j are random numbers uniformly distributed as U[1 − ∆ k , 1 + ∆ k ] (with ±∆ k the published systematic error), while the numbers γ i,j are sampled from the standard Gaussian distribution N [0, 1]. When different systematic-error sources are quantified, ξ k,j is the combination of all the contributions. According to Ref. [26], σ exp i for the TAPS data set includes a ±5% point-to-point systematic error added in quadrature to the statistical error of each individual data point.
The minimization is performed after a complete cycle, and the output for the fitted values of the polarizabilities is stored. Repeating the bootstrap cycle a very large n R number of replicas (we choose n R = 10000), we are fi-nally able to reconstruct the probability distributions for every fitted parameter. In order to obtain a cross-check for our fitting method, we first assumed as fit parameters only α E1 − β M1 , using the constraint from the Baldin's sum rule for the polarizability sum, with the TAPS value α E1 + β M1 = (13.8 ± 0.4) · 10 −4 fm 3 [26] 2 , fixing the leading-order spin polarizabilities to the central values of Ref. [8], all the other static polarizabilities as well as the residual functions f α (ω) and f β (ω) to subtracted DRs. This configuration with only one free parameter allows the gradient method to converge, thus providing a benchmark both for our new minimization algorithm and for the theoretical framework. The validation of the method should pass the following tests: • The one-parameter fit with our strategy (the combination of subtracted DRs, LEX, multipole expansion and bootstrap technique, labeled as fit 1) should be consistent with the fit using the complete DR calculation (without multipole expansion and LEX), and the gradient method for the minimization of the χ 2 function (labeled as fit 2); • The results of the fit 1 should be consistent with the results obtained from the combination of subtracted DRs, LEX, multipole expansion and gradient method (labeled as fit 3).
For the purposes of this test, the bootstrap sampling procedure was performed after fixing to unity all ξ k,j parameters in Eq. (9). The results from the three fits using the FULL data set are shown in Table I. They are all consistent with each other and also agree with the PDG values [7]. In addition, the errors evaluated using fit 1 are Gaussian-distributed, in agreement with the statistical expectations (see, for instance, [47]). The bootstrap solution does not significantly change when the uncertainty in the Baldin's sum rule value is taken into account by randomly generating, for each cycle, a different α E1 + β M1 value sampled from the N [13.8, (0.4) 2 ] distribution. Similarly, the results of the fit parameters change at most by 1% when the values of the leading-order spin polarizabilities are varied within the uncertainties quoted in Ref. [8].
All this gives us confidence in both the theoretical framework, based on the LEX and multipole expansion, and the statistical tools based on the bootstrap technique. In summary, the main advantages of the adopted technique are: • The straightforward inclusion of systematic errors in the minimization procedure, as shown in Eq. (9). This feature allows us to reduce the overall number of fit parameters with respect to the extended- χ 2 procedure, with a normalization factor for each data set left as free parameters [4]; • The fact that any error distribution of the experimental data can be easily implemented. Moreover, the probability distributions of the fitted parameters are not assumed a priori to be Gaussian, but are directly evaluated from the probability distributions assigned to the experimental data; • The possibility to automatically take into account the effects of the differential cross section systematics in the error bars of the LEX coefficients of the DDPs.
The application of this strategy to the eight-parameter fit gave probability distributions with very broad and asymmetric tails in all cases, except for α E1 . For instance, the 68% confidence interval of α E1,ν was found to be α E1,ν = (2.81 +8.14 −6.44 ) · 10 −4 fm 5 . This is a further confirmation that the quality of the present experimental database is too poor to allow any meaningful estimate of the higher-order coefficients in the LEX of the scalar DDPs, considering the very low sensitivity to them of the differential cross section below pion-production threshold. For this reason we reduced the number of free parameters to three: the polarizability difference α E1 − β M1 and the dispersive polarizabilities α E1,ν and β M1,ν [29]. The remaining five parameters were fixed using the Baldin's sum rule value, smeared according to its resolution, for the polarizability sum α E1 + β M1 , and DRs for the four parameters in the residual functions f α (ω) and f β (ω). In particular, the coefficients of the residual functions were fixed to describe the full energy-dependence of the scalar DDPs predicted from DRs (see Fig. 1).

IV. RESULTS AND DISCUSSIONS
In this section we discuss the results from the new fitting strategy with three fit parameters. The results for the LEX coefficients of the scalar DDPs in both the cases of the FULL and TAPS data set are shown in Table II. The corresponding probability distributions are still Gaussian distributed, as displayed in Fig. 2.
In the three-parameter fit, also the gradient method showed convergence (without the inclusion of systematic errors), but the covariance matrix was forced to be positive definite. This feature casts doubts on the validity of the procedure even if the fitted parameter estimates turned out to be close to those given in Table II Some comments are in order: • It is likely that the FULL data set includes inconsistent data, but the small number of experimental points in each data subset does not allow us to perform consistency checks without introducing biases; • The central values of the static dipole polarizabilities α E1 and β M1 from the fit with FULL and the TAPS data sets are different, but still compatible within the errors and in fairly good agreement with the PDG values.
This difference could also be due to the correlation between the different angular distributions of the two data sets and the varying sensitivity to α E1 and β M1 in different angular regions. In Fig. 3 we show the two-dimensional joint probability distributions for α E1 , β M1 , α E1,ν and β M1,ν . Apart from the strong correlation between α E1 and β M1 , which is due to the constraint of the Baldin's sum rule, we observe significantly strong correlations also for all the other distributions. This is especially true in the case of the dispersive polarizabilities α E,1ν and β M1,ν , that have very strong negative correlation coefficients with α E1 and β M1 , respectively.
As already noticed before, this behavior is mainly a consequence of low sensitivity of the existing data to the magnetic polarizabilities.
Even if this effect could also partially be due to inconsistent data subsets, as also recently discussed in [20], only a relevant progress in both the quality and the quantity of the existing data set will allow to significantly improve this situation.
Predictions for the dispersive polarizabilities have been obtained within unsubtracted [1] and subtracted [30] DRs, and baryon chiral perturbation theory [16]. Our extraction of the dispersive polarizabilities is, within the large error bars, consistent with these theoretical predictions and can not discriminate between them. We pointed out that our fitting method provides a realistic probability distribution of the fitted parameters; another advantage is that it allows us to straightforwardly compute the error bands for the DDPs including all the correlations between the parameters. In Fig. 4, we show our fit results for the scalar DDPs, extracted from the FULL and TAPS data set, as function of the cm energy ω and with the corresponding 68% and 95% confidence level (CL) uncertainty bands. They are compared with the subtracted DR predictions, obtained with the values of the static dipole polarizabilities from the fit to the FULL and TAPS data sets in Tab. II. The DR results for both the scalar DDPs are within the 68% confidence area of the fit results for ω 60 MeV. At higher energy, the DR predictions for β M1 (ω) remain within the 95% CL region, while for α E1 (ω) we observe deviations from the fit results in the case of the FULL data set and a very good agreement, within the 68% confidence area, in the case of the TAPS data set. This different behavior can be again a hint of inconsistencies between the two data sets. The larger relative error in the case of β M1 (ω) also reflects the lower sensitivity of the unpolarized RCS data to the magnetic polarizability than to the electric polarizability.

V. CONCLUSIONS
In summary, we have presented a new method based on the parametric bootstrap technique that allows us to extract, for the first time, information on the proton scalar DDPs from RCS data at low energies. This method was never exploited so far to analyze Compton scattering data and has several advantages with respect to the standard χ 2 -minimization technique. For example, it allows one to include in a straightforward way the systematic errors in the minimization procedure, without introducing a large number of additional fit parameters, and provides error distributions of the experimental data which are not assumed Gaussian a priori, but are directly evaluated from the probability distributions assigned to the data. The extraction of the energy dependence of the DDPs turned out to be quite challenging, because of the very low sensitivity of the unpolarized RCS data to the higherorder dispersive coefficients. This gives both large error bands of our estimates, in particular for β M1 (ω), and strong correlations between the fit parameters.
In the present analysis, the theoretical framework is based on dispersion relations, but it can be conveniently adapted to use other inputs, such as effective field theory calculations [4,16] that have been recently employed to extract the static polarizabilities.
Finally, future measurements at MAMI [48] hold the promise to improve the accuracy and the statistics of the data and will help to determine with better accuracy the effects of the leading-order static and dynamical polarizabilities.