Nucleon axial and pseudoscalar form factors using twisted-mass fermion ensembles at the physical point

We compute the nucleon axial and pseudoscalar form factors using three $N_f=$2+1+1 twisted mass fermion ensembles with all quark masses tuned to approximately their physical values. The values of the lattice spacings of these three physical point ensembles are 0.080 fm, 0.068 fm, and 0.057 fm, and spatial sizes 5.1 fm, 5.44 fm, and 5.47 fm, respectively, yielding $m_\pi L$>3.6. Convergence to the ground state matrix elements is assessed using multi-state fits. We study the momentum dependence of the three form factors and check the partially conserved axial-vector current (PCAC) hypothesis and the pion pole dominance (PPD). We show that in the continuum limit, the PCAC and PPD relations are satisfied. We also show that the Goldberger-Treimann relation is approximately fulfilled and determine the Goldberger-Treiman discrepancy. We find for the nucleon axial charge $g_A$=1.245(28)(14), for the axial radius $\langle r^2_A \rangle$=0.339(48)(06) fm$^2$, for the pion-nucleon coupling constant $g_{\pi NN} \equiv \lim_{Q^2 \rightarrow -m_\pi^2} G_{\pi NN}(Q^2)$=13.25(67)(69) and for $G_P(0.88m_{\mu}^2)\equiv g_P^*$=8.99(39)(49).


I. INTRODUCTION
The nucleon axial form factors are important quantities for weak interactions, neutrino scattering, and parity violation experiments.There are currently a number of neutrino scattering experiments that require knowledge of the axial form factors.At Fermi Lab, the two neutrino experiments, NOνA and MINERνA [1], share the same neutrino beam.The former is designed to study neutrino oscillations and the latter to perform high-precision measurements of neutrino interactions on a wide variety of materials, including helium, carbon, iron, and lead.The MicroBooNE experiment, also at Fermi Lab, aims at measuring low-energy neutrino cross sections, investigating the low-energy excess events observed by the Mini-BooNE experiment, and studying neutrinos produced in supernovae.The T2K experiment at KEK in Japan and the CNGS experiment in Europe investigate neutrino flavor changes.The upcoming experiment DUNE will be the next-generation flagship experiment on neutrino physics.
These experimental efforts need to be matched by theoretical investigations.Computing reliably the nucleon axial form factors provides crucial input for these experiments.However, the theoretical extraction of these form factors is difficult due to their non-perturbative nature.Phenomenological approaches include chiral perturbation theory that provides a non-perturbative framework suitable for low values of Q 2 up to about 0.4 GeV 2 [2][3][4].Other models used include the perturbative chiral quark model [5], the chiral constituent quark model [6] and light-cone sum rules [7].Lattice QCD provides the ab initio non-perturbative framework for computing such quantities using directly the QCD Lagrangian.Early studies of nucleon axial form factors were carried out within the quenched approximation [8,9], as well as, using dynamical fermion simulations at heavier than physical pion masses [10].Only recently, several groups are computing the axial form factors including simulations generated directly at the physical value of the pion mass [11][12][13][14][15][16][17][18][19][20][21][22][23].This work is the first to use solely simulations performed at physical values of the pion mass to take the continuum limit, avoiding a chiral extrapolation.
The nucleon matrix element of the isovector axialvector current A µ is written in terms of two form factors, the axial, G A (Q 2 ), and the induced pseudoscalar, G P (Q 2 ).The axial form factor, G A (Q 2 ), is experimentally determined from elastic scattering of neutrinos with protons, ν µ +p → µ + +n [24][25][26], while G P (Q 2 ) from the longitudinal cross-section in pion electro-production [4,27,28].At zero momentum transfer the axial form factor gives the axial charge g A ≡ G A (0), which is measured in high precision from β-decay experiments [29][30][31][32].The induced pseudoscalar coupling g * P can be determined via the ordinary muon capture process µ − + p → n + ν µ from the singlet state of the muonic hydrogen atom at the muon capture point, which corresponds to momentum transfer squared of Q 2 = 0.88m 2 µ [33][34][35][36][37], where m µ is the muon mass.We also study the nucleon matrix element of the isovector pseudoscalar current that determines the pseudoscalar form factor G 5 (Q 2 ) and from it the pion-nucleon coupling constant g πN N .
In this work, we use three ensembles generated at physical quark masses of the light, strange, and charm quarks and at three values of the lattice spacing, namely a = 0.080 fm, a = 0.068 fm, and a = 0.057 fm.This same setup has been used in the calculation of the electromagnetic form factors [38] and transversity form factors [39].This allows us to directly take the continuum limit of the axial and pseudoscalar form factors using, for the first time, only simulations performed at the physical pion mass.This is a major achievement since it avoids chiral extrapolation which, for the baryon sector, may introduce an uncontrolled systematic error.Such simulations at the physical pion mass can be used to check important relations, such as the partially conserved axial-vector current (PCAC) relation that at form factor level connects G A (Q 2 ) and G P (Q 2 ) with G 5 (Q 2 ).At low Q 2 and assuming pion pole dominance (PPD) one can further relate G A (Q 2 ) to G P (Q 2 ) and derive the Goldberger-Treiman relation.These relations have been studied within lattice QCD and will be discussed in detail in this paper.
The remainder of this paper is organized as follows: In section II we discuss the decomposition of the nucleon matrix elements of the axial-vector and pseudoscalar operators in terms of form factors and the PCAC and Goldberger-Treiman relations and the pion pole dominance.In section III we give the details on the parameters of the twisted mass fermion ensembles analyzed and in section IV we discuss the extraction of the form factors from the two-and three-point correlators including the renormalization procedure.In section V we present the methods we employ for the identification of excited states and the extraction of the ground state matrix element, as well as the various fits we perform and the model averaging procedure.In section VI, we discuss our procedure of fitting the q 2 -dependence of the form factors and taking the continuum limit, and in section VII, we give the results on the axial form factor, G A (Q 2 ), in the continuum limit.In section VIII we present the analogous analysis for the induced pseudoscalar, G P (Q 2 ), and pseudoscalar, G 5 (Q 2 ), form factors.We also investigate the PCAC and Goldberger-Treiman (GT) relations and evaluate the GT discrepancy.In section IX we compare with other recent lattice QCD results and in section X we summarize and provide our conclusions.In the appendix A, we provide values and parameterization of form factors at the continuum limit.

II. DECOMPOSITION OF THE NUCLEON AXIAL-VECTOR AND PSEUDOSCALAR MATRIX ELEMENTS
In this work, we consider only isovector quantities and neglect isospin-breaking effects due to QED interactions and u-d quark mass difference.Any corrections arising from such isospin-breaking effects are in fact immaterial as compared to our present accuracy and are expected to become relevant only at better than one percent precision.We summarize here for completeness the various relations using the same notation as that used in our previous work [19].The isovector axial-vector operator is given by where u and d are the up and down quark fields respectively.In the chiral limit, where the pion mass m π = 0, the axial-vector current is conserved, namely ∂ µ A µ = 0.For a non-zero pion mass, the spontaneous breaking of chiral symmetry relates the axial-vector current to the pion field ψ π , through the relation We use the convention F π = 92.9MeV for the pion decay constant.In QCD, the axial Ward-Takahashi identity leads to the partial conservation of the axial-vector current (PCAC) where P is the pseudoscalar operator and m q = m u = m d is the light quark mass for degenerate up and down quarks.Using the PCAC relation, it then follows that the pion field can be expressed as The nucleon matrix element of the axial-vector current of Eq. ( 1) can be written in terms of the axial, G A (Q 2 ), and induced pseudoscalar, G P (Q 2 ), form factors as where u N is the nucleon spinor with initial (final) 4momentum p (p ′ ) and spin s (s ′ ), q = p ′ − p the momentum transfer, q 2 = −Q 2 and m N the nucleon mass.
The axial form factor is commonly parameterized as where are the axial charge and radius, respectively.A quantity of interest for the induced pseudoscalar form factor is the induced pseudoscalar coupling determined at the muon capture point [40], namely with m µ = 105.6 MeV the muon mass.It was computed in chiral perturbation theory in Ref. [41].
The nucleon pseudoscalar matrix element is given by where P = ūγ 5 u − dγ 5 d is the isovector pseudoscalar current.The PCAC relation at the form factors level relates the axial and induced pseudoscalar form factors to the pseudoscalar form factor via the relation Making use of Eq. ( 4), one can connect the pseudoscalar form factor to the pion-nucleon form factor G πN N (Q 2 ) as follows Eq. ( 12) is written so that it illustrates the pole structure of G 5 (Q 2 ) and the preferred usage of m q G 5 (Q 2 ), which is a scale-independent quantity unlike G 5 (Q 2 ).Substituting m q G 5 (Q 2 ) in Eq. (11), one obtains the Goldberger-Treiman relation [10,42] (13) The pion-nucleon form factor G πN N (Q 2 ) at the pion pole gives the pion-nucleon coupling which can be computed using Eq. ( 12) to obtain lim Equivalently, g πN N can be computed using Eq. ( 13), where the pole on the right-hand side of Eq. ( 13) must be compensated by a similar pole in Additionally, close to the pole, the following relation holds due to pion pole dominance (PPD).Inserting it in Eq. ( 12) we obtain the relation which relates in Eq. ( 13) we obtain the well-known relation [43] m which means that G P (Q 2 ) can be expressed as [44] G close to the pion pole.From Eq. ( 19), the pion-nucleon coupling can be expressed as where the latter holds in the chiral limit, m π = 0.The deviation from Eq. ( 21) due to the finite pion mass is known as the Goldberger-Treiman discrepancy, namely and it is estimated to be at the 2% level [45] in chiral perturbation theory.The Goldberger-Treiman discrepancy is related to the low-energy constant d18 [46,47] via Given the above relations, we define the following ratios to test whether our lattice results satisfy these relations.
The first is based on the PCAC relation in Eq. (11).Since PCAC is an exact operator relation, it provides a stringent test of our analysis on the form factor level.The second and third relations assume pion pole dominance and use Eqs.(20) and (18), respectively, and they are only expected to be unity near the pion pole.We note that we can use the PCAC relation in Eq. ( 11) to write Using the parameterization of G A (Q 2 ) in Eq. ( 6) to evaluate G A (−m 2 π ) we obtain that near the pion pole the ratio  I. Parameters for the N f = 2 + 1 + 1 ensembles analyzed in this work.In the first column, we give the name of the ensemble, in the second the lattice volume, in the third β = 6/g 2 with g the bare coupling constant, in the fourth the lattice spacing, in the fifth the pion mass, and in the sixth the value of mπL.Lattice spacings and pion masses are taken from Ref. [54].
at leading order in m 2 π , ∆ GT and Q 2 .Using the latter in Eq. ( 27) we obtain [48] and therefore a deviation from unity in r PPD,2 (Q 2 ) can be related to the Goldberger-Treiman discrepancy.

