Nucleon axial and pseudoscalar form factors from lattice QCD at the physical point

We compute the nucleon axial and induced pseudoscalar form factors using three ensembles of gauge configurations, generated with dynamical light quarks with mass tuned to approximately their physical value. One of the ensembles also includes the strange and charm quarks with their mass close to physical. The latter ensemble has large statistics and finer lattice spacing and it is used to obtain final results, while the other two are used for assessing volume effects. The pseudoscalar form factor is also computed using these ensembles. We examine the momentum dependence of these form factors as well as relations based on pion pole dominance and the partially conserved axial-vector current hypothesis.


I. INTRODUCTION
A central aim of on-going experimental and theoretical studies is the understanding of the structure of the proton and the neutron arising from the complex nature of the strong interactions. The electron scattering off protons is a well developed experimental approach used in such studies. An outcome of the multi-years experimental programs in major facilities has been the precise measurement of the electromagnetic form factors, see e.g. [1][2][3][4][5][6][7][8]. However, despite many years of experimental effort, new features are being revealed by performing new more precise experiments as, for example, the measurement of the proton charge radius [9][10][11]. Experimental efforts are accompanied by theoretical computations of such quantities [1,[12][13][14][15][16]. However, the theoretical extraction of such form factors is difficult due to their non-perturbative nature. The lattice formulation of Quantum Chromodynamics (QCD) provides the nonperturbative framework for computing non-perturbative quantities from first principles. Lattice QCD computations using simulations at physical parameters of the theory of Electromagnetic form factors is a major recent achievement [17][18][19][20][21].
While the electromagnetic form factors are well mea-sured and are being used to benchmark theoretical approaches, the nucleon axial form factors are less well known. The axial form factors are important quantities for weak interactions, neutrino scattering and parity violation experiments. Neutrinos can interact with nucleons via the neutral current of weak interactions, exchanging a Z 0 boson or via the charged current of weak interactions exchanging a W ± boson. The nucleon matrix element of the isovector axial-vector 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 [22][23][24], while G P (Q 2 ) from the longitudinal cross section in pion electro-production [25][26][27]. 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 [28][29][30][31]. 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 µ [32][33][34][35][36], where m µ is the muon mass.
Besides experimental extractions, phenomenological arXiv:2011.13342v1 [hep-lat] 26 Nov 2020 approaches are being applied to study the axial form factors. Chiral perturbation theory provides a nonperturbative framework suitable for low values of Q 2 up to about 0.4 GeV 2 [27,37,38]. Other models used include the perturbative chiral quark model [39], the chiral constituent quark model [40] and light-cone sum rules [41].
As already mentioned, lattice QCD provides the ab initio non-perturbative framework for computing such quantities using directly the QCD Lagrangian. Early studies of the nucleon axial form factors were done within the quenched approximation [42,43], as well as, using dynamical fermion simulations at heavier than physical pion masses [44]. Only recently, several groups are computing the axial form factors using simulations generated directly at the physical value of the pion mass [20,21,[45][46][47][48][49]. Such simulations at the physical pion mass can check important phenomenological 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 the pseudoscalar G 5 (Q 2 ) form factor. 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 this paper. The computation of the form factors is performed using one ensemble of mass degenerate up and down quarks, and a strange and a charm quark (N f = 2 + 1 + 1) with masses tuned to their physical values, referred to as physical point. In addition, we present results for two ensembles of N f = 2 light quarks tuned to the physical pion mass. They have the same lattice spacing a but different volumes in order to check for finite size effects. Final results are given for the N f = 2 + 1 + 1 ensemble where high statistics are used and systematic errors due to excited states are better controlled.
The remainder of this paper is organized as follows: In Section II we discuss the PCAC and PPD relations and in Sec. III the parameterization of the Q 2 dependence. In Sec. IV, we explain in detail the lattice methodology to extract the axial and pseudoscalar form factors. The renormalization of the operators is discussed in Sec. V. In Sec.VI, we explain how we extract the energy of the excited state and in Secs. VII and VIII we show results for the nucleon state matrix elements of the axial-vector and pseudoscalar currents. We compare our results of the three ensembles in Sec. IX and present the final results in Sec. X. A comparison with other studies is undertaken in Sec. XI. Finally, we conclude in Sec. XII.

II. DECOMPOSITION OF THE NUCLEON AXIAL-VECTOR AND PSEUDOSCALAR MATRIX ELEMENTS INTO THE FORM FACTORS AND THEIR RELATIONS
In this work we will consider the isovector axial-vector operator given by where u and d is the isospin double of the up and down quark fields. 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 MeV for the pion decay constant. the In QCD the axial Ward-Takahashi identity leads to the partial conservation of the axial-vector current (PCAC) where 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 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) momentum p(p ) and spin s(s ), q = p − p the momentum transfer and q 2 = −Q 2 . The nucleon pseudoscalar matrix element is given by where P 5 =ūγ 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. (8) is written so that it illustrates the pole structure of G 5 (Q 2 ). Substituting G 5 (Q 2 ) in Eq. (7), one obtains the Goldberger-Treiman relation [44,50] The pion-nucleon form factor G πN N (Q 2 ) at the pion pole gives the pion-nucleon coupling g πN N ≡ G πN N (Q 2 = −m 2 π ). In the limit Q 2 → −m 2 π , the pole on the right hand side of Eq. (9) must be compensated by a similar one in G P (Q 2 ), since G A (−m 2 π ) is finite. Therefore, if we multiply Eq. (9) by (Q 2 +m 2 π ) and take the limit towards the pion pole we have and, thus, one can extract g πN N from the induced pseudoscalar form factor too. Close to the pole, pion pole dominance means that G P (Q 2 ) = 4m N F π G πN N (Q 2 )/(m 2 π + Q 2 ). Inserting it in Eq. (9) we obtain the well known relation [51] which means that G P (Q 2 ) can be expressed as [52] From Eq. (11), the pion-nucleon coupling can be expressed as g πN N = m N G A (−m 2 π )/F π . In the chiral limit, lim mπ→0 G A (−m 2 π ) → g A and we have that The deviation from Eq. (13) due to the finite pion mass is known as the Goldberger-Treiman discrepancy, namely and it is estimated to be at the 2% level [53].

III. Q 2 -DEPENDENCE OF THE AXIAL AND PSEUDOSCALAR FORM FACTORS
For the parameterization of the Q 2 -dependence of the axial and pseudoscalar form factors typically two functional forms are employed, the dipole Ansatz and the model independent z-expansion [54,55].
The dipole Ansatz is given by with m the dipole mass. In the case of the axial form factor G A (Q 2 ), its value for Q 2 = 0, gives the axial charge g A ≡ G A (0) and the dipole mass m is the axial mass m A . Customarily, one characterizes the size of a hadron probed by a given current by the root mean square radius (r.m.s) defined as r 2 . The radius of the form factors can be extracted from their slope as Q 2 → 0, namely Combining Eq. (15) and Eq. (16) one can show that the radius is connected to the dipole mass as In the case of the z-expansion the form factor is expanded as, where imposing analyticity constrains, with t cut the particle production threshold. For t cut , we use the three-pion production threshold, namely t cut = (3m π ) 2 [55]. The coefficients a k should be bounded in size for the series to converge and convergence is demonstrated by increasing k max . Since the possible large values of the a k for k > 1 can lead to instabilities, we use Gaussian priors centered around zero with standard deviation w max(|a 0 |, |a 1 |) [56], where w controls the width of the prior. The value of the form factor at zero momentum is G(0) = a 0 , while the radius is given by In the case of the axial form factor, a 0 and a 1 should have opposite signs leading to positive radii. By comparing Eq. (20) to Eq. (17), we define the corresponding mass determined in the z-expansion to be In the case of G P (Q 2 ) and G 5 (Q 2 ), the pion pole is first factored out and thus (Q 2 +m 2 π ) G P,5 (Q 2 ) could be fitted using the dipole and z-expansion functions.

IV. LATTICE METHODOLOGY
In this section we describe the lattice QCD methodology to extract the form factors, presenting the construction of the appropriate three-and two-point correlation functions, the procedure to isolate the ground state and the details about the ensembles used.

A. Correlation functions
The extraction of the nucleon matrix elements involves the computation of both three-and two-point Euclidean correlation functions. The two-point function is given by (22) where with x 0 is the source and x s the sink positions on the lattice where states with the quantum numbers of the nucleon are created and destroyed, respectively. The interpolating field is where C is the charge conjugation matrix and Γ 0 is the unpolarized positive parity projector Γ 0 = 1 2 (1 + γ 0 ). By inserting the unity operator in Eq. (22) in the form of a sum over states of the QCD Hamiltonian only states with the quantum numbers of the nucleon survive. The overlap terms between the interpolating field and the nucleon state |N j as Ω|J N |N j are terms that 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 signal-to-noise ratio decays exponentially with the Euclidean time evolution. To accomplish ground state dominance, we apply Gaussian smearing [57,58] to the quark fields entering the interpolating field where the hopping matrix is given by The parameters a G and N G are tuned [17,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. For the construction of the three-point correlation function the current is inserted between the time of the creation and annihilation operators giving (26) where Γ k = iΓ 0 γ 5 γ k . The Euclidean momentum trasfer squared is given by Q 2 = −q 2 = −(p − p) 2 , and from now on we will use p = 0.

B. Treatment of excited states contamination
The interpolating field in Eq. (23) creates a tower of states with the quantum numbers of the nucleon. Gaussian smearing helps to reduce them but we still need to make sure that we extract the nucleon matrix element that we are interested in and that any contribution from nucleon excited states and/or multi-particle states are sufficiently suppressed.
In order to cancel the Euclidean time dependence of the three-point function and unknown overlaps of the interpolating field with the nucleon state, we construct an appropriate ratio of three-to a combination of two-point functions [61][62][63][64], Without loss of generality, we take t s and t ins relative to the source time t 0 , or equivalently t 0 is set to zero. The ratio in Eq. (27) is constructed such that in the limit of large time separations (t s − t ins ) a and t ins a, it converges to the nucleon ground state matrix element, namely How fast we ensure ground state dominance depends on the smearing procedure applied on the interpolating fields, as well as on the type of current entering the threepoint function. In order to check for ground state dominance we employ three methods as summarized below: where ∆E is the energy gap between the nucleon first excited state and the ground state. Assuming that the exponential terms in Eq. (29) are small we can extract the first term that gives the matrix element of interest by looking for a range of t ins for a given t s for which Eq. (27) is time-independent (plateau region) and fit to a constant (plateau value). We then increase t s until the plateau values converge. The converged plateau values determine the ground state nucleon matrix element of the current considered. Summation method: The insertion time, t ins , of the ratio in Eq. (27) can be summed leading to [65,66] R summ Although we also take into account only the lowest state, the contributions from excited states decay faster as compared to the plateau method. Since t ins is taken around t s /2 the summation method may be considered equivalent to the the plateau method with about twice t s . If e −∆Ets is sufficiently suppressed in Eq. (30) the slope gives the ground state matrix element. We probe convergence by increasing the lower value of t s , denoted by t low s entering in the linear fit. The disadvantage of the summation method is that one needs to do a linear fit with two parameters instead of one as for the plateau method. This leads to an increased statistical error. Two-state fit method: In this approach one considers explicitly the contribution of the first excited state. Namely, the two-point function is taken to be C( p, t s ) = c 0 ( p)e −E0( p)ts + c 1 ( p)e −E 2pt 1 ( p)ts (31) and the three-point function where contributions from states beyond the first excited state are neglected. As will be discussed in detail in the following sections, we allow the first excited state in the three-point function to be in general different from that of the two-point function. The coefficients of the exponential terms of the two-point function in Eq. (31) are overlap terms given by where spin indices are suppressed. The i-index denotes the i th nucleon state that may also include multi-particle states. The terms A i,j appearing in the three-point function in Eq. (32) are given by where N i ( 0)|A µ |N j ( p) is the matrix element between i th and j th nucleon states. Multi-particle states are volume suppressed and are typically not observed in the two-point function. However, if they couple strongly to a current they may contribute in the three-point function. As pointed out in Refs. [67,68], this may happen for the case of the axialvector current considered here. In order to include the possibility that multi-particle states contribute to the three-point function, we perform the following types of fits: M1: We assume that the first excited state is the same in both the two-and three-point functions. In this case, we first fit the two-point function extracting c 1 ( p) and E 1 ( p) and then use them when fitting the ratio of Eq. (27). We also fit the zero momentum two-point function to determine the nucleon mass and then use the continuum dispersion relation E 0 ( p) = m 2 N + p 2 to determine the nucleon energy for a given value of momentum. The continuum dispersion relation is satisfied for all the momenta considered in this work as can be seen in Fig. 2. We will refer to this as fit M1.
M2: We allow the first excited state to be different in the two-and three-point functions. In this case, the first excited energy of the three-point function is left as a fit parameter. We will refer to this as M2 fit.
In Fig. 1 we show the energies extracted from the nucleon two-point function as well as the two-particle noninteracting πN energies computed as the sum of the pion and nucleon energies. We show these energies for both the charged and neutral pions. As can been seen, the first excited state E 2pt 1 (p) extracted from two-point function coincide with that of the Roper resonance with at the same momentum. The lowest two-particle states are not visible in the two-point functions, although they are much lower than the energy of the Roper. This is expected since they volume suppressed. We note that the energies of the π + N and π 0 N system are consistent within errors. More details on these two fit approaches are given in Sec. VI.

C. Extraction of the axial and induced pseudoscalar form factors
While the pseudoscalar matrix elements lead directly to the G 5 (Q 2 ) as given in Eq. (B4), the matrix element of the axial-vector current in general contributes to both axial and induced pseudoscalar form factors, as given in Eqs. B1) and (B2. A procedure to extract the two form factors is to minimize χ 2 given by where w µ (Γ k , q; t s , t ins ) is the statistical error of the ratio R µ (Γ k , q; t s , t ins ) of Eq. (27) and F (Q 2 ; t s , t ins ) is a two component vector of the axial form factors The definition of the coefficient matrix G µ (Γ k , q) that has the kinematical factors is given in Eq.(B3). Minimization of the χ 2 defined in Eq. (35) is equivalent to a singular value decomposition (SVD), where U is a hermitian N × N matrix with N being the number of combinations of µ, k and components of q that contribute to the same Q 2 . Σ is the pseudo-diagonal N × 2 matrix of the singular values ofG and V is a hermitian 2 × 2 matrix since we have two form factors. Typically, N 2 for finite momenta. In our analysis, we use the SVD to extract the form factors since it does not need any minimization algorithm that might depend on the initial parameters. In addition, using the SVD approach for a relatively small matrix is much faster than using minimization algorithms.
In the following sections, results are presented for the ratios of G A and G P as described by Eq. (36).

D. Parameters of the gauge configuration ensembles
In this work we analyze an N f = 2 + 1 + 1 twisted mass clover-improved fermion ensemble. The parameters are given in Table I. In addition, we analyze two N f = 2 ensembles with the same light quark action, referred to as cA2.09.48 and cA2.09.64 ensembles. They have the same lattice spacing and two different volumes to check for finite volume effects. The physical volume of the cB211.072.64 ensemble is in between the volume of the two N f = 2 ensembles. Results on the axial form factors for the cA2.09.48 ensemble have been presented in Ref. [45] but are reanalyzed in this study and the results are used for the volume comparison. For both N f = 2 + 1 + 1 and N f = 2 ensembles the lattice spacing is determined using the nucleon mass. More details on the lattice spacing determination are given in Refs. [17,59,69,70].
The gauge configurations were produced by the Ex-tended Twisted Mass Collaboration (ETMC) using the twisted mass fermion formulation [71,72] with a clover term [73] and the Iwasaki [74] improved gauge action.
Since the simulations were carried out at maximal twist, we have automatic O(a) improvement [71,72] for the physical observables considered in this work.  [69] and the two N f = 2 ensembles, cA2.09.48 [75] and cA2.09.64. cSW is the value of the clover coefficient, β = 6/g where g is the bare coupling constant, N f is the number of dynamical quark flavors in the simulation, a is the lattice spacing, V the lattice volume in lattice units, mπ the pion mass, mN the nucleon mass, and L the spatial lattice length in physical units.

E. Three-point functions and statistics
Since in this work we study only isovector combinations, only connected contributions are needed. For their evaluation we employ standard techniques, namely the so-called fixed-sink method using sequential propagators through the sink. In this method, changing the sinksource time separation t s , the momentum, the projector or the interpolating field at the sink requires a new sequential inversion. We, thus, fix the sink momentum p = 0 and use four projectors, namely the unpolarized Γ 0 and the three polarized projectors Γ k . In the case of the cB211.072.64 ensemble, we perform the analysis using in total seven sink-source time separations, t s , in the range 0.64 fm to 1.60 fm. In order to better isolate the contribution from excited states, we need to compute the three-point functions at similar statistical accuracy. However, the signal-to-noise ratio drops rapidly with t s and, thus, we increase statistics as t s increases, keeping approximately the statistical error constant. The number of configurations analyzed for the N f = 2+1+1 ensemble is kept at 750 for all values of t s . Statistics are increased by increasing the number of source positions per gauge configuration after checking that the error continues to scale as expected for independent measurements. The statistics used for all the three ensembles for the computation of the connected contribution per t s are shown in Table II.

V. RENORMALIZATION FUNCTIONS
Matrix elements computed in lattice QCD need to be renormalized in order to extract physical observables. The renormalization functions, or Z-factors, for the N f = 2 ensembles have been computed previously [45]. A detailed description about our procedure can be found in Ref. [76]. Here we present a summary on the evaluation of the Z-factors for the N f = 2 + 1 + 1 cB211.072.64 ensemble. For this work in the twisted mass formulation, we need the renormalization functions Z S used for the renormalization of the pseudoscalar form factor G 5 (Q 2 ), Z P used for the renormalization of the bare quark mass and Z A used to renormalize the axial-vector current.
We employ the Rome-Southampton method or the socalled RI scheme [77], and compute the quark propagators and vertex functions non-perturbatively. This scheme is mass-independent, and therefore, the Z-factors do not depend on the quark mass. However, there might be residual cut-off effects of the form a 2 m 2 q [78] and, for the scale dependent renormalization functions Z S and Z P , the RI-MOM Green functions have also a dependence on m 2 q /µ 2 . This is why the RI-MOM renormalization functions must be explicitly defined in the chiral limit. If not, the scheme would not be mass-independent. To eliminate any systematic related to such effects, we extract the Z-factors using multiple degenerate-quark ensembles. We use five N f = 4 ensembles generated exclusively for the renormalization program at the same β value as that of the cB211.072.64 ensemble. These are generated with quark mass which is less than half of the strange mass, in order to suppress the m 2 q /p 2 for to scale-dependent renormalization functions, and the lattice artifacts O(a 2 m 2 q ). These ensembles are generated at different pion masses in the range of [366-519] MeV and a lattice volume of 24 3 × 48 in lattice units. Having five pion masses enables us to perform the chiral extrapolation to eliminate from the Z-factors any residual cut-off effects. It should be noted, that the extrapolation in the Z-factors does not have an impact on the nucleon matrix elements, which are calculated directly at the physical point.
In this study, we use the operators written in the twisted (χ,χ) and physical basis (ψ,ψ) with ψ and χ the u and d doublet and τ b are the three Pauli matrices. In the chiral limit, the renormalization functions become independent of the isospin index b, and can be dropped. We use the combinationūΓd, which is extracted fromτ ≡ τ 1 +iτ 2
We note that the PCAC relation in the twisted basis is given by where the axial-vector current A b µ = Z Aχ γ µ γ 5 τ b χ, the pseudoscalar operator P b = Z Pχ γ 5 τ b χ and the scalar S 0 = Z Sχ χ. For the isovector flavor combination b = 3 and at maximal twist where the PCAC mass m PCAC is tuned to zero, Eq. (44) reduces to where m q the renormalized quark mass determined from the twisted light quark mass parameter µ as m q = µ/Z P . The aforementioned operators are renormalized multiplicatively with Z O , using the condition where S L (p) and Γ L (p) are the quark propagator and amputated vertex function, respectively, while S Born (p) and Γ Born are their tree-level values. The trace is taken over spin and color indices and the momentum p is set to be the same as the RI renormalization scale µ 0 . For the non-perturbative calculation of the vertex functions we use momentum sources [79] that allow us to reach per mil statistical accuracy with O(10) configurations [80,81]. High statistical precision means that one has to sufficiently suppress systematic errors. We choose momenta in a democratic manner, namely (a p) ≡ 2π where n t ∈ [2,10], n x ∈ [2, 5] and T /a(L/a) the temporal(spatial) lattice extent. The momenta are chosen in the aforementioned ranges with the constraint [78] to suppress non-Lorentz invariant contributions. These constraints are chosen to suppress O(a 2 ) terms in the perturbative expansion of the Green's function and are expected to have non-negligible contributions from higher order in perturbation theory [76,80,81]. We subtract such finite lattice spacing effects by explicitly computing such unwanted contributions to one-loop in perturbation theory and all orders in the lattice spacing. These finite a artifacts appear in both the S L (p) and Γ L (p) functions. This improvement of non-perturbative estimates using perturbation theory, significantly improves our estimates, as can be seen in the plots of this section.
Let us first discuss our results on Z A , which is scheme and scale independent. In order to eliminate cut-off effects in Z A , we perform a linear fit with respect to (am π ) 2 (equivalently am q ), for every value of the renormalization scale. In Fig. 3 we show the mass dependence for a specific value of the RI scale. We find a slope that is compatible with zero, as expected from our previous studies [76]. In order to eliminate the residual dependence on the initial scale due to lattice artifacts, we perform an extrapolation to (aµ 0 ) 2 → 0. In Fig. 4, we show the linear extrapolation in (aµ 0 ) 2 . In the plot we show the purely non-perturbative values of Z A , as well as the improved values obtained after subtracting the lattice artifacts calculated perturbatively. Such a subtraction procedure improves greatly the estimates for Z-factors, as it captures the bulk of lattice artifacts. Indeed, a linear fit in (aµ 0 ) 2 in the improved subtracted data yields a slope that is consistent with zero within uncertainties. The Z P and Z S renormalization factors are scheme and scale dependent. Therefore, after the extrapolation (am π ) 2 → 0, we convert to the MS-scheme, which is commonly used in experimental and phenomenological studies. The conversion procedure is applied on the Zfactors at each initial RI scale (a µ 0 ), with a simultaneous evolution to a MS scale, chosen to be µ=2 GeV. For the conversion and evolution we employ the intermediate Renormalization Group Invariant (RGI) scheme, which is scale independent and connects the Z-factors between the two schemes: with O = P, S. Therefore, the appropriate conversion factor to multiply Z RI O is .
(50) The quantity ∆Z S O (µ 0 ) is expressed in terms of the β-function and the anomalous dimension, γ S , of the operator under study and may be expanded to all orders of the coupling constant. The superscript S denotes the scheme of choice. The expressions for the scalar and pseudoscalar operators are known to three-loops in perturbation theory and can be found in Ref. [76] and references therein. In Fig. 5 we present our results for Z P and Z S . We collect our results for the renormalization functions in Table. III. We note that the errors given are statistical. A full analysis of systematic errors is ongoing and will be presented in an upcoming publication. It is expected that systematic errors will mostly affect the errors on Z P and Z S and will not have any significant effect on the results presented here.

VI. EXTRACTION OF EXCITED ENERGIES
In this section we discuss the details for the identification of the nucleon matrix elements. As mentioned in Sec. IV B, we apply two procedures, referred to as M1 and M2. Our fit procedure is illustrated for the case of the N f = 2 + 1 + 1 cB211.072.64 ensemble but the same procedure is carried out for the two N f = 2 ensembles. The two-state M1 fit has been used in previous analyses of form factors, including G P (Q 2 ). However, as pointed out in Ref. [67], the πN state that is suppressed in the two-point function may become dominant in the threepoint function in the case of G P (Q 2 ) that is dominated by the pion pole at low Q 2 values. Therefore, we allow the energy of the first excited state to be different in the two-and three-point functions, as done in the type M2 fit. As suggested in Ref. [46], one can use the temporal component of the axial vector current, A 0 , which is very precise, in order to determine the first excited energy. The temporal component has not been used in past studies [20,45,47,56,82,83], since it has been found to suffer from large excited state contributions.
In Figs. 6 and 7 we show, respectively, the results when using the two-state M1 and M2 fit types. We use the ratio constructed with the three-point function of the temporal axial-vector current. We perform a simultaneous fit on several sink-source time separations, t s , excluding the three smallest t s to ensure no contamination from higher excited states. As can be seen, the M2 fit describes better the data as reflected by the better χ 2 /d.o.f.
In Fig. 8 we show the energy of the first excited state extracted from fitting the two-point and the three-point function of the temporal axial-vector current. We observe that the first excited energy as extracted from the two-point function is in agreement with the energy of the Roper. This is a different behavior from what is observed in the two recent studies [46,49], where the first excited state extracted from the two-point function is much higher. Moreover, the energy of the first excited state extracted from the three-point function, is in general in agreement with the energy of the non-interacting two-particle states of N (0) + π(− p) and N ( p) + π(− p).  We do not observe states with energies lower than the non-interacting state energies unlike what was found in Ref. [46]. in particular the effect of the excited states. We first consider the pseudoscalar matrix element, since it is only connected to one form factor, G 5 (Q 2 ), as described in Eq. (B4) and thus the simplest to extract.
For the identification of the nucleon matrix element we apply the three approaches discussed in Sec. IV B in order to analyze contributions from excited states. In Figs. 9 and 10, we demonstrate how excited state contributions are identified for the two smallest Q 2 . In particular, we show the ratio of Eq. (27) for all the available values of t s . In the construction of the ratio we use two-point functions computed at the same source positions as the corresponding three-point functions to exploit their correlation that results in a reduction in the the error. As t s increases we see a significant increase in the values of the ratio pointing to a sizeable excited states contamination. In the same figure, we show also the plateau values for the two largest time separations obtained by discarding 7 time slices from source and sink or the midpoint (t ins = t s /2) of the ratio for the smaller time separations.
In Fig. 9, we include results when using both the M1 and M2 fits for the two-state approach as well as the summation method. We note that, while for the G P (Q 2 ) form factor there is a chiral perturbation theory support for the dominance of the lowest πN state in the threepoint function [67], such a theoretical argument is not presented in the case of G 5 (Q 2 ). However, we empirically use the M2 fit in order to examine if the excited states in the lattice data can be described with such a fit function. As can been seen, both M1 and M2 describe well the data with M2 providing a better agreement for the larger values of t s . Increasing t low s does not change the results extracted from the two-state fits, which shows that including an excited state captures well the time dependence of the ratio. This is unlike the summation method, for which we observe an increase with increasing t low s . We use as our final values the one determined from the two-state fit at a value of t low s that is consistent with the summation values in some range. The final value is larger for the case of M2. This is expected since for these momentum transfer the exited energy extracted from the three-point function is lower as compared to the one extracted from the two-point function. However, this increase is not as large as observed in studies of Refs. [46,48]. Comparing the behavior of the excited states at the second smallest Q 2 in Fig. 10 we find the same conclusions as for the lowest Q 2 value. In both cases our final value is the one from the two state fit at t low s =1.12 fm as discussed in Sec. VI. This is what we use for all the Q 2 values. In this section we discuss the analysis of the correlators of the axial-vector current from which the axial and induced pseudoscalar form factors are extracted. For the determination of the two form factors we follow the procedure discussed in Sec. IV C. As explained in Secs. VI and VII, the dominance of two-particle states is expected only to enter the determination of the induced pseudoscalar form factor. For G A (Q 2 ) no such strong coupling is expected. Therefore, only the M1 fit is applied for the extraction of G A (Q 2 ).
In Fig. 11, we present the analysis of the effect of excited states for the ratio leading to the axial form factor G A (Q 2 ). We show results at the smallest Q 2 value and at some intermediate Q 2 value to give the general behavior as Q 2 increases. We observe that there is a faster convergence as compared to the case of G 5 (Q 2 ). It is interesting that, while for the smaller values of Q 2 the effect of suppressing excited states is to increase the value of G A (Q 2 ), for higher momenta we find that the effect is to decrease it. Comparing the values of G A (Q 2 ) extracted from the summation and the two-state fits, we find agreement.
In Fig. 12, we present the analysis of the effect of ex-cited states for the ratio for the induced pseudoscalar FIG. 11: Excited states analysis for the ratio of the three-point correlator for the extraction of GA(Q 2 ) for Q 2 = 0.057 GeV 2 (top) and Q 2 = 0.271 GeV 2 (bottom). The notation is the same as that in Fig. 9. For the middle panel, the plateau values are used. The two-state fit analysis is done only with the type M1 fit. In this case we use t low s /a = 8 because it does not suffer from the issues discussed for G5(Q 2 ).
form factor for the smallest Q 2 . What we observe is that the effect of excited states is similar to what is observed in the analysis of G 5 (Q 2 ) in Fig. 9. G 5 (Q 2 ) and G P (Q 2 ) have the same pion pole behavior and therefore such sim-ilarities are expected. As in the case of G 5 (Q 2 ) we carry out the M2 fit in addition to M1. The notation is the same as that in Fig. 9.

IX. COMPARISON OF RESULTS USING THE THREE ENSEMBLES
We perform a similar analysis as for the N f = 2 + 1 + 1 cB211,072.64 ensemble also for the two N f = 2 ensembles. In Fig. 13, we compare results from the three ensembles for G A (Q 2 ). In particular, comparing the results between the two N f = 2 ensembles we do not observe any finite volume effects in the range m π L ∈ [3,4]. In The M1 approach has been used for this case. and 15 we compare our results for G 5 (Q 2 ) and G P (Q 2 ) for the three ensembles, using the M2 fit. We observe a very good agreement among the results for the three ensembles. Like for the case of G A (Q 2 ), comparison between the results of the two N f = 2 ensembles does not show any finite volume effects in the range m π L ∈ [3,4].
In Figs. 16 and 17 we show a comparison between the M1 and M2 fits for G 5 (Q 2 ) and G P (Q 2 ). For G 5 (Q 2 ) we include the prediction when using the PCAC and PPD relations given by Eqs. (7) and (12) Comparing the results extracted using the two-state M1 to M2 fits to extract the nucleon matrix elements, we find that the latter approach yields higher values for Q 2 < 0.2 GeV 2 . Despite the increase, however, results for G 5 (Q 2 ) predicted from PCAC deviate significantly in the low Q 2 region from those extracted directly from the nucleon matrix element of the pseudoscalar operator, in contrast to what has been observed in Refs. [46,49]. This different behavior can be traced to the fact that the authors of Refs. [46,49] find a higher energy for the first excited state from their two-point functions as compared to us. Also the energy of the first excited state extracted from the three-point function of the temporal axial-vector current in Ref. [46] is lower than what we find and lower than the corresponding non-interacting energy.
The authors of Ref. [49] on the other hand find an exited state that is closer to the non-interacting energy as we do, although a direct comparison is not possible since only results for a heavier than physical pion mass are shown. These observations also hold for G P (Q 2 ), as shown in Fig. 15. It is interesting to examine the breaking of the PCAC and PPD relations as a function of Q 2 . We define two ratios, one checking the PCAC and one the PPD relation as follows and These ratios are unity if PCAC and PPD hold, respectively. In Fig. 18, we concentrate on the results for the The notation is as in Fig. 16. Open red circles are results from the pion pole dominance prediction. cB211.072.64 ensemble since the results using the other two ensembles behave similarly. As can be seen there is a sizeable deviation for both ratios at small Q 2 even though we use the M2 fit. In our view, further investigation is needed to understand the deviations from the PCAC and PPD relations. Therefore, in what follows we use the results of G A (Q 2 ) to extract both G P (Q 2 ) and G 5 (Q 2 ) from Eqs. (12) and (52). Also we only discuss our results extracted using the cB211.072.64 ensemble, since they are more precise and are computed for a lattice volume that is in between the two lattice volumes used for checking for volume effects, for which we see no effects.

X. RESULTS
All results given in this section are extracted using the N f = 2 + 1 + 1 cB211.072.64 ensemble. In Fig. 19 we show our results for the axial form factor. The value of the form factor at zero momentum transfer gives the axial charge, g A ≡ G A (0). We find g A = 1.283 (22) [85]. In order to fit the form factor, we use both the dipole and z-expansion (see Sec. III for details). Since for G A (Q 2 ) the value for zero momentum transfer is directly accessible, we use G A (0) in the jackknife fits. This reduces the number of fit parameters in each jackknife bin. The consequence is that the error on the determination of the radius is smaller. In the case of the z-expansion, we use k max = 5, where we check that this is large enough to ensure convergence. The width coefficient, w, of the Gaussian priors is chosen to be w = 5. We provide a systematic error taken as the difference in the mean values when using w = 5 and when w = 20. Comparing the dipole Ansatz with the z-expansion we find excellent agreement for all Q 2 values. Therefore, we conservatively quote as final values (Table IV) those from the z-expansion, since it is model independent although they typically carry larger statistical uncertainties. The axial mass and the radius are determined from the parameters of the z-expansion as given in Eqs. (21 and 20) correspondingly.
In Fig. 20 we show our results for the induced pseudoscalar form factor extracted using G A (Q 2 ) and the PPD of Eq. (12). The induced pseudoscalar coupling determined at the muon capture [86] is determined as with m µ = 105.6 MeV the muon mass. The pion-nucleon coupling constant given in Eq. (10), and the Goldberger- Treiman discrepancy given in Eq. (14) can also be extracted from the induced pseudoscalar form factor. We tabulate the extracted values in Table IV. The error on both g * P and g πN N due to using a different fit Ansatz as well as the maximum value of Q 2 used in the fits is negligible compared to the statistical error. The Goldberger-Treiman discrepancy is determined to high precision since it uses the precise values of the axial form factor.
Finally, our results for the pseudoscalar form factor G 5 (Q 2 ) extracted using Eq. (52) are shown in Fig. 21. In principle, the pion nucleon coupling can be extracted from this form factor but since we use the PCAC and PPD relations for both G P (Q 2 ) and G 5 (Q 2 ), one would obtain the same value as the one extracted from the G P (Q 2 ) form factor.
We tabulate our results for the three form factors as a function of Q 2 in the Appendix A. A , and the r.m.s r 2 A , the induced pseudoscalar coupling determined at the muon capture [86], the pion nucleon coupling as in Eq. (10), and the Goldberger-Treiman discrepancy as in Eq. (14). The first error is statistical and the second a systematic taken as the difference in the mean values when using w = 5 and w = 20. mA

XI. COMPARISON WITH OTHER STUDIES
While there are a number of lattice QCD studies on the isovector axial and pseudoscalar form factors using simulation with heavier than physical pion masses, we restrict our comparison here with results obtained using ensembles at the physical point. We summarize below the setup used by other groups to compute the isovector axial and pseudoscalar form factors: • The PNDME collaboration [46] used a hybrid action with N f = 2+1+1 HISQ configurations generated by the MILC collaboration with lattice spacing a 0.0871 fm, lattice volume 64 3 × 128 and m π = 130 MeV in the sea (referred as a09m130W) and clover improved valence quarks with m π = 138 MeV. Three-point functions were computed from three sink-source time separations in the range of [1-1.4] fm. They performed the two-state analysis using both the M1 and M2 fits discussed in Sec. IV B. In what follows we show their results extracted using the M2 fit since they considered them as their final values (referred in their work as S A4 type fit). No improvement of the currents used is discussed in order to eliminate O(a) cut-off artifacts, which would imply that they have larger finite lattice spacing effects as compared to our formulation.
• The RQCD collaboration [49], analyzed 37 CLS ensembles, but only two of these were simulated using physical pion masses. The ensembles were generated using a tree-level Symanzik-improved gauge action and N f = 2 + 1 clover-improved fermions.  Table IV, on the other hand, will be done for all three ensembles.
In Fig. 22, we compare our results for G A (Q 2 ) using the N f = 2 + 1 + 1 ensemble with the aforementioned lattice QCD studies. Overall, there is a very good agreement among all results, which indicates that lattice artifacts are small. PACS results [20] are available for very small Q 2 values since their lattice spatial extent is approximately twice as compared to the other lattices for which we show results.
In Fig. 23, we compare results for the isovector induced pseudoscalar form factor. The results from PACS are extracted using the plateau method at their largest time separation. The results from the PNDME and RQCD collaborations, were extracted using a two-state M2 fit. Results from this work using the cB211.072.64 ensemble are shown with red circles, from the PNDME collaboration [46] with green squares, from the RQCD collaboration [49] with blue upward-pointing triangles and from the PACS collaboration [20] with brown down-pointing triangles. Our results are determined using G A (Q 2 ) and Eq. (12) and are in agreement with those from the PNDME and RQCD collaborations. While results from PACS are lower than the others at small Q 2 values, their G P (Q 2 ) has been determined using the plateau fits at relatively small value of the the source-sink separations. Their values are higher as compared to what we find at the same time separation for the direct extraction of the G P (Q 2 ). This is something that needs to be further investigated.
In Fig. 24, we compare results for G 5 (Q 2 ). Results from PNDME are omitted since they show only bare results and no renormalization factor is provided. Results from RQCD are omitted because they give only results multiplied by m q /m N and they do not provide the renormalized value of m q . Comparing our results with those from PACS we observe agreement. This is interesting since the PACS results are extracted using the plateau method at a relatively small source-sink time separation. However, their results, unlike what we find directly from the three-point function of the pseudoscalar current using the M2 fit, show the correct pion pole behavior. Whether the reason is because they use a large volume has to be further investigated. We plan to do such a comparison in the future when an ensemble using a large volume becomes available.
In Fig. 25, we compare our results for the isovector m A and r 2 A with results from other lattice QCD studies and with phenomenological analyses using experimental data. Our results from the three ensembles are in agreement with those using the N f = 2 + 1 + 1 ensemble being the most precise. That value of m A agrees with the value reported by the MiniBooNE collaboration [87] as well as the one from the MINOS Near detector [88] and Ref. [23]. Comparing with other lattice QCD results we find that our values are compatible with the ones from the PACS and RQCD collaborations.
We compare our values on muon capture coupling constant, g * P , pion-nucleon coupling g πN N and the Goldberger-Treiman discrepancy, ∆ GT , with other lattice QCD groups, experimental results and phenomenology in Figs. 26 and 27. Our results using the three ensembles are in agreement with the values from the N f = 2 + 1 + 1 ensemble being the most precise. They are also in agreement with other lattice QCD results, although the errors on some lattice QCD results are large. Phenomenological results are in general much more precise for g πN N and ∆ GT . On the other hand, experimental results on g * P from ordinary muon capture are compatible with lattice QCD results but carry large errors, while the result from chiral perturbation theory [89], is as precise 0.6 0. 8 [46] (blue leftpointing triangle), from the RQCD collaboration [49] (purple right-pointing triangle) when using the z-expansion, and from the PACS collaboration [20] (brown rhombus). Inner error bars are statistical errors while outer errors bars include systematic errors. The black crosses are results from phenomenology. From top to bottom we show results from the MiniBooNE experiment using charged-current muon neutrino scattering events [87], from νµ-iron interactions using the MINOS Near Detector [88], from Ref. [23] using world data from neutrino-deuteron scattering and the z-expansion for the fit, and two very accurate results from world averages, one is from (quasi)elastic neutrino and anti-neutrino scattering experiments [89] and the other from charged pion electroproduction experiments [89]. as our value from the cB211.072.64.  The results for g * P . The notation for the lattice QCD results is the same as that in Fig. 25. Black crosses are results from experimental analyses for ordinary muon capture from Refs. [86,[89][90][91][92] and the precise result at the top of the figure is from chiral perturbation theory [89].
In Fig. 27 we compare our results for g πN N and ∆ GT . The only other lattice QCD results on g πN N and ∆ GT are from the RQCD collaboration [49]. As can be seen, our value for g πN N has smaller error since it is determined from Eq. (13) unlike the value by RQCD that does not use g A but instead uses the G P (Q 2 ) form factor and Eq. (10). Analyses of experimental results of pionnucleon scattering yield very precise values. We can determine ∆ GT precisely, extracting a value that is in agreement with the one obtained from the recent analysis of π−N elastic scattering data [53]. Results using QCD sum rules [93], heavy baryon chiral perturbation theory [94] and an older analysis of experimental data [95] are spread around our value. The notation for the lattice QCD results is the same as that in Fig. 25. We also show phenomenological results with the black symbols. For gπNN , these are taken from Refs. [96][97][98][99] and are results from analyses of experimental data on pionnucleon scattering cross-sections. For the case of ∆GT , these are from Refs. [53,95] , from baryonic QCD sum rules [93], and from heavy baryon chiral perturbation theory [94].

XII. CONCLUSIONS
Results on the axial and pseudoscalar form factors are presented using an N f = 2 + 1 + 1 ensemble directly at the physical point avoiding chiral extrapolation that may introduce uncontrolled systematic errors in the nucleon sector. Using N f = 2 ensembles with spatial extent 4.5 fm and 6 fm no detectable finite volume effects are observed within the range of these two volumes. Given that the analysis of the N f = 2 + 1 + 1 ensemble uses more statistics and allows for a better investigation of excited states effects, we quote as our final results those obtained using the N f = 2 + 1 + 1 ensemble.
Our results for the axial form factor, G A (Q 2 ), are the most accurate compared to those from other recent lattice QCD studies. The axial charge G A (0) ≡ g A is in agreement with the experimental value. Fitting the Q 2dependence of G A (Q 2 ), we extract precisely the axial mass m A and r.m.s radius given in Table IV. Our value for m A agrees with the value reported by the MiniBooNE collaboration [87] as well as the one from the MINOS Near detector experiment [88] and Ref. [23].
The analysis of the lattice data that yield the induced pseudoscalar G P (Q 2 ) and pseudoscalar G 5 (Q 2 ) form factors is performed using two approaches. In the first approach we take the excited energies extracted from the nucleon two-point function to coincide with those entering the three-point correlators, and in the second, we allow them to be different. While we obtain different excited energies from the three-point correlators, the difference is not as large as observed in two recent studies [46,49]. The reason is that the first excited state extracted from our two-point function is already lower as compared to what these other two studies find. The consequence is that the effect on the low Q 2 -dependence is smaller and, thus, the G P (Q 2 ) and G 5 (Q 2 ) do not fulfill the PCAC and the pion-pole relations. It is interesting to note that the analysis by the PACS collaboration that uses a significantly larger volume but extracts G P (Q 2 ) assuming ground state dominance (via the plateau approach), finds almost agreement with pion pole dominance. Therefore, in our view, further investigation is needed to settle the pion dominance behavior of both G P (Q 2 ) and G 5 (Q 2 ). In the future, we plan to perform an analysis on a larger twisted mass ensemble of spatial extent ∼ 7.7 fm, which is currently under production by ETMC.
Using the axial form factor G A (Q 2 ) and PCAC and pion-pole dominance, we extract the values of the pion nucleon coupling constant g πN N , the Goldberger-Treiman deviation from chiral symmetry ∆ GT and the muon capture coupling constant g * p , all of which are in agreement with other recent lattice QCD studies, with our results being more accurate. These are also consistent with phenomenological extractions. This agreement is a success of lattice QCD in being now in a good position to compute from first principles these quantities.
For the case of the pseudoscalar matrix element we have In the above expressions, E is the energy and m the mass of the nucleon. The kinematic factor C is given by