Solving the PCAC puzzle for nucleon axial and pseudoscalar form factors

It has been observed in multiple lattice determinations of isovector axial and pseudoscalar nucleon form factors, that, despite the fact that the partial conservation of the axialvector current is fulfilled on the level of correlation functions, the corresponding relation for form factors (sometimes called the generalized Goldberger-Treiman relation in the literature) is broken rather badly. In this work we trace this difference back to excited state contributions and propose a new projection method that resolves this problem. We demonstrate the efficacy of this method by computing the axial and pseudoscalar form factors as well as related quantities on ensembles with two flavors of improved Wilson fermions using pion masses down to 150 MeV. To this end, we perform the $z$-expansion with analytically enforced asymptotic behaviour and extrapolate to the physical point.


Introduction
The axial nucleon structure is central for the description of weak interactions and plays a prominent role in long-baseline neutrino experiments, where it is important for a precise determination of the neutrino flux and the cross section for nuclear targets [1][2][3].The nonperturbative information encoded in the corresponding form factors is manifold.For example, the isovector axial form factor is linked to the flavor asymmetry in the difference between helicity aligned and anti-aligned quark densities in impact parameter space [4].While the axial coupling g A is measured quite precisely in β-decay, the nucleon form factors that encode the spatial structure are much less well known.The axial form factor, G A (Q 2 ), and the induced pseudoscalar form factor, GP (Q 2 ), enter the description of (quasi-)elastic neutrino-nucleon scattering [5][6][7][8] and exclusive pion electroproduction [9][10][11][12] (e.g., e − p → π − pν).They can also be measured in muon capture [13][14][15][16].For reviews see, e.g., Refs.[16,17].
Apart from experimental measurements, there are various tools available that can constrain form factors from the theory side.At small virtualities chiral perturbation theory can be used to obtain valuable constraints (see, e.g., Refs.[12,17,18]).At intermediate and large Q 2 the form factors can be calculated (up to some systematic uncertainty of ∼15%) using light-cone sum rules [19,20].Alternatively, one can use functional renormalization group methods [21].
However, the cleanest method for the determination of hadron form factors is lattice QCD.Various determi-nations of the nucleon couplings and form factors using a wide variety of lattice actions and analysis methods can, e.g., be found in Refs. .Recent calculations that have precisely determined the axial, the pseudoscalar, and the induced pseudoscalar form factors separately from lattice data yield an unexpected result: The relation between these form factors inferred from the partial conservation of the axialvector current (denoted the PCAC FF relation in the following) is broken rather badly [46][47][48][49][50]. 1 Adding to the confusion, one should note that the PCAC relation itself is still fulfilled quite well on the level of the correlation functions, which leads to the conclusion that either discretization effects or excited state effects are responsible for the observed discrepancy.However, the former have been ruled out as (the sole) explanation in [46], while it was found in [48] that even a 3-state fit cannot resolve the issue.Let us note in passing that simply enforcing the PCAC relation on the form factor level might lead to uncontrollable systematic effects.
In this work we demonstrate that the largest part of the deviation from the PCAC FF relation is indeed due to excited states in the temporal axialvector (A 0 ) and pseudoscalar (P ) channels.These excited states are, however, so strongly enhanced relative to the ground state that the usual multistate fit ansatz is bound to fail for any feasible time distances between the source, the sink, and the current insertion.While this is directly visible in A 0 correlation functions, which are therefore usually omitted in the analysis when extracting G A (Q 2 ) and GP (Q 2 ), up to Table 1: Details of the ensembles used in the analysis, including the inverse lattice coupling β, the hopping parameter κ, the lattice geometry, the pion mass mπ, and the spatial lattice extent L = aNs in units of m −1 π .The finite volume pion masses were determined in Ref. [51] and the errors include an estimate of both the systematic and statistical uncertainty.The lattice spacings and renormalization factors are listed in Table 2.
Ens.If the reader is now eager to learn the details of the method, he or she should skip the usual description of simulation parameters and analysis methods provided in Section 2, and directly jump to Section 3 where the problem and its resolution will be explained in detail.The results are then presented in Section 4, before we conclude.