III. GAUGE ENSEMBLES AND STATISTICS
We employ the twisted-mass fermion discretization scheme [49,50], which provides automatic O(a)improvement [51].The bare light quark parameter µ l is tuned to reproduce the isosymmetric pion mass m π = 135 MeV [52,53], while the heavy quark parameters, µ s and µ c are tuned using the kaon mass and an appropriately defined ratio between the kaon and D-meson masses as well as the D-meson mass, following the procedure of Refs.[52,53].The action also includes a clover term that reduces isospin-breaking effects due to the twisted-mass fermion discretization.The values of the parameters of the ensembles analyzed in this work can be found in Table I.The lattice spacings and pion masses are taken from Ref. [54].The values of the lattice spacing are determined both in the meson and nucleon sectors.We quote the ones from the meson sector which are compatible with the values determined from the nucleon mass in Ref. [55].
The resulting tuned pion masses, given in Table I, deviate by up to 4% from the isosymmetric pion mass.This deviation is comparable with the mass difference between charged and neutral pion.Thus, we expect any correction on the form factors arising from such a deviation to be of the same order of magnitude as isospin-breaking effects and, thus, immaterial as compared to our present accuracy.
The nucleon matrix elements of the axial-vector and pseudoscalar operators are obtained via appropriate combinations of three-and two-point nucleon correlation functions, as will be explained in more detail in the following section.In Table II, we give the statistics used for computing the two-and three-point functions in terms of the number of configurations analyzed and the number of point sources employed per configuration.The statis-  In each table, we provide the sink-source separations used in lattice units (first column) and physical units (second column) and the number of source positions per configuration (third column).For each ensemble, the bottom row indicates the number of source positions used for the twopoint functions.
tics of the three-point functions are increased at increasing source-sink separation such that the errors are kept approximately constant among all the time separations.For the twisted mass formulation employed here, the disconnected quark loop contributions are of order a 2 and, thus, vanish in the continuum limit [49].For this reason, we can safely neglect them in the present work.