Lattice setup
In this work we analyse ensembles with two flavors of nonperturbatively improved clover fermions (also known as Sheikholeslami-Wohlert fermions [52]) and the Wilson gauge action at three different β values corresponding to lattice spacings in the range of 0.060 fm to 0.081 fm.The spacing was set using the Sommer parameter r 0 = 0.5 fm determined in Ref. [53] (see also Ref. [54]).The pion masses range from 490 MeV down to an almost physical value of 150 MeV.The ensemble details are provided in Table 1 and Fig. 1 visualizes the landscape of available pion masses and volumes.The multiple volumes for a pion mass of ∼290 MeV, ranging from Lm π = 3.4 to 6.7, enable an investigation of finite volume effects.
The isovector form factors can be extracted from connected three-point functions.The latter have been computed as part of a previous study of the nucleon isovector charges [32] using the traditional sequential source  Lmπ plotted against the pion mass for our ensembles listed in Table 1.The color coding for the ensembles is used throughout this work.
method [55], where the insertion time t ins of the local current is varied, while the sink and source times, t snk and t src , are fixed.To increase statistics two measurements of the three-point functions are performed at different source positions per ensemble.Auto-correlations between configurations are taken into account by binning with a binsize of 10.
In order to minimize excited state contributions to the three-point functions (and to the two-point functions that are also required), spatially extended source and sink interpolating operators are constructed using Wuppertal smearing [57] with APE-smeared [58] gauge links.In Ref. [32] the smearing was optimized such that ground state dominance was observed in the nucleon two-point functions at around the same physical time, t = t snk − t src ≈ 0.8 fm, for different pion masses and lattice spacings.For key ensembles, labelled III, IV, and VIII, corresponding to pion masses of 420, 290, and 150 MeV, respectively, at the lattice spacing a = 0.071 fm, multiple source-sink separations for the three-point functions were generated to enable an investigation of remaining excited state contamination.The separations correspond to t ≈ (1.1, 1.2) fm for ensemble III, t ≈ (0.5, 0.6, 0.8, 0.9, 1.1, 1.2) fm for ensemble IV and t ≈ (0.6, 0.9, 1.1) fm for ensemble VIII.The analysis of the isovector charges indicated that, for the smearing applied, the ground state contribution could be reliably determined for t ≳ 1 fm and this (single) sourcesink separation was employed for the remaining ensembles.
The simultaneous fits to the two-and three-point functions for the present analysis, taking into account the leading excited state contribution, are discussed in Section 2.3.

Correlation functions and form factor decompositions
We analyse the two-and three-point functions where t = t snk − t src , τ = t ins − t src , O N = (u T Cγ 5 d)u is the Wuppertal-smeared interpolating current for the nucleon with the charge conjugation matrix C, and P + = (1 + γ 0 ) 2 projects onto positive parity for zero momentum.We choose Γ to be P i + = P + γ i γ 5 , i = 1, 2, 3, in our analysis.In our actual simulations, we restrict the kinematics to Inserting complete sets of states, the correlation functions in Euclidean time can be expanded in terms of hadronic matrix elements.For the two-point function this yields where only the ground state contribution is given and m N denotes the nucleon mass.The excited state corrections are discussed in Sect.2.3.The normalization Z ⃗ p is smearing-dependent and encodes the overlap of the ground state with the interpolating operators at the source and the sink, where u β ⃗ p,σ is a nucleon spinor.An analogous spectral decomposition of the three-point functions with an operator insertion O ∈ {P, A µ } corresponding to the isovector pseudoscalar and axialvector currents, leads to with J[O] is defined by the form factor decomposition For the different channels it reads where q = p ′ − p and  The axial Ward identity yields a partial conservation of the axialvector current, ∂ µ A µ = 2i m q P , known as the PCAC relation.On the lattice this relation can be broken by discretization effects.For the nucleon matrix elements it implies: Using the definitions ( 9) and ( 10) together with the equations of motion one can deduce the corresponding relation for the form factors (PCAC FF ): Eqs. ( 11) and ( 12) should be satisfied to a similar degree once the ground state matrix elements have been extracted reliably.

Excited states analysis
In three-point functions the signal-to-noise ratio decreases exponentially with the time distance between the source and the sink.Hence, within typical separations t ≲ 1.5 fm a sufficient suppression of excited states may not be achieved.Including the leading excited state contributions to Eqs. ( 3) and ( 6) gives: where ∆E ⃗ p denotes the energy gap between the first excited state and the ground state.The excited state coefficient Z depends on the the nucleon interpolator, its smearing, and the momenta, while B 10 , B 01 , and B 11 also depend on the current O and on the projector Γ.  20)) utilizing the PCAC relation (19) for the ensembles VI (diamonds) and VIII (circles).The result is given for various initial momenta ⃗ p in units of 2π L, while the final momentum is always fixed to ⃗ p ′ = ⃗ 0 in our kinematics.Indeed, the PCAC relation is valid on the three-point function level, up to O(a 2 ) effects.
For illustrative purposes we define the ratio which eliminates the leading order time dependence as well as the overlap factors.Our ground state energies are well described by the continuum dispersion relation as shown in Fig. 2 and we assume this functional form in our fitting analysis.Excited states will in general contain more than one hadron and hence we make no such assumption for ∆E ⃗ p .We remark that the spectrum includes N π, N ππ, and higher states and that this spectrum becomes more dense as the pion mass decreases.At zero momentum the lowest multiparticle excited states are P -wave N π and S-wave N ππ.
For a given momentum transfer q 2 , a simultaneous fit of the form of Eqs. ( 13) and ( 14) is performed to the relevant two-and three-point functions, including all available momentum directions, hadron polarizations, and sourcesink separations.The form factors enter the fit directly as parameters, substituting the amplitudes B ⃗ p ′ , ⃗ p Γ,O utilising Eqs. ( 7)- (10) and the dispersion relation.For ensembles with only one value of t, the parameter B 11 is set to zero.The fit range is chosen to be 2a ≤ τ ≤ t − 2a resulting in reasonable values of χ 2 d.o.f..

Renormalization and O(a) improvement
The renormalization factors have been calculated nonperturbatively in an RI ′ -MOM scheme using the Rome-Southampton method [59] and then converted into the MS scheme using three-loop continuum perturbation theory.A detailed discussion can be found in [56].
The isovector currents are multiplicatively renormalized using where the relevant Z X , listed in Table 2, contain both the nonperturbative renormalization and the conversion to the MS scheme at a scale of µ = 2 GeV.The one-loop improvement coefficients b X have been calculated perturbatively in Refs.[60,61] and are close to unity.The numeric values for our lattices are provided in Ref. [32].Within the Symanzik improvement program [62,63] also the currents themselves have to be O(a)-improved.For the axialvector current this yields where we used the improvement coefficient c A , nonperturbatively determined in Ref. [64], and ∂ µ denotes the symmetrically discretized derivative.Note that for the pseudoscalar current P imp = P .

Quark mass from nucleon correlators
Since the partial conservation of the axialvector current, is an operator relation, it has to hold on the correlation function level such that, using the three-point functions defined by Eq. ( 2), the PCAC quark mass can be determined as independent of any spectral analysis.As no significant O(a 2 Q 2 ) effects are observed within the statistical error, we determine the quark mass by fitting to several Q 2 simultaneously (see Fig. 3).Note that the spatial derivatives can be calculated using the formula as long as one considers only states ⃗ p⟩ for which the external three-momenta have been fixed via an appropriate Fourier transform in the correlation function (2).The time derivative has to be calculated explicitly, since the energy in the correlation function is not fixed.
Usually the PCAC mass is obtained from ratios of pion two-point functions, where one can achieve very small statistical errors.Up to discretization effects, one would expect these values to agree with those obtained from ratios of baryonic three-point functions via Eq.( 20).This is indeed the case, as depicted in Fig. 4 using ensembles I, IV, and XI (which have different lattice spacings, but very similar pion masses and volumes).Note that despite the unphysical pion mass the extrapolated ratio compares reasonably well with the N f = 2 value in the physical limit: m q m 2 π = 0.198(11) GeV −1 [65].