IV. EXTRACTION OF NUCLEON MATRIX ELEMENTS
To evaluate the nucleon matrix elements of the operators given in Eqs. ( 5) and (10), we compute three-and two-point correlation functions.The two-point function is given by where x 0 is the source, x s is the sink positions on the lattice, and Γ 0 is the unpolarized positive parity projector Γ 0 = 1 2 (1 + γ 0 ).States with the quantum numbers of the nucleon are created and destroyed by the interpolating field where C is the charge conjugation matrix.By inserting the unit operator in Eq. ( 29) in the form of a sum over states of the QCD Hamiltonian only states with the quantum numbers of the nucleon survive.The overlaps between the interpolating field and the nucleon state |N ⟩,  such as ⟨Ω|J N |N ⟩, need to be canceled to access the matrix element.It is desirable to increase the overlap with the nucleon state and reduce it with excited states so that the ground state dominates for as small as possible Euclidean time separations.This is because the signalto-noise ratio decays exponentially with the Euclidean time evolution.To accomplish ground state dominance, we apply Gaussian smearing [56,57] to the quark fields entering the interpolating field where the hopping matrix is given by The parameters α G and N G are tuned [58,59] in order to approximately give a smearing radius for the nucleon of 0.5 fm.For the links entering the hopping matrix, we apply APE smearing [60] to reduce statistical errors due to ultraviolet fluctuations.In Table III, we give the APE and Gaussian smearing parameters used for each ensemble.
For the construction of the three-point correlation function, the current is inserted at time slice t ins between the time of the creation and annihilation of the states with the nucleon quantum numbers, t 0 and t s , respectively.The expression for the three-point function is given by where Γ k = iΓ 0 γ 5 γ k and j µ is either the axial-vector current A µ needed for computing the matrix elements in Eq. ( 5) or P for computing the pseudoscalar form factor in Eq. (10).The Euclidean momentum transfer squared is given by Q 2 = −q 2 = −(p ′ − p) 2 .The connected three-point functions are computed using sequential propagators inverted through the sink, i.e. using the so-called fixed-sink method.This requires new sequential inversions for each sink momentum.Therefore, we restrict to ⃗ p ′ = 0, meaning the source momentum ⃗ p is determined via momentum conservation by the momentum transfer as ⃗ p = −⃗ q and in the following we drop the usage of ⃗ p ′ .Without loss of generality, we also take, in the following, t s and t ins relative to the source time t 0 , or equivalently t 0 is set to zero.
A. Excited states contamination and large time limit The interpolating field in Eq. ( 30) creates a tower of states with the quantum numbers of the nucleon.The spectral decomposition of the two-and three-point functions are given respectively by A i,j µ (Γ k , ⃗ q)e −Ei( ⃗ 0)(ts−tins)−Ej (⃗ q)tins . ( The coefficients of the exponential terms in the two-point function of Eq. ( 34) are overlap terms given by where spin indices are suppressed.The i-index denotes the i th state with the quantum numbers of the nucleon that may also include multi-particle states.The coefficients A i,j appearing in the three-point function of Eq. (35) are given by where ⟨N i ( ⃗ 0)|A µ |N j (⃗ p)⟩ is the matrix element between i th and j th states.In practice, one truncates the sums in Eqs.(34) and (35) up to some state N st .Finally, E i (⃗ p) is the energy of the i th state carrying momentum ⃗ p.For the ground state, we use the dispersion relation to obtain where the nucleon mass m N is determined from the zero momentum projected two-point function.
To extract the nucleon matrix element that we are interested in, any contribution from nucleon excited states and/or multi-particle states has to be sufficiently suppressed.How fast ground state dominance is achieved, depends on the smearing procedure applied on the interpolating fields and the current type entering the threepoint function.Since the noise increases exponentially with increasing t s , establishing from the data convergence to the asymptotic ground state matrix element is very difficult.For this reason, we employ a multi-state analysis by fitting the explicit contribution of the first N st − 1 excited states.Our fitting strategy is described in Sec.V and aims at determining reliably the values of c 0 and A 0,0 µ .In order to cancel unknown overlaps of the interpolating field in Eq. ( 30) with the nucleon state, one commonly constructs an appropriate ratio of three-to a combination of two-point functions [61][62][63][64], The ratio in Eq. ( 39) is constructed such that it converges to the nucleon ground state matrix element in the limit of large time separations ∆E(t s − t ins ) ≫ 1 and ∆E t ins ≫ 1, where ∆E is the energy difference between the first excited state and ground state, namely By substituting Eqs.(34) and (35) into Eq.( 39) we obtain The ground-state matrix elements are extracted from the ratio of the amplitude A 0,0 µ (Γ k , ⃗ q) of the three-point function and the amplitudes c 0 ( ⃗ 0) and c 0 (⃗ q) of the two-point functions.In this work, we determine these amplitudes from a simultaneous fit to two-and three-point functions as described in Sec.V, so that the energy spectrum is the same for two-and three-point functions.For visualization purposes only, we use the following variant of the ratio in Eq. (39), which has the same large time-separation limit as the ratio of Eq. (39) when t ins = t s /2 while avoiding potential excited state contaminations in the two-point functions for small values of t ins .

B. Analysis of nucleon correlators
The ground state matrix elements, Π µ , are decomposed into form factors.In the following, we provide their decomposition in Euclidean space and for ⃗ p ′ = 0.In the case of the matrix element of the axial-vector current we have for the case that µ = i.For the temporal direction, the corresponding expression is One can then form a 2 × 2 matrix of kinematical coefficients multiplying G A (Q 2 ) and G P (Q 2 ), given by where the first row of the matrix is for µ = 0 and the second row for µ = i, while the first column gives the kinematic coefficients multiplying G A (Q 2 ) and the second column those multiplying G P (Q 2 ).For the case of the matrix element of the pseudoscalar current we have In the above expressions, E N is the energy of the nucleon and K is a kinematic factor given by Given the above momentum-dependence of the decomposition, we can average over all momentum components for a given Q 2 value, namely where p 2 runs over the possible q 2 k values, the symbol Σ stands for the average and we indicate above the symbol the indices of the sum, which are always spatial, and below the symbol the conditions to be satisfied such that values are included in the sum.We note that, while G P (Q 2 ) and G 5 (Q 2 ) can be extracted directly from Π P and Π 5 , respectively, is the only form factor accessible at zero momentum transfer, while all others need to be extrapolated to Q 2 = 0. Our strategy for extracting the three form factors is to perform a combined fit of the Πs at fixed Q 2 and express the ground state matrix elements in Eq. ( 41) in terms of the above linear combinations of form factors.

C. Renormalization
Matrix elements computed in lattice QCD need to be renormalized in order to relate to physical observables.In the twisted mass fermion formulation, we need the renormalization functions Z S for the renormalization of the pseudoscalar form factor G 5 (Q 2 ), Z P for the renormalization of the bare quark mass and Z A for the renormalization of the axial-vector current.We note that we do not use Z S since G 5 (Q 2 ) is evaluated in the scale-independent and ultra-violet finite combination m q G 5 (Q 2 ).In Figs. 7 and 11 where G 5 (Q 2 ) is shown without m q , it is only done for visualization purposes.In those cases, we use Z S computed as Z P /(Z P /Z S ) with Z P computed in the RI ′ scheme.This is because a direct evaluation of Z S in RI ′ is more difficult than for Z P due to increased hadronic contamination effects observed in the case of Z S .
We use methods based on Ward identities or on the universality of renormalized hadronic matrix elements, which are often referred to as hadronic methods, in order to compute ultra-violet finite renormalization factors, such as Z A and Z P /Z S .Hadronic methods are fully nonperturbative and require no gauge fixing, unlike the RI ′ scheme.For more details, we refer to Appendix B of Ref. [54], where this approach is used to extract the renormalization constants for the ensembles employed here.This approach is preferred to the usual RI ′ scheme because it provides much more accurate results on Z A and Z P /Z S .The RI ′ scheme is employed for the determination of Z P , as discussed in Ref. [55].For completeness, the values of the renormalization constants used in this work are collected in Table IV.
In what follows we will denote by G A (Q 2 ) and G P (Q 2 ) the renormalized form factors obtained by multiplying the lattice three-point functions of the axial-vector cur-rent by Z A .For G 5 (Q 2 ) we consider the combination m q G 5 (Q 2 ) that renormalizes with µZ S /Z P , involving only the ratio Z S /Z P that is determined accurately from hadronic matrix elements.The light bare quark mass µ takes values aµ = 0.00072, 0.00060, and 0.00054 for the cB211.72.64, cC211.60.80, and cD211.54.96 ensembles, respectively.

V. EXTRACTION OF FORM FACTORS
As described in Sec.IV B, bare form factors at each value of Q 2 are extracted from combined fits to the values of the two-and three-point functions, after we construct the averages given in Eqs. ( 48)- (51).Two-point functions are available for all source-sink separations, t s , while three-point functions are measured at selected values of t s , listed in Table II, and available for all t ins ∈ [0, t s ].Since the optimal fit range in t s and t ins may vary for each case, as well as the number of states needed to describe the correlation functions, we explore a wide parameter space in the fitting ranges and number of excited states included.Results are then combined using model averaging as described below.Specifically, at each value of Q 2 , we use the following fitting approach: N st : We perform either two-or three-state fits of all quantities, cutting the sum in Eqs. ( 34) and ( 35), to a maximum of i max = N st − 1 with N st ∈ {2, 3}.
t 2pt, min : We vary t 2pt, min , the lower bound in the fit of the two-point functions.The upper bound is taken to be the source-sink separation where the correlator becomes compatible with zero within 5σ.This upper maximum value varies from approximately 2.5 fm at t 3pt, min : We vary t 3pt, min , the smallest value of t s used for fitting the three-point functions.We fit to all t s ≥ t 3pt, min available.
t ins, 0 and t ins, S : We vary the number of insertion time slices from the source and the sink kept in the fit, using t ins ∈ [t ins, 0 , t s − t ins, S ].We only allow for t ins, 0 ≥ t ins, S since the energy gap at the source, where we have momentum, is expected to be smaller than the energy gap at the sink, where there is no momentum.At Q 2 = 0 we fix t ins, 0 = t ins, S .
N O : We vary the number of exponential terms when we perform three-state fits to the three-point functions, since certain overlaps may be sufficiently suppressed.The suppression rate is ordered according to the energy gaps of the first and second excited state energies.Beyond the ground state A 0,0 µ , the suppression increases for the terms containing the overlaps µ .We use either the first four, six, or all parameters, namely, we take N O ∈ {4, 6, 9}.N O = 4 corresponds to a full two-state fit.
In summary, N st and N O affect the number of parameters in the fit, while t 2pt, min , t 3pt, min t ins, 0 , and t ins, S the number of data used in the fit.We fit together the data for Q 2 = 0 and for the lowest non-zero value of Q 2 obtained when the momentum transfer in one spatial direction is 2π/L.After performing the model averaging for the zero and for the lowest non-zero value of Q 2 , we extract m N and use it as a prior to fit independently each larger Q 2 value.We note that the summation method is not used since taking into account only the ground state fails to describe the results in the case of G P (Q 2 ) and G 5 (Q 2 ), where excited state effects are large.We find that for the summation method to converge one would require source-sink time separations larger than 2 fm which are not available.

A. Model average
Results obtained using the different fit approaches are averaged using the Akaike Information Criterion (AIC) and we refer to Ref. [65,66] for a detailed introduction to the method.In the following, we summarize the practical aspects of our implementation.To each fit i, we assign a weight w i , defined as where N dof = N data − N params is the number of degrees of freedom, given as the difference between the number of data, N data , and the number of parameters, N params , used in the corresponding fit.We use correlated fits and, therefore, the χ 2 is defined as where, for each fit i, C i is the covariance matrix between the selected data ⃗ y i and ⃗ r i is the residual computed using the selected fit approach f i evaluated on the selected data range x i .From the weights in Eq. ( 52), we define the probability The model-averaged value of an observable O is given as where Ōi and σ i are, respectively, the central value and the error of the observable O measured using the parameters of the i th fit.The horizontal bands spanning both the left and right panels are the results of the model average among all fits using two states (red points and red band) and three states (blue points and blue band).The most probable fit is depicted with open symbols.In the left panel, we show for all ensembles, the result of the most probable fit using two states (red curve) and three states (blue curve) over the range used in the fit.Panels from top to bottom are for the cB211.72.64, cC211.60.80, and cD211.54.96 ensembles, respectively.

B. Selection of data and fits
We first illustrate our fitting procedure by considering the zero-momentum nucleon two-point function.In Fig. 1, we show the nucleon effective mass for each ensemble.We observe an impressive agreement among the data using the three ensembles showing very mild cut-off effects and compatible excited state contamination.This confirms that maintaining the radius constant of the Gaussian smearing, as shown in Table III, is a good strategy.
In Fig. 2, we show the nucleon effective mass separately for each of the three ensembles, as well as the values of the nucleon mass obtained via fits to two-and three-point functions keeping only those fits with model probability ≥ 1%.Note that since we perform a combined fit of the nucleon two-and three-point functions, for Q 2 = 0 and the lowest non-zero value of Q 2 as described at the beginning of this section, the values depicted in the figure are not obtained by fitting only the nucleon effective mass data shown in these figures.We observe a good distribution of the probabilities of the fits with the most probable fit having a probability between 10% to 50%.
Similarly, in Figs. 3 and 4, we show, respectively, results for Q 2 = 0 and the lowest momentum transfer for the axial form factor.The data shown in these figures are those obtained from the ratio Π AP (Q 2 , 0) defined in Eq. ( 49).We note that the results for the lowest non-zero value of the momentum transfer also have information from Π 0 (Q 2 ) and Π AP (Q 2 , (2π/L) 2 ) given in Eqs. ( 48) and (49).Results on the latter two are shown in Figs. 5  and 6, respectively, from which G P (Q 2 ) is extracted for the lowest non-zero value of Q 2 .In Fig. 7 we show the corresponding results for G 5 (Q 2 ) for the lowest non-zero We note that in most cases the model average gives a central value and error that is compatible with the corresponding ones of the most probable fit.In the few cases where there is a discrepancy, as e.g. in the three-state fit results for G A,1 in Fig. 4, the outcome of the model average has a Gaussian distribution and thus the model average procedure, given in Eq. ( 55), properly accounts for the systematics arising from the model choices.We demonstrate this by depicting in Fig. 8 the cumulative distribution of the Gaussians associated with each fit and the resulting Gaussian distribution outcome of the model average.We observe an excellent agreement.

C. Determination of axial charge and radius
Before presenting the analysis of the Q 2 -dependence of form factors, we perform fits to the zero and the lowest non-zero momentum transfers.For Q 2 = 0 only G A (Q 2 ) can be extracted, yielding the isovector axial charge g A .Computing the slope using the values of G A (Q 2 ) at these two values yields the radius ⟨r 2 A ⟩, namely where Q 2 1 is the lowest non-zero momentum transfer squared.Results obtained using two-or three-state fits are analyzed separately.The results on g A and r 2 A extracted after model averaging for each ensemble are collected in Table V and depicted in Fig. 9.We also include  42) that yields gA versus tins − ts/2 (left column) and versus ts for tins = ts/2 (middle column).In the header of the figure, we give the symbols used to denote the various values of ts/a.In the right column, we show the value of the nucleon isovector axial charge obtained via fits, as in Eq. ( 41), versus the fit probability, using the notation of Fig. 2. In the left and middle panels, the curves correspond to the fit results, which have the largest probability among all two-(blue) or three-(red) state fits.the results obtained for m N .We perform a linear extrapolation in a 2 to the continuum limit.We observe a very good agreement at the continuum limit between the -0.FIG. 6. ΠAP (Q 2 ) for the lowest non-zero value of Q 2 using the notation of Fig. 3.
results from fits using two-and three-states for all three quantities.On the other hand, there are slight deviations at finite lattice spacings between two-and three-state fits.
For this reason, we analyze separately all quantities using two-and three-state fits.We take as our mean value the one extracted using the two-state fit and we give as a systematic error the difference between the central values A ⟩, and nucleon mass mN for each ensemble and extrapolated to the continuum limit using a linear function in a 2 .These results are referred to as obtained via the "direct approach" in the text.Results are given separately for values extracted from a two-and a three-state fit analysis.
using two-and three-state fits.We find