Uncovering the ground state contribution
In a number of lattice simulations it has been observed that, even though the PCAC relation is fulfilled on the correlation function level up to discretization effects (in our case of O(a 2 ), cf.Figs. 3 and 4), the equivalent equation for the nucleon form factors, Eq. ( 12), seems to be broken rather badly.Note that this problem cannot be solved just by using the PCAC mass obtained from the ratio of nucleon three-point functions discussed above.
In the following we will demonstrate that the breaking of the PCAC relation on the form factor level is a consequence of very large excited state effects that cannot be resolved by the standard excited state analysis described in Section 2.3.Using Eq. ( 10) together with ( p−m N )u ⃗ p,σ = 0 one finds for the nucleon ground state: where p = 1 2 (p ′ + p).From Fig. 5 one can easily verify that this equation is violated by the data, demonstrating the presence of large excited state contaminations.This is mostly due to the A 0 correlation function, which has an almost linear dependence on the insertion time, cf. the left panel of Fig. 6.
One possible (and up to date the most widely-used) workaround to this problem is the exclusion of A 0 from the analysis. 2A more satisfactory approach is to consider the current which fulfills p µ A ⊥ µ = 0 by construction.Due to Eq. ( 22) one can be sure that the subtraction only removes excited state contributions, while leaving the ground state contribution unchanged.For the new current the three-point function has the expected behaviour, as shown in Fig. 6.
The connection between the axial and pseudoscalar channels via the PCAC relation suggests that similar excited state contributions are present in the pseudoscalar channel and it is advantageous to construct the combination 2 For comparison we will show results obained using this method in Figs.7-9 Figure 5: Ratio of correlation functions (cf.Eq. ( 15)) corresponding to the l.h.s. of Eq. ( 22) on ensemble VIII, showing the dominant excited state effects.
such that ∂ µ A ⊥ µ = 2i m q P ⊥ .Again, Eq. ( 22) ensures that the ground state matrix element is not affected by this subtraction.We use the PCAC mass m q obtained from baryon three-point functions, cf. the discussion in Section 3.1, which we found leads to smaller discretization effects compared to employing the PCAC mass from pion two-point functions.The effect of Eq. ( 24) is illustrated in the right panel of Fig. 6: While the original data looked less conspicuous than for the A 0 channel, the resulting shift is very significant.
To conclude this section, we remark that the subtraction constructed above should not be viewed as an operator improvement, as it depends on the external momentum.Instead, one should interpret it as a method to systematically construct combinations of correlation functions that suffer less from excited states.In principle, this method can be used wherever an exact relation for the   15)) for the axialvector (left, middle) and pseudoscalar currents (right) with three different source sink separations t a = 9 (blue), 12 (green), 15 (red) for ensemble VIII with mπ ≈ 150 MeV.In both cases we compare results from the standard current (open symbols) and our modified current (filled symbols).The shown results correspond to an average over all relevant combinations of polarizations and momenta with ⃗ p = 2π L and ⃗ p ′ = ⃗ 0 (i.e., Q 2 = 0.073 GeV 2 ).For A 0 the problem is clearly visible.The signal for A ⊥ 0 has a significantly smaller statistical error and only shows mild excited state contributions (middle, zoomed), which are resolvable with the multiexponential ansatz given in Eqs. ( 13) and (14).In contrast, the extent of the excited state contaminations to the data for the pseudoscalar current P is not so obvious.However, subtracting the same excited states causing the problem in the axialvector channel (by using P ⊥ ), one finds that the true ground state plateau lies much higher.The yellow bands indicate the ground state contributions extracted from the fits for O ∈ {A ⊥ 0 , P, P ⊥ }. ground state matrix elements exists (in this case Eq. ( 22)).However, the analogous constraint for the vector current, is fulfilled almost exactly by the data such that the method does not lead to an improvement.

Restoration of PCAC on the form factor level
We define the ratio (cf.also Ref. [46]) where deviations from r PCAC = 1 quantify the violation of the PCAC FF relation (12).Fig. 7 demonstrates that using the method described in Sect.3.2 all ensembles, in particular the ones with small pion mass that previously exhibited the largest deviations, now fulfill the PCAC FF 0.2 0.4 0.6 0.8 1.0 relation reasonably well.Even for large Q 2 ≈ 1 GeV 2 we see a significant improvement, although small deviations of ∼5% remain.This residual violation can be attributed to O(a 2 ) discretization effects of Eq. (19).
In absence of better information, the induced pseudoscalar form factor is often estimated by usually called the pion pole dominance (PPD) assumption.Fig. 8 demonstrates that this approximation does not describe the data, especially not at small Q 2 .This is true for both, the original and the improved data (which reaffirms the findings of Refs.[32,46,49]).Hence, the observed disagreement with the PPD ansatz is certainly not caused by the same excited state effects that have been responsible Note that G P is scale-dependent; here it is plotted for the MS scale 2 GeV.In the upper panels the results using the excited state subtraction explained in Section 3.2 (filled symbols) are compared to those obtained without this method (open symbols).The largest effect is found in the pseudoscalar form factor, while the others are almost not affected.In the lower panels the yellow band shows the final result for the form factors extrapolated to the physical point using the z-expansion described in Section 4.2.In these plots the transparent points show the subtracted data itself, while the solid ones are parallel transported to the physical point in order to give an intuitive grasp on how good these fits actually describe the data points.The color coding follows Fig. 1.
for the PCAC FF violation.Since all data for r PPD collapse onto an almost universal function of m 2 π +Q 2 it seems highly unlikely that the deviation is due to discretization, volume, or quark mass effects.

Parametrization of the form factors
We parametrize the form factors using the z-expansion [66,67], which automatically imposes analyticity constraints.This corresponds to an expansion of the form factors in the variable where t cut = 9m 2 π is the particle production threshold and t 0 is a tunable parameter. 3Isolating the pion pole in the pseudoscalar channels (cf.Ref. [43]) one obtains To enforce the correct scaling in the asymptotic limit, 68], one has to implement the four constraints 3 Varying t 0 between 0 and tcut 2 has no significant impact on the result.Therefore we have simply set it to zero in our analysis.
which can be incorporated by fixing the first four coefficients according to for k = 0, 1, 2, 3, such that one is left with N − 3 free coefficients.In the following we will show the fit results with 3 free coefficients (called a z 3+4 fit in the literature), since the z 2+4 fits failed to describe the data at low momentum transfer, while z 4+4 fits did not yield further improvement.
In order to extrapolate to the physical point (m π → m phys π , a → 0, L → ∞), we use the parametrization a This allows us to perform a combined fit to all ensembles with 15 fit parameters for each form factor.

Form factors, charges, and radii
Our results for the form factors are shown in Fig. 9, where the band shows the value extrapolated to the physical point using the combined fit to all ensembles described in the previous section.Since we know that the PCAC FF relation is fulfilled, we use G P (0) = G A (0)m N m q to obtain a data point for G P in the forward limit (for each ensemble) in order to stabilize the continuum extrapolation of this form factor.In Table 3 we list the z-expansion coefficients corresponding to our central values at the physical point.From the form factors we obtain the charges and the corresponding mean squared radii r 2 = −6G ′ (0) G(0) as The value of g A is in perfect agreement with the experimental result g A g V = 1.2724( 23) [69].Our relatively large result for the squared axial radius still agrees within errors with z-expansion fits to experimental νd [8] and muon capture [16] data, but lies higher than other lattice results in the range 0.2 − 0.4 fm 2 [43][44][45][46], cf.Fig. 7 in Ref. [16].
For the induced pseudoscalar coupling defined at the muon capture point [17] we obtain gP = m µ 2m N GP (0.88m 2 µ ) = 2.8( 7) , where m µ = 105.6MeV is the muon mass.Note that the small value obtained for gP is consistent with the strong violation of the pion pole dominance assumption at small momentum transfer shown in Fig. 8.The experimental value gP = 8.06(48)(28) from muon capture [15] is consistent with PPD, but not with our data.