D. Energy spectrum and dispersion relation
As customarily done in similar studies [12,18,19,23], we analyze the first excited state at the source and sink, E 1,⃗ p and E 1, ⃗ 0 , respectively.These are obtained using two-state fits to the three-point functions.Our results are depicted in Fig. 10 for the three ensembles, where we also depict the dispersion relation for the nucleon energy E N .We observe the following:  55).For each panel, all distributions are normalized with the maximum value of the cumulative distribution and we have doubled the height of the colored curves for visualization purposes.
• Since the dispersion relation is included using priors when fitting each Q 2 value larger than the lowest non-zero it is not surprising that we see excellent agreement between the extracted energy and the dispersion relation; • The first excited state at zero momentum transfer is compatible with the Roper; • E 1, ⃗ 0 for low values of the momentum is compatible with the lowest energy of the πN system in the rest frame, namely π and N moving with momentum p back to back; The dependence on the momentum is due to the strong enhancement of the πN excited state due to the pion pole [18]; As the momentum grows the lowest energy of πN becomes larger than the mass of the Roper and then E 1, ⃗ 0 becomes approximately constant somewhat above the mass of the Roper which is expected since in a two-state fit the first excited energy is contaminated by higher states; • E 1,⃗ p is compatible with N (0)+π(⃗ p) for all non-zero values of ⃗ p 2 ≤ 0.6 GeV 2 .After that, the energy of the Roper denoted by R(⃗ p) becomes smaller and the results tend to be in between the energy of the Roper and the energy of the N (0) + π(⃗ p) system.These are, indeed, the two lowest one-particle and two-particle excited state energies; • Excited state contamination is similar for all three ensembles and this is in line with the observation of mild cut-off effects for the nucleon mass.

E. Comparison of results extracted with two-and three-state fits
The renormalized form factors obtained using two-and three-state fits are compared in Fig. 11, where we also depict the difference ∆ between the results extracted using two-and three-state fit normalized such that errors are unity, namely with σ δ the jackknife error on the difference δ.We observe a very good agreement between the results extracted using two-and three-state fits with all differences ∆ lying within three standard deviations.We note that we only perform three-state fits to extract the form factors up to Q 2 ≃ 0.5 GeV 2 despite the fact that the stability of the three-state fits improves because time separations become more dense as the lattice spacing decreases. .These are determined by performing a two-state fit to the two-and three-point functions.To guide the eye, we depict the expected dispersion relations: the energy of a boosted nucleon EN (top), a non-interacting πN state with either total momentum zero, N (−⃗ p) + π(⃗ p) (middle), or with momentum ⃗ p, N ( ⃗ 0) + π(⃗ p) (bottom); and the energy or the first excited state of the nucleon, the Roper with mR = 1.44 GeV, either at rest (R(0), middle) or moving with momentum ⃗ p (R(⃗ p), bottom).
Since we observe consistency between the results from two-and three-state fits and the two-state fits do not suffer from instabilities at large values of Q 2 , we opt to take the results from two-state fits up to Q 2 ≃ 1 GeV 2 .The three-state fits are only considered up to Q 2 ≃ 0.5 GeV 2 , which is the largest Q 2 where the three-state fits for the cB211.72.64 are stable.

VI. Q 2 -DEPENDENCE AND CONTINUUM LIMIT
We will first discuss the parameterization of the Q 2 dependence of form factors that are free from the pion pole such as the axial form factor. Typically two functional forms are employed, the dipole Ansatz and the model independent z-expansion [68,69].The dipole Ansatz, given by has two parameters, the charge g and the dipole mass m.
In this case, the radius defined in Eq. ( 8) is given by We parameterize cut-off effects of the charge and of the radius using a linear function in a 2 , namely and obtain for the dipole Ansatz the following combined (Q 2 , a 2 )-dependence which we will use to fit all form factors after factoring in any pion pole dependence for a given lattice spacing.
In the case of the z-expansion, the form factor is parameterized as, where with −t cut < t 0 < ∞ and t 0 an arbitrary number and t cut the particle production threshold.For t cut , we use the three-pion production threshold, namely t cut = (3m π ) 2 [69] with m π = 0.135 GeV.For t 0 we use a vanishing value such that the charge is given by a 0 and the radius is proportional to the ratio a 1 /a 0 , namely We introduce the dependence on the lattice spacing by writing where c k = a k /a 0 and The coefficients c k can be further constrained by requiring that the z-expansion converges smoothly to zero at infinite momentum, namely [70]  From left to right we show results for the cB211.72.64, cC211.60.80, and cD211.54.96 ensembles.From top to bottom, we show results for the axial, induced pseudoscalar, and pseudoscalar form factors.In the lower panel of each plot, we include the difference ∆ defined in Eq. ( 58) between the results extracted using two-and three-state fits.The difference ∆ is normalized such that errors are unity and the dashed line represents a three standard deviation difference.Three-state fits are stable for Q 2 < 0.465 GeV 2 for all three ensembles.For larger values the fits become unstable.We display these points by using lighter blue color and we thus do not use them in the extraction of form factors.
This suggests that priors centered around zero should be used to help enforce this condition at various orders with a width that falls like 1/k [25].Additionally, an examination of the explicit spectral functions and scattering data [69] motivates the bound of |c k | ≤ 5. We, therefore, use the following Gaussian priors where w ≤ 5 is a fitting parameter that we vary together with the order of the expansion k max ∈ [1,4].
In both the dipole and z-expansion fits that follow, we will refer to one-and two-step fits.In two-step fits we first fit the Q 2 dependence for each lattice spacing separately and then take the continuum limit of the parameters, while in the one-step fits the three ensembles are fitted together.The one-step approach provides for a global χ 2 .

VII. THE AXIAL FORM FACTOR
We first present the analysis of the axial form factor, which at Q 2 = 0 yields the axial charge already discussed in Sec.II.12.The axial form factor obtained on each of the three ensembles using two-state fits (blue circles, orange downwards pointing triangles, and green triangles).The continuum limit form factor (red curve and band) and value for the axial charge (red cross) are obtained via dipole fits within the onestep approach.Also shown are the form factor curves obtained at the three values of the lattice spacing (blue, orange, and green curves, respectively).

A. Dipole Ansatz
An example fit using the dipole Ansatz is shown in Fig. 12, where we depict the results of using Eq. ( 62) to perform a combined fit of the form factors for the three ensembles, i.e. following a one-step approach.The values for G A (Q 2 ) shown in this figure are obtained using twostate fits in the range 0 As already mentioned, an alternative to the two-step approach is to perform a simultaneous fit to the Q 2 dependence for all three ensembles.To demonstrate that taking a global fit in a one-step approach is equivalent to performing a two-step approach, we show in Fig. 13 the continuum limit of the axial charge and the radius extracted from the one-and two-step approaches.In Table VI, we give the corresponding values extracted when using the one-and two-step approaches, including their reduced χ 2 .As can be seen, the continuum values are in perfect agreement both in terms of the central value and in terms of the error.The reduced χ 2 for the twostep procedure only refers to the linear extrapolation to the continuum limit.The one-step approach provides for a single value of χ 2 that reflects the quality of the fit to the combined Q 2 -and a 2 -dependence and it is thus more practical to compute the relative weights between the various fits when carrying out our model averaging.Therefore, from now on we will proceed with the one-step approach.
Using the one-step approach, we perform dipole fits to results obtained using two-and three-state fits to the correlators.We vary the largest Q 2 included in the fits, and for the two-state fit results, which are more precise, we also repeat the fits omitting the result at Q 2 = 0.The reasoning is that at Q 2 = 0 only G A (0) survives, which can affect the determination of the energy extracted for Values for the nucleon isovector axial charge gA, radius rA, and reduced χ 2 obtained using a dipole Ansatz to fit each ensemble (first three rows).We also give the continuum limit values using the one-step (second to last row) and two-step (last row) approaches.
the first excited state, as already shown in Fig. 10.A comparison of the results obtained using these variations is shown in Fig. 14.We perform a model average of the results separately for the case of using two-and threestate fits on the correlators.We find (34) (3-state) Since the values are compatible, we opt to quote the model-averaged values obtained from data that were extracted using two-state fits to the correlators.We then give as a systematic error the difference between the central values of the model-averaged results obtained from data extracted using two-and three-state fits to the correlators.We find g A = 1.196 (24) Results for the axial charge and radius extracted using the dipole Ansatz, versus Q 2 max , i.e. the maximum value of Q 2 included in the fit.We show separately dipole fits to results obtained from two-(red squares) and three-(blue diamonds) state fits to the three-and two-point correlation functions.For the two-state fit case, we include a variation in which the values of the form factors at Q 2 = 0 are not included in the fit (open squares).

B.
z-expansion

First order z-expansion
We repeat the same procedure using a first-order zexpansion that has the same number of parameters as the dipole Ansatz, and where no priors are employed.In Fig. 15 we demonstrate the one-step approach as an example, with the same notation as Fig. 13; in Fig. 16 and Table VII we similarly demonstrate that our one-step approach is equivalent to the two-step approach; and in Fig. 17 we depict the results as a function of Q 2 max for the case of data extracted using two-and three-state fits to the correlators.After model averaging we find g A = 1.283 (31) (2-state) = 1.249 (37) (3-state) Again, we observe that the values are compatible whether we use results from the two-or three-state fits to the correlators and thus we quote the model-averaged value for the case of using two-state fits and take as systematic error the difference between the central values obtained when using two-and three-state fits to the correlators.We find 1.32 a = 0, 2-step 1.291(24) 0.431 (12) 0.29 TABLE VII.Values for the nucleon isovector axial charge gA, radius rA, and reduced χ 2 obtained using a first-order zexpansion fit on each ensemble and extrapolated to the continuum limit using our one-or two-step approach.