Summary
We have presented a method to identify (and subtract) excited state contributions that spoil the PCAC relation on the form factor level.This mainly affects correlation functions involving the P and A 0 currents, which have much larger coupling to the pion at small momentun transfer than A i .After our subtraction, PCAC FF is fulfilled up to small deviations at large momentum transfer, which can be interpreted as lattice artifacts.In spite of this improvement, we find that the pion pole dominance assumption that relates the axial and the induced pseudoscalar form factor is still strongly broken at small momentum transfer, which confirms the findings of Refs.[32,46,49].A recent calculation within chiral perturbation theory [70] (along the lines of Refs.[71,72] using interpolating currents from Ref. [73], cf. also Refs.[74,75]) indicates that this deviation at small momentum transfer might be due to additional large excited state contributions in the induced pseudoscalar form factor, which is unaffected by the subtraction method described in Section 3.2.However, our data do not indicate any presence of these excited states.
We parametrize the form factors using the z-expansion and extrapolate them to the physical point by a combined fit to all ensembles at hand.We find that these fits provide a very reasonable description of the data and we use them to extract the values of the charges g A and g P .The result for the axial form factor exhibits a rather steep slope, which is not related to our improvement, but a consequence of the applied z 3+4 parametrization.It corresponds to an axial dipole mass M A = √ 12 r A = 0.77 (8) GeV that is in rough agreement with phenomenological values around 1 GeV.The disagreement of our result for the induced pseudoscalar charge gP with experiment could point to persistent problems in the form factor GP and warrants further investigation.All our physical results (extrapolated to physical quark masses, to infinite volume, and to the continuum) currently come with relatively large errors that are mainly caused by the narrow range of available lattice spacings.Reducing the systematic error at the physical point (e.g., by the analysis of CLS ensembles with finer lattice spacings [76,77]) will be of high priority in the future.
Figure1: Lmπ plotted against the pion mass for our ensembles listed in Table1.The color coding for the ensembles is used throughout this work.

Figure 2 :
Figure 2: Lattice data of the ground state nucleon energy, normalized to the continuum expectation (16).The color coding follows Fig. 1.

Figure 3 :
Figure 3: Quark mass obtained from the ratio (see Eq. (20)) utilizing the PCAC relation(19) for the ensembles VI (diamonds) and VIII (circles).The result is given for various initial momenta ⃗ p in units of 2π L, while the final momentum is always fixed to ⃗ p ′ = ⃗ 0 in our kinematics.Indeed, the PCAC relation is valid on the three-point function level, up to O(a 2 ) effects.

259. 5 Figure 4 :
Figure 4: The continuum extrapolations of the PCAC mass determined from pion two-point functions (orange band) and nucleon three-point functions (green band) are consistent within the statistical errors.Both determinations show the leading quadratic behaviour characteristic for O(a)-improved currents.Dividing the quark mass by the pion mass squared accounts for the leading pion mass dependence given by the Gell-Mann-Oakes-Renner relation.

Figure 6 :
Figure 6: Correlation function ratios (cf.Eq. (15)) for the axialvector (left, middle) and pseudoscalar currents (right) with three different source sink separations t a = 9 (blue), 12 (green), 15 (red) for ensemble VIII with mπ ≈ 150 MeV.In both cases we compare results from the standard current (open symbols) and our modified current (filled symbols).The shown results correspond to an average over all relevant combinations of polarizations and momenta with ⃗ p = 2π L and ⃗ p ′ = ⃗ 0 (i.e., Q 2 = 0.073 GeV 2 ).For A 0 the problem is clearly visible.The signal for A ⊥ 0 has a significantly smaller statistical error and only shows mild excited state contributions (middle, zoomed), which are resolvable with the multiexponential ansatz given in Eqs.(13) and (14).In contrast, the extent of the excited state contaminations to the data for the pseudoscalar current P is not so obvious.However, subtracting the same excited states causing the problem in the axialvector channel (by using P ⊥ ), one finds that the true ground state plateau lies much higher.The yellow bands indicate the ground state contributions extracted from the fits for O ∈ {A ⊥ 0 , P, P ⊥ }.

Figure 8 :
Figure 8: Violation of the pion pole dominance ansatz (26) for the induced pseudoscalar form factor.The filled and open data points are with and without the excited state subtraction described in Sect.3.2, respectively, and agree within errors.The color coding follows Fig. 1.

Q 2 [GeV 2 ]Figure 9 :
Figure9: Results for the form factors.Note that G P is scale-dependent; here it is plotted for the MS scale 2 GeV.In the upper panels the results using the excited state subtraction explained in Section 3.2 (filled symbols) are compared to those obtained without this method (open symbols).The largest effect is found in the pseudoscalar form factor, while the others are almost not affected.In the lower panels the yellow band shows the final result for the form factors extrapolated to the physical point using the z-expansion described in Section 4.2.In these plots the transparent points show the subtracted data itself, while the solid ones are parallel transported to the physical point in order to give an intuitive grasp on how good these fits actually describe the data points.The color coding follows Fig.1.