Convergence of the z-expansion
In order to check the stability of the z-expansion fits, we study the convergence of the z-expansion as a function of the order and the amplitude of the priors used in Eq. (69).Priors are used to ensure convergence of the z-expansion [70] being centered around zero with a width that decreases with 1/k.We observe convergence for k max ≥ 3 for all cases.We vary the amplitude of the prior using w ∈ [1,5].In Fig. 18  axial charge and radius as a function of the prior width and Q 2 max , the largest Q 2 used in the fit.We observe that the results obtained by changing the width of the priors are all consistent.The result of the model average is g A = 1.245 (28) (2-state) = 1.231 (34) (3-state) quoting again the value from the data extracted from two-state fits with a systematic error in the difference between the central values of the model-averaged results when using data from two-and three-state fit to the correlators g A = 1.245(28)( 14) Following a similar procedure to the one-step approach we perform the fits in two steps, namely we first perform the fits and model average for each ensemble and then take the continuum limit.This two-step approach yields results that are compatible with the one-step approach, as depicted in Fig. 19 with results provided in Table VIII.

C. Final results
Having presented the variations used to extract the axial charge and radius at the continuum limit, we continue here to discuss the consistency among them and how we choose our final values.To summarize the variations, we have used: i) in Sec.V C a direct determination using the matrix element at Q 2 = 0 and, for the radius the  (25) 0.58 a = 0, 1-step 1.245(28) 0.339 (48) 0.69 a = 0, 2-step 1.256(27) 0.333 (51) 0.05 TABLE VIII.Values for the nucleon isovector axial charge gA and radius rA, and reduced χ2 of the best fit obtained using a model average over third-order z-expansion fits on each ensemble and extrapolated to the continuum limit using the one-or two-step approaches.
matrix element at the lowest non-zero Q 2 value yield-ing the results of Eq. ( 57); ii) in Sec.VII A the dipole Ansatz to describe the Q 2 -dependence resulting in the values given in Eq. ( 71); iii) in Sec.VII B 1 the first-order z-expansion to describe the Q 2 -dependence resulting in the values given in Eq. ( 73); and iv) in Sec.VII B 2 the higher-order z-expansion with the resulting values given in Eq. ( 75).In Table IX, we collect these values and depict them in Fig. 20.In all cases, the central value and first error are obtained via model averaging of the two-state fit data, while the second error is a systematic obtained as the difference between the central values when model averaging the two-or three-state fit data.We note the following • We opt to take as our final values the results extracted from using the two-state fits because they are in general more stable, having fewer fit parameters and smaller errors and, thus, allowing a better study of the convergence to the ground state.
We reiterate that we observe a very good agreement between results extracted using data determined by fitting correlators to two-or three-states, as depicted in Fig. 11.
• All results agree with the value extracted using the z k -expansion with k = 3 within error bars.Results using the dipole Ansantz and the z 1 -expansion yield respectively smaller or larger values as compared to those using the z 3 -expansion.This observation is compatible with what has been found in another study [23].
• Furthermore, the results using the z 3 -expansion are completely consistent with those from the direct approach.Since the direct approach uses the matrix element at Q 2 = 0 for g A and for ⟨r 2 A ⟩ the slope using also the lowest non-zero value of Q 2 , it does not depend on any Ansatz used to fit the Q 2dependence of the form factor.The fact that the z 3 -expansion yields the same results shows that indeed it provides a model-independent approach to extract the same information on these two quantities.The error when using the z 3 -expansion is smaller compared to the errors when using the direct approach since the z-expansion makes use of more information.
Given the above observations, we quote as our final values the results from the z k -expansion that has shown convergence for k = 3 and is model-independent.Thus, we take as our final values for the axial charge and radius g A = 1.245 (28) (14) [31] ⟨r 2  A ⟩ = 0.339(48)(06) [48] fm 2 , (final value) (76) where in the square brackets we have combined quadratically the two errors.We also collect all our final values in the conclusions in Eq. ( 97).In Figs.21 and 22 we show the model-averaged results as a function of Q 2 when using the z 3 -expansion for either two-or three-state fits to the correlators, respectively.As can be seen, in Fig. 11 the data extracted using two-and three-state fits are compatible.However, small statistical fluctuations and the lack of data in the case of the threestate fit analysis for Q 2 > 0.5 GeV 2 can affect the fits of the Q 2 -dependence and thus the continuum limit, given that we only have three lattice spacings.We remind the reader that we perform simultaneously fits to the Q 2dependence for each ensemble and at the same time take the continuum limit.We compare the resulting G A (Q 2 ) in the continuum limit for these two cases in Fig. 23, where we give the continuum fits only.The fit parameters of the curves corresponding to the standard form of the z-expansion in Eq. ( 63) with (77) As can be seen, the resulting curve for the three-state fit case is in agreement with that for the two-state fit for Q 2 ≤ 0.5 GeV 2 .Also, the parameters of the two fits are in good agreement.Since however, the three-state fits become unstable for Q 2 > 0.5 GeV where we have used the correlation matrix of the parameters from two-state fit data.More information on the form factors at the continuum limit is provided in Appendix A.
We include the resulting G A (Q 2 ) when we assign this systematic error to the parameters of the continuum fit in Fig. 23.We consider the values of G A (Q 2 ) including the systematic uncertainty as our final results.Our final results for the form factor G A (Q 2 ) are given in Table XII.We will adopt the same strategy for the analysis of the other two form factors and for checking the PCAC and PPD relations.FIG.23.Results on GA(Q 2 ) at the continuum limit when fitting data extracted from the two-(red band) and three-(blue band) state fit analysis of the correlators.The darker blue curve indicates up to which Q 2 we had data for the three-state fit analysis.The yellow band is when we added systematic errors to the parameters that define the red curve as discussed in the text.The parameters of the fit are given in Eq. (78).
We perform a similar analysis to determine G P (Q 2 ) and G 5 (Q 2 ) to the one discussed in detail above for G A (Q 2 ).The additional complication in the case of π , which needs to be removed before proceeding to apply similar fit functions to the ones applied for G A (Q 2 ).Therefore, before proceeding with the Q 2 -dependence analysis of these form factors, we present a detailed study of pion pole dominance.

A. Pion pole dominance (PPD)
The pion pole dominance (PPD) hypothesis introduced in Sec.II can be tested by forming two ratios of form factors, one of which is arising from Eq. ( 20) and the second r PPD,2 derived in Eq. ( 26) assuming a non-zero Goldberger-Treiman discrepancy.Using the results for the form factors from the two-state fits to the correlators we find the ratios depicted in Fig. 24.We indeed observe for both ratios a linear dependence in Q 2 , as expected from Eq. ( 79) and Eq. ( 27), respectively.We also observe clear cut-off effects for the first ratio whereas for the second the results from the three ensembles are consistent among them.To capture the a-dependence we fit the ratios using the functional form where we include the leading order a dependence to both the intercept and the Q 2 -slope.We also perform fits  80) to take into account cut-off effects.The red band is the continuum extrapolation after performing the model average as described in the text.The dashed line shows the expected value of the ratios if PPD close to the pion pole is satisfied, namely a slope of 1/4m 2 N from Eq. ( 79) for the first ratio and unity for the second.
where we set b 2 and c 2 to zero to account for the fact that the second ratio r PPD,2 shows no detectable cut-off effects to the accuracy of our data.We then perform a model average over all fits where we both include and exclude cut-off effects as well as change the largest Q 2 value, Q 2 max , used in the fit.The resulting continuum fits are shown in Fig. 24.
The conclusions drawn from Fig. 24 are • The pion pole dominance hypothesis is satisfied for both ratios at the pole since we obtain the expected value at • For the first ratio, G A (Q 2 )/G P (Q 2 ), we find a slope and intercept in the continuum limit consistent with the PPD hypothesis.Namely, from the intercept at Q 2 = −m 2 π , we find a value of the pion pole mass at the continuum limit of m pole π = 0.141 (20) GeV ≈ 0.135 GeV (82) compatible with the physical pion mass; and from the Q 2 -slope, whose value is expected to be 1/4m Pion pole mass extracted from the ratio GA(Q 2 )/GP (Q 2 ) for each ensemble, compared to the simulated unitary pion mass and the Osterwalder-Seiler (OS) pion mass.
we find a nucleon mass of m N = 0.9401(39) GeV ≈ 0.938 GeV , when the fit is done with Q 2 max = 0.3 GeV 2 .Thus, we conclude that the ratio G A (Q 2 )/G P (Q 2 ) satisfies the PPD relation close to the pole.
• Examining the pion pole mass determined from the fits to a given ensemble at finite lattice spacing, we find a significantly different pion pole mass.It is well-known that the twisted mass fermion formulation has significant cut-off effects in the pion mass [71] but much milder ones in other quantities such as other hadron masses [72] or hadronic operator matrix elements [71].In Table X, we give the pion masses that we find per ensemble as well as the values of the unitary pion mass, denoted by m TM π and the Osterwalder-Seiler (OS) pion mass [73,74], denoted by m OS π .We would like to clarify that the difference between the charged pion mass m TM,+/− π and its neutral counterpart m TM,0 π is an O(a 2 ) artifact due to the breaking of isospin symmetry in the clover twisted-mass fermion lattice action formulation.We find that the mass difference between charged and neutral unitary pions m TM,−/+ π − m TM,0 π is of order 20-40 MeV at the lattice spacings employed here.The uncertainty in the determination of this mass difference arises due to the statistical error of the disconnected quark contribution entering in the computation of m TM,0 π .While, the larger difference between m TM,+/−,0 π and m OS π is also an O(a 2 ) cutoff effect of the OS mixed action, the coefficient is different.Technically, the difference between the neutral m TM,0 π and m OS π can be traced back to the presence in the two-point correlator of the neutral TM pion of quark-disconnected contributions of O(a 2 ) that are absent in the two-point correlator of the OS pion.
We observe that the pion pole mass that we find is very close to m OS π .This is because, in our evaluation of G P (Q 2 ), we use the flavor diagonal isovector current and neglect the noisy O(a 2 ) quark disconnected contributions, which in turn corresponds to computing the three-point correlators in the OS mixed action formulation.We, thus, obtain a neutral pion pole with mass given by the OS pion mass, m OS π .In this way, indeed, we can understand the large cut-off effects observed for form factors that have a pion pole behavior within the twisted mass formulation with OS-type valence quarks in contrast to other fermion discretization schemes that observe similar cut-off effects in both G A (Q 2 ) and G P (Q 2 ), as when e.g. using Clover Wilson fermions [23], although in their formulation the form factors have O(a) discretization errors.
• We demonstrate that cut-off effects are due to the pion pole in Fig. 25, where we show results for the ratio r PPD,1 defined in Eq. ( 25) removing the pole using either the unitary or the OS pion mass.Data obtained using the OS pion mass shows, indeed, very mild cut-off effects and yield a value of r PPD,1 close to unity as expected by PPD and as observed by other groups using clover fermions.
• The deviation from unity observed for the ratio, r PPD,2 , is connected to the Goldberger-Treiman discrepancy as given in Eq. ( 28) where we followed a similar analysis to that of Ref. [48].Using the slope of our fit at the continuum limit to extract ∆ GT from Eq. ( 28) and d18 from Eq. ( 23), we find We note that in extracting these values we use our final values for g A and ⟨r 2 A ⟩ given in Eq. (76).For both ∆ GT and the low energy constant d18 we find values that are compatible with chiral perturbation theory, which predicts for ∆ GT ∼ 2% and for −1.40(24)GeV −2 < d18 < −0.78(27) GeV −2 depending on the type of fit used in the determination [46].Our determination gives a more precise value and provides valuable input for chiral perturbation theory.
• The mild cut-off effects observed for the ratio r PPD,2 in Fig. 24 is understood by the fact that this ratio involves G P (Q 2 ) and G 5 (Q 2 ) both of which have the same pion pole mass dependence, thus canceling the cut-off effects.
B. Parametrizations for the fits at the continuum of GP (Q 2 ) and G5(Q 2 ) In the previous section, we made three important observations, namely: i) G P (Q 2 ) and G 5 (Q 2 ) have the same pion pole mass dependence at each lattice spacing; ii) the pion pole mass obtained using valence OS quarks in the mixed action twisted fermion mass formulation shows significant cut-off effects, much larger than the mass splitting between the unitary charged and neutral pion; and iii) in the continuum limit we obtained a pion pole consistent with the physical pion mass of m π = 0.135 GeV.Based on these observations, we use the following functional form , where for G res we used the z k -expansion and repeat the same analysis presented for The combination m q G 5 (Q 2 ) is scale-independent and renormalizes with Z S /Z P , which is accurately determined.Furthermore, scaling by 1/m 2 π takes into account the slight difference in the unitary pion mass for each ensemble (see Table I) and by m N makes the whole combination dimensionless.Since the pion pole mass is the same for both G P (Q 2 ) and G 5 (Q 2 ), we in addition, perform a combined fit of both form factors taking the parameter "b" of the pole to be a common fit parameter.Since, as demonstrated in the previous section, the PPD relation is satisfied at the continuum limit, we perform fits where we enforce the value of g πN N extracted from both form factors to take the same value at the continuum limit.This is implemented by using a z 3 -expansion with t 0 = −m 2 π and, therefore, since, according to Eq. ( 64), z(−m 2 π ) = 0 and a 0 is a fit parameter as given in Eq (63).The coupling constant g πN N is then the same for both form factors if namely by making a 0 a common fit parameter for both G P (Q 2 ) and G5 (Q 2 ).26.Induced pseudoscalar coupling, g * P , and pionnucleon coupling, gπNN from a z 3 -expansion fit as a function of the largest Q 2 used in the fit, Q 2 max , and the width of the priors.For each Q 2 max we depict five points having prior width of w = 1, 2, 3, 4, 5.The points are shifted to the right as w increases with an increasing symbol size.

C. Convergence of the z-expansion
As for the analysis performed for G A (Q 2 ), we study the convergence of the z k -expansion as a function of the order k, the width of the priors used, and the largest Q 2 employed in the fit.We first discuss results when we fit separately G P (Q 2 ) and G5 (Q 2 ) without enforcing to have the same pole or the same value of g πN N and we monitor convergence by looking at the values it takes fitting G P ( 2 ), and G5 (Q 2 ).We also monitor the value of g * P , which is extracted from G P (Q 2 ) using Eq. ( 9).In Fig. 26, we show results on these quantities when we use data determined from the two-and three-state fit analysis to the correlators.We observe convergence for k max ≥ 3, stability in the values we extract as a function of Q 2 max and the width of the priors for both data from the twoand three-case analysis fits.After model averaging we find All values determined using data from the two-and These results are in agreement with those where we did not enforce the value of g πN N and have smaller uncertainties thanks to the combined fit approach.Since we have demonstrated that the PPD relation is satisfied at the continuum limit, we opt to quote these as our final results for these quantities.We follow the same strategy as for G A (Q 2 ) and quote the model-averaged value determined from using the data from the two-state fit analysis of the correlators and take as a systematic error the difference between the model-averaged central values of the data from the two-and three-state fits.We find the following values g * P = 8.99(39) (49)[63] g πN N = 13.25 (67) (69) [96] (final value) , where in the square brackets we have summed in quadrature the two errors in parentheses.In Figs.27 and 28 we depict results on G P (Q 2 ) obtained after taking the model average using the data from the two-or three-state fit analysis, respectively.In Figs. 29 and 30 we show the corresponding results for G 5 (Q 2 ).
In Table XI, we quote the values of the pion pole masses per ensemble as extracted from the individual or combined fit of G P and G5 from two-or three-state fit data.We observe an overall good agreement within errors, and the values confirm the agreement already discussed in Sec.VIII A with the OS pion mass reported in Table X.
One can drastically reduce cut-off effects by considering for G P (Q 2 ) and G 5 (Q 2 ) the following modified expression ) using the z 3 -expansion to fit the Q 2 -dependence of the data determined from the three-state fit analysis of the correlators up to Q 2 = 0.47 GeV 2 .The notation is the same as that for Fig. 27.
In Fig. 31, we show the improved expressions for form factors per ensemble.We observe that upon using the improved expression defined in Eq. ( 92), the results per ensemble are compatible with each other and with those obtained in the continuum limit by extrapolating G wpole (Q 2 ).These findings further confirm the interpretation that the sizable cut-off artifacts in G P (Q 2 ) and G 5 (Q 2 ) stem from the cutoff effects in using the OS pion mass for the pole.Since at finite a we neglect disconnected O(a 2 ) terms in our form factor computation this is indeed the expected behavior and fully justifies our fit Ansatz in Eq. ( 85) for the continuum extrapolation of the data when using G wpole .
D. Continuum results for GP (Q 2 ) and G5(Q 2 ) We follow the same procedure as the one for G A in Sec.VII C to arrive at the Q ) as defined in Eq. ( 86) for each ensemble (blue band for the cB211.72.64, orange band for the cC211.60.80 and green band for the cD211.54.96 ensemble) and at the continuum limit (red band) using the z 3 -expansion to fit the Q 2 -dependence of the data determined from the twostate fit analysis up to Q 2 = 1 GeV 2 .The inner panel shows a zoom-in of the region marked by the red square.) as defined in Eq. ( 86) using the z 3 -expansion to fit the Q 2 -dependence of the data determined from the three-state fit analysis of the correlators up to Q 2 = 0.47 GeV 2 .The notation is the same as that of Fig. 29 ) (open symbols) compared to those when using G improved (Q 2 ) (filled symbols) of Eq. ( 92) by correcting for the pole OS pion mass.The continuum limit form factors (red band) are those determined in Fig. 27 and Fig. 29  ) GP (Q 2 ) at the continuum limit when fitting data extracted from the two-(red band) and three-(blue band) state fit analysis of the correlators.The darker blue curve indicates up to which Q 2 we had data for the three-state fit analysis.The yellow band is when we added systematic errors to the parameters that define the red curve as discussed in the text.The parameters of the fit are given in Eq. ( 94).
G P (Q 2 ) and G 5 (Q 2 ).In particular, in Fig. 32, we show the corresponding results for G P (Q 2 ) as those shown for G A (Q 2 ) in Fig. 23 showing the comparison of results obtained when using data from the two-and three-state fit analysis after removing the pole, namely we show results for (Q 2 + m 2 π )G P (Q 2 ).As in the case of G A (Q 2 ), the data from the three-state analysis of the correlators are in agreement with those from the two-state fit analysis.However, after the continuum extrapolated results using the data from the three-state fit analysis yield systematically smaller values for G P (Q 2 ) for higher Q 2 values.As we already pointed out, the three-state fit analysis becomes unstable for Q 2 > 0.5 GeV affecting the fits to the Q-dependence.Since cut-off effects are larger for G P (Q 2 ) the slope in linear a 2 extrapolation is larger and thus more affected by small fluctuations in the data given that we also only have three lattice spacings.This explains why the continuum results in the two cases differ by up to a standard deviation at large Q 2 while the lattice data for the three ensembles are compatible.Having higher statistics will enable us to extract more reliable results using the three-state fit procedure and having more lattice spacings will better control the continuum extrapolation, something that we plan to do in the future when more computational resources are available.
(93) As can be seen, the parameters are in agreement albeit some carry large statistical errors and thus, we follow the same strategy as for G A (Q 2 ) for determining the best The values of G P (Q 2 ) that result from this parametrization are given in Table XIII of the Appendix.
Repeating the same analysis for G5 (Q 2 ), we find the results shown in Fig. 33, after removing the pole, namely we show results for (Q 2 + m 2 π ) G5 (Q 2 ).The behavior of the continuum limit results is the same as that observed for G P (Q 2 ) since both have similar cut-off effects due to the pion pole dominance.The fit parameters of the twoand three-state fit data curves are given by The values of G5 (Q 2 ) that result from this parametrization are given in Table XIV of the Appendix A, where we also provide more information on the form factors at the continuum limit.

E. Continuum limit of the PCAC and PPD relations
Having determined the three form factors G A (Q 2 ), G P (Q 2 ), and G 5 (Q 2 ), we can check the PCAC relation at the continuum limit.We use the values of the fit parameters of the z 3 -expansion to the Q 2 -dependence after taking the model average for each ensemble.We use the form factors extracted from the two-state fits to correlators.We also repeat using the three-state fits correlators.In both cases, we also take the continuum limit of the parameters determined at each lattice spacing, as previously discussed.In Fig. 34, we depict the resulting r PCAC as a function of Q 2 using data from the twoand three-state fit analysis, upper and lower panels, respectively.As can be seen, in both cases the PCAC relation is recovered in the continuum limit.In addition, we obtain the PCAC ratio in the continuum limit by using the final parameterizations of the form factors that take into account the systematic uncertainty due to how we treat excited states, i.e. difference of central values when we use two-or three-state fits, namely the results shown by the yellow band of Figs.23

, 27 and 29 for G
, and (m 2 π + Q 2 ) G5 (Q 2 ), respectively.As expected, the PCAC relation is recovered but the systematic error due to the treatment of excited states increases the error band.For comparison, we plot in Fig. 35 in the same format, the results for the ratio r PPD,1 .It is no surprise that it also fulfills the PPD dominance in the continuum limit, as already discussed in relation to Fig. 24.As in the case of r PCAC we show both the continuum limit curve extracted using the data from the two-state fit analysis and the one when we include the systematic uncertainty difference between the central values of the fit parameters determined by using data from to the two-and three-state fit analysis.

IX. COMPARISON WITH OTHER RESULTS
Before comparing with other lattice QCD studies, we compare in Fig. 36 our older results [19] where only the cB211.72.64 ensemble was used to the ones obtained in this work.As can be seen, the central values are in agreement showing that cut-off effects are mild for these quantities.The error on g A increases after taking the continuum limit, while the error on the axial radius is approximately the same.The fact that the errors on g * P and g πN N are much smaller is a combination of two things: i) taking the continuum limit and ii) in our previous work we used the PCAC and PPD relation and lattice QCD data on G A (Q 2 ) which is more precisely determined.The reason was that with one lattice spacing, we could not account for the large cut-off effects on G P (Q 2 ) and G 5 (Q 2 ) leading to a violation of the PCAC relation.In this work, g * P and g πN N are determined directly from our data on G P (Q 2 ) and G 5 (Q 2 ), although, as shown in this work, in the continuum limit the PCAC relation ) for the cB211.72.64, cC211.60.80 and cD211.54.96 ensembles, respectively.The fits are done using the z 3 -expansion to fit the Q 2 -dependence of the data determined from the two-state fit analysis up to Q 2 = 1 GeV 2 .The red curve and band show the results at the continuum limit.The yellow band shows the errors on the curve when systematic effects are taken into account, i.e. using the final parameterization of the form factors given in Eqs. ( 78), (94), and (96) for GA( , respectively.Bottom: Same as top panel, but using three-state fit data and the respective z 3 -expansion to fit the Q 2 -dependence up to Q 2 ∼ 0.5 GeV  25).The notation is the same as that in the top panel of Fig. 34.
holds and could be used to determine them.We note  [23], the green triangles by RQCD [18,75], the yellow squares by NME [48], the gray diamonds by the Mainz group [22] and the magenta square by CalLat [14].The cyan circle shows the FLAG21 average of lattice results published at the time of the report [76].
that the trend that we observe of errors becoming larger in a number of quantities highlights the importance of having results using ensembles with smaller lattice spacings.However, as it is well known, simulations for lattice spacing a < 0.05 fm become difficult due to the increase in the autocorrelation time.There are ongoing efforts to address the critical slowing down of Hybrid Monte Carlo simulations [77][78][79].
The nucleon axial charge and radius as well as the coupling constants g * P and g πN N are compared to other recent lattice QCD results in Fig. 37.We selected studies that provide results at the continuum limit and at the physical pion mass, either computed directly like ours or via a combined chiral and continuum extrapolation.There is a nice agreement among all lattice QCD results on these quantities which are defined either at Q 2 = 0 or at the limit Q 2 → 0.
In Fig. 38 we compare our final parameterization of G A (Q 2 ) given in Eq. ( 78) with fits to experimental data on G A (Q 2 ) and with fits to data computed by other lattice QCD groups.When compared to experimental data, our results fall off slower than the fits to experimental data.While our results are within two standard deviations as compared to the recent results from the Minerνa experiment [1], they show more tension with the fit to the deuterium bubble-chamber data [25].In addition, our results are in good agreement with the results by the 0.00 0.25 0.50 0.75 1.00 1. 25  ) determined within this work using the parameters a2-state (red solid line and band) and when including the systematic uncertainty as the difference in the central values of a2-state and a3-state (yellow band).We compare to the fit to the deuterium bubble-chamber data [25] shown by the green dashed line with error band and with the fit to the recent MINERνA antineutrino-hydrogen data [1] shown by the blue dot-dashed line with error band.Bottom: GA(Q 2 ) determined within this work compared to two recent lattice QCD calculations: i) by the Mainz group [22] shown with the gray dashed line with its error band, and ii) by PNDME [23] shown with the blue dashed line with its error band.
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 This work PPD FIG.39. Results on GP (Q 2 ) determined directly from our lattice data using the parameters a2-state (red solid line and band) and when including the systematic uncertainty as the difference in the central values of a2-state and a3-state (yellow band) are compared to those obtained by using PPD given in Eq. ( 20) and to our data on GA(Q 2 ) (dark blue and light blue when including the systematic error).
Mainz group [22] and close to the results by PNDME [23] and NME [48,80].We comment below on some aspects of the lattice QCD calculations: • The results of this work are the only ones that are extrapolated to the continuum limit using only ensembles simulated directly with physical pion mass.
• The rest of the collaborations combined results extracted using ensembles simulated with larger than physical pion masses to extrapolate to the physical pion mass and to the continuum limit.Specifically, NME [48] uses no physical point ensembles, the Mainz group [22] uses one with their physical point results having large errors and RQCD [18] and PNDME [23] use two physical pion mass ensembles.
• In the case of the RQCD [18] and PNDME [23], results using the physical pion mass ensembles have larger statistical errors as compared to those of ensembles with heavier than physical mass.This means that results at the physical point weigh less in the extrapolation.Additionally, in both studies, the form factors are computed using the physical pion mass ensembles only for Q 2 ≲ 0.3 GeV 2 and information at high Q 2 is provided by a subset of the ensembles.We instead compute the form factors up to Q 2 = 1 GeV 2 .
• The PNDME collaboration [23] uses a N f = 2+1+1 mixed action of clover fermions on a staggered sea.They have employed thirteen ensembles simulated at four values of the lattice spacing, three values of the pion mass (135, 220, and 310 MeV), and volumes with 3.7 ≤ m π L ≤ 5.5.Their axial-vector current is unimproved which means they have O(a) cut-off effects.Nevertheless, they observe mild dependence on the lattice spacing.They also do not observe any significant lattice volume dependence, while they do see a stronger dependence on the pion mass.In this work, we use ensembles with approximately the same volume, namely with 3.6 ≤ m π L ≤ 3.9.Given the volume study by PNDME [23], we expect finite size effects on our results to be small.PNDME also carried out an elaborate study of excited states highlighting the effects of πN states and concluded that an approach compatible with the one employed in this work is the most suitable.Namely, they propose performing a combined fit of all matrix elements at the same Q 2 using common fit parameters for the excited states.However, they do not include a systematic error due to excited states and this explains why their results have a smaller error band.
• The Mainz group [22] has used fourteen CLS N f = 2 + 1 ensembles simulated with clover-improved Wilson fermions at four values of the lattice spacing, pion masses in the range 130 MeV ≤ m π ≤ 350 MeV, and volumes with 3.9 ≤ m π L ≤ 5.9.Their current is O(a)-improved and they observe a mild dependence on the lattice spacing and the lattice volume, while a stronger dependence on the pion mass, which requires the inclusion of higher order corrections that are not considered by PNDME.
The Mainz group also includes a systematic error in a similar way to what we do.
• The RQCD collaboration [18] uses the same CLS ensembles as the Mainz group, but thirty-seven of them, having five values of the lattice spacing, pion masses in the range 130 MeV ≤ m π ≤ 420 MeV and volumes with 3.5 ≤ m π L ≤ 6.4.Their physical point limit also involves a limit to the physical strange quark mass that is not required in the set of ensembles used by the Mainz group.They perform a thorough study of excited states including the effect of πN states.They also observe a strong dependence on the pion mass.They have provided results using dipole fits (labeled !2P ) or using zexpansion (labeled !z 4+3 ), without selecting one of the two as the final value.For this reason, we report two sets of points in Fig. 37 for RQCD19.In their recent work [75], they have provided a new value for the axial charge extracted from the analysis of matrix elements at zero momentum transfer and from an analysis of ten additional ensembles having four at the physical pion mass.
• NME [48] has used seven N f = 2 + 1 Wilsonclover fermions ensembles simulated at five values of the lattice spacing, with pion masses in the range 166 MeV ≤ m π ≤ 285 MeV and volumes with 3.9 ≤ m π L ≤ 6.2.Their axial-vector current is not O(a)-improved and they observe strong lattice spacing effects and pion mass dependence.The analysis of the excited states is compatible with the one carried out by PNDME and they include an elaborated study of excited states using priors around the πN excited states.For this case, we only show the values of the coupling constants and axial radius in Fig. 37.
• The PACS collaboration computed G A (Q 2 ), G P (Q 2 ) and G 5 (Q 2 ) using N f = 2 + 1 cloverimproved fermions and a large spatial volume of length L = 8.1 fm [17] and pion mass of 146 MeV and lattice spacing of a = 0.085 fm.However, they only have time separations up to t s ∼ 1.3 fm and only perform plateau fits to individual form factors.They also have results for these form factors for lower Q 2 values up to ∼ 0.25 GeV 2 .Since their results are given only at one lattice spacing they are not included in the comparisons.
While G A (Q 2 ) is determined by a number of lattice QCD collaborations, there are scarce results on G P (Q 2 ) and, to our knowledge, this work is the first to compute G 5 (Q 2 ) at the continuum limit.Experimental studies also probe G A (Q 2 ) and one could use PCAC and PPD to estimate G P (Q62).In Fig. 39, we show our results from a direct evaluation of G P (Q 2 ) in comparison with the ones obtain using our data on G A (Q 2 ) and the Eq. ( 20) to extract G P (Q 2 ).As can be seen, the results are in perfect agreement with the uncertainties.Therefore, one would be justified to use Eq. ( 20) and the experimental data on G A (Q 2 ) to estimate G P (Q 2 ).

X. CONCLUSIONS
In this work, we present results on the axial, induced pseudoscalar, and pseudoscalar form factors in the continuum limit.This study is performed using three N f = 2 + 1 + 1 twisted-mass ensembles with all quark masses tuned to their physical value and simulated at three values of the lattice spacing.Our analysis is also done up to Q 2 = 1 GeV 2 as compared to some other lattice QCD studies where, for physical point ensembles, only smaller values of Q 2 were accessible.Our final values for the nucleon axial charge and radius as well as the coupling constants g where the central values and the first error in the parenthesis are obtained from an analysis of data extracted from the two-state fits to the correlators, the second error is the systematic error due to the excited states computed as the difference between the central values from using the data extracted from the two-and three-state fit analysis of the correlators and the third error in the square brackets is the total error obtained by summing in quadrature the first two.In Appendix A, we provide the final parameterization and values of the form factors at the continuum limit with and without systematic uncertainties due to the excited states.
From our analysis of the ratio r PPD,2 defined in Eq. ( 27) we also determine the values of the Goldberger-Treiman discrepancy and the low-energy constant d18 ∆ GT = 2.13(38)% d18 = −0.73(13)GeV −2 . (98) Our results on G A (Q 2 ) are in good agreement with other recent lattice QCD studies.Having taken the continuum limit using only ensembles at the physical point mass we avoid a chiral extrapolation that in the nucleon sector can lead to an uncontrolled systematic error.An advantage of this setup is that it allows us to directly access cut-off effects.We find that for G A (Q 2 ), cut-off effects for the range of lattice spacings used are mild, ranging from not detectable within our errors at low Q 2 to slightly positive at high Q 2 .On the other hand, the induced pseudoscalar and the pseudoscalar form factors exhibit similar large cut-off effects that can be traced back to the known O(a 2 ) artifacts on the pion mass pole.Such cut-off effects are expected in the twisted mass fermion formulation used in this work for the computation of these form factors and as such they can be conveniently parameterized in our continuum extrapolation fits.Alternatively, they can be also substantially reduced by considering modified expression for the pole of the form factors.As shown in this work, the important conclusion is that in the continuum limit, all cutoff effects are safely eliminated as expected.In particular, both the pion pole dominance close to Q 2 = −m 2 π , with m π = 135 MeV, and the fundamental PCAC relation that follows from QCD chiral Ward identities, are fully recovered.Regarding finite volume effects, we would like to point out that the PNDME collaboration [23] and the CLS-Mainz group [22] investigated finite volume effects for these quantities and found that for the range of m π L > 3.5 of our three ensembles, no detectable finite volume effects were observed within the accuracy of their lattice data, that is similar to ours.

FIG. 1 .
FIG.1.Nucleon effective mass using the three physical point ensembles.The dashed line is the value of the nucleon mass mN = 0.938 GeV.

kmax k=0 c k d n z k dz n z=1 = 0 FIG. 11 .
FIG. 11. Results for the three form factors obtained using a two-(red squares) or a three-state (blue crosses) fit analysis.From left to right we show results for the cB211.72.64, cC211.60.80, and cD211.54.96 ensembles.From top to bottom, we show results for the axial, induced pseudoscalar, and pseudoscalar form factors.In the lower panel of each plot, we include the difference ∆ defined in Eq. (58) between the results extracted using two-and three-state fits.The difference ∆ is normalized such that errors are unity and the dashed line represents a three standard deviation difference.Three-state fits are stable for Q 2 < 0.465 GeV 2 for all three ensembles.For larger values the fits become unstable.We display these points by using lighter blue color and we thus do not use them in the extraction of form factors.

FIG. 22 .
FIG.22.Continuum limit of GA(Q 2 ) using the z 3 -expansion and data from the three-state fit analysis of the correlators up to Q 2 = 0.47 GeV 2 for the three ensembles with the symbols as indicated in the header of the figure.
FIG.26.Induced pseudoscalar coupling, g * P , and pionnucleon coupling, gπNN from a z 3 -expansion fit as a function of the largest Q 2 used in the fit, Q 2 max , and the width of the priors.For each Q 2 max we depict five points having prior width of w = 1, 2, 3, 4, 5.The points are shifted to the right as w increases with an increasing symbol size.

2 )FIG. 33 .
FIG.33.Results on (Q 2 + m 2 π ) G5(Q 2 ) at the continuum limit when fitting data extracted from the two-(red band) and three-(blue band) state fit analysis of the correlators.The darker blue curve indicates up to which Q 2 we had data for the three-state fit analysis.The yellow band is when we added systematic errors to the parameters that define the red curve as discussed in the text.The parameters of the fit are given in Eq. (96).

TABLE III .
The number of Gaussian smearing iterations nG and the Gaussian smearing coefficient αG used for each ensemble.We also provide the number of APE-smearing iterations nAPE and parameter αAPE applied to the links that enter the Gaussian smearing hopping matrix.The resulting source r.m.s.obtained is given in the last column, where the error is due to the uncertainty in the lattice spacing.

TABLE IV .
[55]es of the scheme-independent renormalization constants ZA and ZP /ZS taken from Ref.[54]and of the scheme-dependent ZP given in MS at µ ref = 2 GeV computed in Ref.[55].
(top)and first excited state E 1, ⃗ 0 (middle) at sink and E 1,⃗ p (bottom) at source as a function of ⃗ p 2

TABLE VI .
13G.13.The axial charge (top) and radius (bottom) obtained via a dipole fit to each ensemble (blue circles, orange downwards pointing triangles, and green triangles).The filled red asterisks and corresponding bands show the continuum limit using the two-step approach, i.e. via fits linear in a 2 to the axial charge and radius obtained at each value of a 2 .At a 2 = 0 we compare with results obtained from a one-step approach (open asterisks), as described in the text.
we depict the extracted Results for the axial charge and radius as a function of Q 2 max obtained from using the z 3 -expansion to parameterize the Q 2 -dependence.For each Q 2 max we depict five points having prior width w = 1, 2, 3, 4, 5.The points are shifted to the right as w increases with an increasing symbol size.
19G. 18.FIG.19.The axial charge (top) and radius (bottom) obtained via third-order z-expansion fits and model average for each ensemble.The notation is the same as that in Fig.13.

TABLE IX .
Results for the isovector axial charge and radius extracted using four different approaches as described in the text.The values are depicted in Fig.20.
2, we take a 2-state as our central values and the difference between the central values of a 2-state and a 3−state as the systematic error to account for systematics due to excited states.Our final FIG.21.Continuum limit of GA(Q 2 ) using the z 3 -expansion and data from the two-state fit analysis of the correlators up to Q 2 = 1 GeV 2 for the three ensembles with the symbols as indicated in the header of the figure.

TABLE XI .
Pion pole masses per ensemble extracted from the individual or combined fit of GP and G5 from two-or three-state fit data.three-statefitanalyses are in good agreement with each other.The values of g πN N extracted from G P (Q 2 ) and G5 (Q 2 ) are also in agreement within error bars.If we enforce the pion pole and g πN N to have the same value as determined from G P (Q 2 ) and G 5 (Q 2 ) and perform the same analysis, we obtain 27G.27.Results on GP (Q 2 ) for each ensemble (blue band for the cB211.72.64, orange band for the cC211.60.80 and green band for the cD211.54.96 ensemble) and at the continuum limit (red band) using the z 3 -expansion to fit the Q 2 -dependence of the data determined from the two-state fit analysis up to Q 2 = 1 GeV 2 .The inner panel shows a zoomin of the region marked by the red square. .
using the data for G w pole . 2.