Quark flavor decomposition of the nucleon axial form factors

We present results on the isoscalar form factors including the disconnected contributions, as well as on the strange and charm quark form factors. Using previous results on the isovector form factors, we determine the flavor decomposition of the nucleon axial form factors. These are computed using an ensemble of $N_f=2+1+1$ twisted mass fermions simulated with physical values of quark masses. We investigate the SU(3) flavor symmetry and show that there is up to 10\% breaking for the axial and up to 50\% for the induced pseudoscalar form factors. By fitting the $Q^2$-dependence, we determined the corresponding root mean square radii. The pseudoscalar coupling of the $\eta$ meson and the nucleon is found to be $g_{\eta NN}=3.7(1.0)(0.7)$, and the Goldberger-Treiman discrepancy for the octet combination about 50\%.


I. INTRODUCTION
Axial form factors play a key role in the interactions of nucleons with the W and Z bosons, the carriers of the weak force. They also provide insights into the structure of the nucleon that in turn can affect our ability to compute cross sections that may aid us into revealing new physics. Neutron beta decay and other charged current weak interaction processes like ν µ + n → p + µ − are sensitive to the isovector axial form factor G u−d A (Q 2 ). Neutrino elastic scattering on protons is sensitive to the strange axial form factor of the proton G s A (Q 2 ), which for Q 2 = 0 determines the strange quark contribution to the proton spin ∆s. The role of strange quarks is also important for calculating the cross sections for a class of popular cold dark matter candidates [1]. A variety of experiments ranging from nuclear recoil direct-detection experiments to collider indirect-experiments are searching for dark matter candidates that use as input either spin-dependent or spin-independent nucleon cross sections. A first measurement of parityviolating asymmetries in forward elastic electron-proton scattering by HAPPEx [2] combined with data from neutrino and antineutrino-proton elastic scattering cross sections from Brookhaven E734 [3] determined simultaneously the strange vector and axial form factors of the proton at non-zero momentum transfer square Q 2 [4]. Additional parityviolating data from the G0 experiments [5,6] improved the determinations of the strange axial form factors [7]. The MicroBooNE neutrino detector at Fermilab aims to extract the strange axial form factor of the nucleon in the range of momentum transfers of 1 GeV 2 to as low as 0.08 GeV 2 [8,9]. Combining neutrino-proton neutral and charged current scattering cross section measurements with available polarized electron-proton/deuterium cross section data is expected to reduce the experimental uncertainty and allow for the extraction of ∆s with an order of magnitude better accuracy complementing polarized deep inelastic scattering experiments. The axial form factors are the main source of error in the description of neutrino-nucleon interactions. Therefore, a calculation of these form factors within lattice QCD will provide valuable input in experiments such as DUNE [10,11] and Hyper-K [12,13].
Lattice QCD provides the ab initio non-perturbative framework for computing the nucleon axial form factors using directly the QCD Lagrangian. While there are a number of lattice QCD studies of the isovector axial form factors with recent results given in Refs. [14][15][16][17][18][19], only a few studies are done for other flavor combinations [18][19][20]. The reason for this is that the isovector flavor combination is free of quark disconnected contributions. In Ref. [14] we presented our results for the isovector axial form factors, while also investigating finite volume effects. This work focuses on the arXiv:2106.13468v1 [hep-lat] 25 Jun 2021 study of the isoscalar, octet and singlet flavor combination by computing all disconnected contributions, allowing us to perform a flavor decomposition. The computation is performed using one ensemble of N f = 2 + 1 + 1 dynamical quarks with the up, down, strange and charm masses tuned to their physical values, referred to as physical point.
The remainder of this paper is organized as follows: In Section II we discuss the PCAC relation and the parameterization of the Q 2 dependence of the form factors. In Section III we explain in detail the lattice methodology to extract the nucleon axial and induced pseudoscalar form factors. In Section IV we discuss the renormalization and in Section V we show results for the isoscalar combination, where both connected and quark disconnected contributions are presented. The strange and charm form factors are presented in Section VI, in Section VII the flavor singlet and octet combinations are discussed and in Section VIII we provide the results for the form factors for each quark flavor. Final results are quoted in Section IX, comparisons with previous studies are carried out in Section X and in Section XI we conclude.

II. MATRIX ELEMENTS, FORM FACTORS AND Q 2 -DEPENDENCE
In a previous paper [14], we presented results on the isovector axial form factors G u−d A (Q 2 ) and G u−d P (Q 2 ), as well as, the pseudoscalar G u−d 5 (Q 2 ). We refer the reader to that paper for details on the computation of the isovector combination. In this paper, we will describe the flavor combinations where disconnected contributions are involved, such as the isoscalar combination A u+d =ūγ µ γ 5 u +dγ µ γ 5 d. ( Combining the isovector and isoscalar matrix elements one can extract the axial form factors for the up and down quarks. We will also compute the strange and charm form factors and construct SU(3) flavor combinations. Considering the u, d and s flavor triplet we form the flavor singlet combination given by and the flavor octet given by In the SU(3) flavor symmetric limit, the matrix elements of A 8 µ will only have connected contributions. The axial Ward-Takahashi identity that leads to the partial conservation of the axial-vector current (PCAC) is where P 8 is the octet pseudoscalar density. The octet combination of the induced pseudoscalar form factor G u+d−2s P is related to the pseudoscalar coupling between the η−meson and the nucleon. The isosinglet flavor combination, on the other hand, has an anomalous term [21] and it satisfies a modified relation, where P 0 =ūγ 5 u +dγ 5 d +sγ 5 s is the isosinglet pseudoscalar current, Q(x) is the topological density Q(x) = 1 32π 2 µνρσ Tr[F µν (x)F ρσ (x)] and F µν is the field strength tensor of QCD. The anomalous gluonic term is induced by the axial anomaly. Since gluons couple equally to each quark flavor, the anomalous term vanishes only for non-singlet combinations as in Eqs. (4). The anomaly term has the consequence that the axial-vector flavor singlet current is not conserved even for massless quarks. The nucleon matrix element of the axial operators in Eqs. (1), (3) and (2) 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 expression is given in Euclidean space. Note that we have suppressed the index denoting the flavor combination for simplicity. Calculations on the lattice allow us to compute the form factors only at given discrete values of Q 2 . In order to investigate their full Q 2 dependence we use two fit forms, the well known dipole Ansatz, and the model independent z-expansion [22,23]. In the case of the dipole Ansatz we have that In the case of the axial form factor, G(0) gives the axial charge and m the axial mass for the flavor combinations under investigation. The radius is extracted from the slope in the limit Q 2 → 0, namely Combining Eq. (7) and Eq. (8) one can show that the radius is connected to the dipole mass as 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 . In the case of the z-expansion, the form factor is expanded in a series as, where imposing analyticity constrains, with t cut the particle production threshold. We use the three-pion cut, t cut = 9m 2 π [23] for all flavor combinations, although apart from the isovector case, the cut-off might be higher due to heavier decay modes. The coefficients a k appearing in Eq. (10) should have an upper bound, so that the series converges at some value of k. Larger values of a k that could appear for k > 1 can lead to instabilities. Therefore, we employ Gaussian priors, which are centered around zero with a chosen standard deviation w max(|a 0 |, |a 1 |) [19], and with w controlling the width of the prior. The value of the form factor at zero momentum is a 0 and the radius is The coefficients a 0 and a 1 are anticipated to have opposite signs in order to lead to positive values of the radii. If we compare the above equation with the one extracted for the dipole fit of Eq. (9) we can define the corresponding mass determined from the z-expansion to be This relation will allow us to compare the radius extracted from the dipole and the z-expansion.

III. LATTICE METHODOLOGY
This section explains the methodology we use within lattice QCD in order to compute correlation functions and ensure ground state dominance. It also provides details on the gauge ensemble used in the analysis.

A. Correlation functions
For the computation of the correlation functions we use the standard nucleon interpolating field where C = γ 0 γ 2 is the charge conjugation matrix and u(x), d(x) the up and down quark fields. The two-point function in momentum space is then expressed as where with x 0 we denote 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. Γ 0 is the unpolarized positive parity projector Γ 0 = 1 2 (1 + γ 0 ). For the construction of the three-point correlation function the axial-vector current is inserted at a time slice, t ins , between the time of the creation and annihilation of states. The three-point function is given by with Γ ρ is the polarized projector, Γ ρ = iΓ 0 γ 5 γ ρ .

B. Ground state dominance
The interpolating field of Eq. (14) creates the nucleon ground state but also excited states. We apply Gaussian smearing [24,25] to the quark fields entering the interpolating field in order to increase the overlap with the ground state. See Ref. [14] for more details about our smearing procedure. To isolate the matrix element of interest we construct a ratio of three-to a combination of two-point functions [26][27][28][29] Overlap terms and time decaying exponentials cancel in the ratio. In Eq. (17) and from now on, we consider that t s and t ins are expressed relative to the source t 0 i.e. t s − t 0 → t s and t ins − t 0 → t ins . The ratio of Eq. (17) leads to the nucleon matrix element in the large time t s and t ins limits, that is where ∆E is the energy gap between the first excited state and the nucleon state. The rate of convergence to the nucleon state depends, besides the smearing procedure, also on the type of the insertion operator. In order to ensure ground state dominance, we employ three methods, namely a one state fit (plateau method), a two-state fit and the summation method. For a more detailed description about those three methods we refer the reader to Ref [14].
In this analysis we consider the same energy spectrum decomposition in both the two-and three-point functions. We determine the first excited state energy E 1 ( p) for each value of p by fitting the two-point function and use it when fitting the ratio of Eq. (17). We also fit the zero momentum two-point function to extract the nucleon mass and then use the continuum dispersion relation to determine the lowest state energy E 0 ( p) = m 2 N + p 2 for a given value of momentum. As shown in Ref. [14], the continuum dispersion relation is satisfied for all the momenta considered in this work.
From the nucleon matrix element one can determine the axial G A (Q 2 ) and induced pseudoscalar G P (Q 2 ) form factors using the decomposition of Eq. (6). Since there are several combinations of insertion, projector indices and momenta there is an over-constrained system of equations which determines the form factors. Details are given in Appendix B of Ref. [14].

C. Ensemble of gauge configurations
This work is based on the analysis of an N f = 2 + 1 + 1 twisted mass clover-improved fermion ensemble, referred to as cB211.072.64 (see Table I). In Ref. [14], where we studied the isovector form factors, we also analyzed two N f = 2 ensembles with the same light quark action, namely the cA2.09.48 and cA2.09.64 ensembles that have the same lattice spacing but different volumes (see Table I). Results on the axial form factors for the cA2.09.48 ensemble were first presented in Ref. [18]. For all the ensembles the lattice spacing is determined using the nucleon mass. More details are given in Refs. [30][31][32][33]. Finite volume effects have not been detected within the statistical precision obtained for the isovector quantities and therefore we do not consider them here.
These gauge configurations are produced by the Extended Twisted Mass Collaboration (ETMC) using the twisted mass fermion formulation [34,35] with a clover term [36] and the Iwasaki [37] improved gauge action. Since the simulation is done at maximal twist, we have automatic O(a) improvement for the physical observables studied here.  [30] but also the two N f = 2 ensembles cA2.09.48 [38] and cA2.09.64. cSW is the value of the clover coefficient and β = 6/g where g is the bare coupling constant. N f is the number of dynamical quark flavors, the lattice spacing is a and the lattice volume is V . mπ is the pion mass and mN the nucleon mass. L the spatial lattice length in physical units.

D. Disconnected three-point functions and statistics
The three-point function defined in Eq. (16), in general, has two different contributions: i) One in which the insertion operator couples directly to a valence quark in the nucleon, leading to the so-called connected three-point function, and ii) one in which the current couples to a sea quark giving the disconnected three-point function. In the case of the flavor isovector or octet currents, the disconnected contribution vanishes in the SU(3) flavor symmetric mass point and in the continuum limit. In the case of the flavor octet given in Eq. (3), and flavor singlet given in Eq. (2), disconnected contributions are non-zero. Since for the octet combination the disconnected contribution vanishes only in the SU(3) flavor symmetric limit, any non-vanishing contribution can be used to assess the level of SU(3) symmetric breaking. For the evaluation of the connected contributions we employ standard techniques, as discussed in Ref. [14], where we also give the statistics used for computing the connected contributions.
Here we describe our approach to compute the disconnected three-point functions. The disconnected quark loop for the axial-vector current is given by The trace of the all-to-all quark propagator D −1 (x ins ; x ins ), is the most computationally intensive quantity. Inverting from every point to every point on the lattice is computationally impossible for the lattice sizes considered in this work and the resources available. Instead, we combine stochastic methods to estimate the value of the quark loop. A novel method we employ is the combination of hierarchical probing [39] and deflation of low eigenvalues. Hierarchical probing allows for partitioning the lattice up to a distance in a hierarchical manner using the Hadamard vectors as a basis. The partitioning is done through a coloring approach, up to a distance 2 k where 2 d(k−1)+1 Hadamard vectors are needed (d = 4 for a 4-dimensional coloring). The computational cost of the method increases as 2 4 each time one increases the coloring distance. Thus, even if the probing is done in a hierarchical manner allowing to reuse the results from previous distances, the gain is small when we increase further the distance. Contributions from points beyond the probing distance are expected to be suppressed since the quark propagator decays exponentially fast with the distance from the diagonal. We further suppress such contributions using stochastic vectors that have the properties and where N r is the number of stochastic vectors. The off-diagonal contributions are suppressed by 1/ √ N r . The hierarchical probing method was first employed in studies for heavier than physical pion masses [19,40,41] yielding results with unprecedented accuracy. For simulations at the physical point, it was shown [31] that a larger probing distance is required, as expected, since the light quark propagator decays slower due to the smaller quark mass. Instead of increasing the probing distance, which translates to a significant increase in computational cost, we combine hierarchical probing with deflation of the low modes [42]. Namely, for the light quarks we construct the low mode contribution to the quark loops by computing exactly the smallest eigenvalues and corresponding eigenvectors of the squared Dirac operator and combine them with the contribution from the remaining higher modes, which is estimated using hierarchical probing. Additionally, we fully dilute in spin and color and employ the one-end trick [43], that was employed in our previous studies [18,44,45].
The parameters used for the evaluation of the quark loops are collected in Table II. Two hundred low modes of the square Dirac operator are computed in order to reduce the stochastic noise in the computation of the light quark loops.
For the charm quark we use a coloring distance 2 2 in hierarchical probing as compared to 2 3 used for the light and the strange quark loops. To increase the accuracy in the charm quark case we compute 12 stochastic vectors instead of one used for the light and strange quark loops. Nucleon two-point functions are evaluated for two hundred randomly chosen source positions that are sufficient for reducing the gauge noise for the large sink-source time separations of the disconnected three-point functions. Since they are available, we use the same number of two point functions for all sink-source time separations.  For disconnected quantities we are not limited to using p = 0 since no additional inversions are needed. Therefore, we consider several values of p , namely p = (2π/L) n up to n 2 = 2 and for p up to n 2 = 22. This allows us to compute the disconnected parts of the form factors for a higher density of Q 2 values.

IV. RENORMALIZATION FUNCTIONS
In order to relate the matrix elements computed on the lattice to physical observables one needs to renormalize. Here, we summarize our procedure. A more detailed description can be found in Ref. [46]. We employ a massindependent renormalization scheme and analyze five N f = 4 ensembles generated specifically for the determination of the renormalization functions. The value of β is the same as that of the cB211.072.64 ensemble. The pion masses are in the range of [366-519] MeV. These are used to take the chiral limit. The lattice volume is 24 3 × 48 for all N f = 4 gauge ensembles. The Rome-Southampton method, RI scheme [47], is employed where the quark propagators and vertex functions are non-perturbatively determined. For the axial-vector operator we need to renormalize with Z A and, since we consider also disconnected contributions, both singlet and non-singlet renormalization factors are needed.
We impose the following renormalization conditions: and S L (p) and Γ L (p) are the quark propagator and amputated vertex function, respectively, while S Born (p) and Γ Born (p) are the corresponding tree-level values. We note that the trace is meant to be taken over both spin and color indices and the RI renormalization scale is denoted by µ 0 . In order to compute the vertex functions non-perturbatively, we make use of momentum sources [48]. This allows to achieve per mil statistical accuracy on a very small sample of configurations [49,50]. With a high statistical precision, systematic errors need to also be under control. The momenta are chosen isotropic in the spatial direction, that is where n t ∈ [2,10], n x ∈ [2, 5] and T /a(L/a) are the temporal(spatial) lattice extent. The momenta satisfy the [51] in order to suppress the non-Lorentz invariant contributions. These appear in O(a 2 ) terms in the perturbative expansion of the Green's function and is expected to have non-negligible contributions from higher order in perturbation theory [46,49,50].
We improve the non-perturbative estimates by removing lattice artifacts in both Z q an Z A . The artifacts are calculated to one loop lattice perturbation theory [46]. In particular, one extracts the Greens functions of the axial operator using the same lattice action and values of the momentum p entering Eq. (23). For an optimal improvement, we calculate O(g 2 a ∞ ) terms, which cannot be obtained analytically. It should be noted that, the subtraction of the O(g 2 a ∞ ) terms can be done either at the level of the vertex functions Γ L (p) , or on Z A after the trace is taken. We have checked that both procedures lead to compatible results for the improved Z A . For consistency, we employ the subtraction in the final estimates of Z A , as performed in Ref. [14].
The evaluation of the Z-factors for the non-singlet current was presented in Ref. [14]. Here we present the evaluation of the singlet Z-factor, which is more complicated. For the computation of the singlet renormalization function Z s A we follow the same procedure as for the non-singlet case. In this case, in addition to the connected contributions, there are contributions from the disconnected quark loops. We employ the same noise reduction approaches discussed in Sec. III D for the evaluation of these disconnected contributions, namely we use hierarchical probing with 512 Hadamard vectors, the one-end trick and spin color dilution. Deflation is not used in this case since the N f ensembles are generated for heavy pion masses. In addition to the appearance of disconnected loops, a further complication is that, in contrast to the non-singlet case, Z s A is scheme and scale dependent. We express it in the MS-scheme, which is commonly used in experimental and phenomenological studies. The conversion procedure is applied on the Z-factors obtained on each initial RI scale (a µ 0 ), with a simultaneous evolution to a MS scale, chosen to be µ=2 GeV. In particular, we use the conversion factor calculated to 2-loops in perturbation theory [52].
In Fig. 1 we compare the non-singlet and singlet Z-factors. As can be seen, including the disconnected quark loop contributions lowers the value of the renormalization function and increases the error. We find Z A = 0.763(1) for non-singlet and Z s A = 0.753(5) for singlet.

V. ANALYSIS OF THE ISOSCALAR AXIAL FORM FACTORS
For the extraction of the axial and induced pseudoscalar form factors from the correlation functions we use the methodology presented in Sec. III. In order to identify the nucleon ground-state contribution we apply the three approaches discussed in Sec. III B. In Fig. 2, we demonstrate the effect of the excited-states contamination to the connected contribution for the isoscalar axial form factor G u+d A (Q 2 ) for two representative values of Q 2 . In the first column, we show the ratio of Eq. (17) for all the available values of t s . In the construction of the ratio, we use the two-point functions computed with the same source positions as the corresponding three-point functions to exploit their correlation leading to a reduction in the overall error. As t s increases, we observe a decrease in the values of the ratio. In the second column of Fig. 2 we show the values extracted by fitting the ratio to a constant excluding five time slices from the source and sink. This is done for t s /a > 12, yield a good χ 2 /d.o.f, namely in the range 0.7 to 1.1. In the third column of Fig. 2 we show results from the two-state and summation fits. For the two-state method we perform a simultaneous fit to all ratios for which t s ≥ t low s excluding t ins /a = 1, 2, t s − 1, t s − 2 and seek to identify convergence in the extracted value of the matrix element as we increase t low s . The resulting fit bands using these two-state fits are shown in the left panel. We show the prediction of the two-state fit of the time dependence of the ratio in the middle panel when we fix t ins = t s /2. As can be seen, the two-state fit prediction describes well the time dependence of the values extracted from the plateau method. We also observe that the value extracted using the two-state fit with t low s = 8a is consistent with the values extracted for higher values of t low s . The results from the summation method converge to those of the two-state fit for t low s > 1.2 fm i.e. at about half the t s value where the plateau fit yields convergent result. Based on these findings, we adopt as a criterion for the final value the one extracted from the two-state fit for the smallest t low s that shows convergence and is in agreement with the value from the summation method at some higher t low In Fig. 3, we present the excited-states contamination analysis for the case of the connected contributions to the isoscalar induced pseudoscalar form factor. In contrast to G A , suppression of excited states results in larger values for G P (Q 2 ) especially for the smaller Q 2 values. As Q 2 increases, contamination from excited states suppresses, with most of the plateau values being compatible with the two-state fit. We use the same criterion as for G u+d A for the selection of our final values. Therefore, we take the values extracted from the two-state fit for t low s = 8a.
In Fig. 4, we show the analysis to identify excited-state contributions for the disconnected parts contributing to G u+d A (Q 2 ) and G u+d P (Q 2 ). Although in these cases, all the sink-source time separations can be computed without additional cost, in practice, as the time separation t s increases, the errors become very large. Thus, we limit ourselves to t s ∈ [0.48 − 1.8] fm in what follows. As can be seen for both form factors, the disconnected contributions are nonzero. Eliminating excited states by increasing t s leads to more negative values for both axial and induced pseudoscalar form factors. In both cases, results for t s ≥ 16a extracted from the plateau method are in agreement with each other, as well as with those extracted using the two-state and summation methods. We thus opt to perform a weighted average of the converged plateau values to extract the final value.
The final values of the axial form factor G u+d A (Q 2 ) are shown in the left panel of Fig. 5, where we show separately the connected and disconnected contributions as a function of Q 2 . We observe that the connected contribution is positive, while the disconnected is negative. To extract the disconnected part, we combine various values of p ≥ 0 and thus can access a larger number of Q 2 values. We show also G u+d A (Q 2 ) after summing the connected and disconnected parts at the common Q 2 values. Since the disconnected contributions have a larger magnitude at smaller Q 2 values the slope of G u+d A (Q 2 ) at small Q 2 is smaller as compared to its connected part. We fit the Q 2 dependence as shown in the right panel of Fig. 5 using the dipole Ansatz and the z-expansion, as described in Sec. II, where for both fits the value at Q 2 = 0 is not a fit parameter but it is fixed form the forward matrix element yielding G u+d A (0) ≡ g u+d A . We find g u+d A = 0.436 (28) in agreement with our previous study [53]. The small difference is due to fact that in this work we use also p > 0 for the evaluation of the disconnected contributions, but also odd numbers of t s when averaging over the plateau values. Both fit forms describe the Q 2 behavior very well. The extracted values for the axial masses and the radii are given in Table III. The values extracted from the two fits are compatible, with the z-expansion yielding larger uncertainties. We note that by excluding larger values of Q 2 in the fit does not have an impact on the extracted parameters.  In Fig. 6 we show separately the connected and disconnected parts for the isoscalar induced pseudoscalar form factor G u+d P (Q 2 ). The disconnected part is of the same magnitude as the connected but with opposite sign. This has already been observed in previous studies [18,19]. This behavior leads to the cancellation of the sharp rise observed in the connected G u+d P (Q 2 ). Consequently, the isoscalar has an almost flat Q 2 -dependence within uncertainties, unlike the isovector combination where the pion pole gives a rapidly rising form factor at small Q 2 [14]. Also, the fact that the connected and disconnected parts are almost equal but with opposite sign means that G u+d P (Q 2 ) carries larger statistical errors.

VI. ANALYSIS OF THE STRANGE AND CHARM AXIAL FORM FACTORS
The strange and charm form factors receive only purely disconnected contributions. They probe sea quark degrees of freedom in the nucleon and provide us with an insight on their non-perturbative dynamics. Let us first examine how the ratio of Eq. (17) behaves when using the three-point function of the strange axial-vector current. In Fig. 7, we show the results on the ratio for different t s . As can be seen, although there is a trend to more negative values, the plateau region is consistent within the statistical uncertainties as we increase t s . This is also seen in the middle panel where we show the values extracted from plateau fits at various t s values. Furthermore, the summation and two-state fit methods yield results that are consistent with those extracted from the plateau fit for all t low s values. Given that the plateau values show convergence, we take the weighted average over the converged plateau values observed for t s 1.12 fm, resulting in the red band. The weighted average is also in agreement with the results from the two-state and summation fits, as we require to accept the final value.
The corresponding analysis of excited states for the three-point function of the charm axial-vector current is shown in Fig. 8. The three-point function in this case is more noisy and for clarity we only show the ratio for time separations up to 1 fm. As in the case of the strange three-point function, the plateau region of the ratio shows convergence as t s is increased within our current statistical accuracy. The results extracted using the summation method are noisy but yield consistent values. Two-state fits are omitted since, given the accuracy of the data, they are very noisy and thus yield no useful information. We take the weighted average of the converged plateau values to determine the final values on G c A (Q 2 ) and G c P (Q 2 ). The results for the strange axial form factor G s A (Q 2 ) are shown in left panel of Fig. 9. G s A (0) gives the strange axial charge and we find g s A = −0.044 (8) in agreement with the values reported in our previous analysis using the cB211.072.64 ensemble [53]. The small difference in the mean value is well within errors and is due to taking different data sets in the analysis. G s A (Q 2 ) is negative for all Q 2 values up to 1 GeV 2 . Both fits to a dipole form and the z-expansion describe the data well. The value at Q 2 = 0 is used as an input parameter. In Table IV, we give the χ 2 /d.o.f for the fits. The reason for the smaller χ 2 /d.o.f for the z-expansion is that higher order terms are taken into account that are sensitive to the values at larger Q 2 values giving rise to more curvature and thus a somewhat better description of the data.
The extracted values for the strange axial mass m s A and r.m.s. radius (r s A ) 2 are given in Table IV. Although the z-expansion fit has a steeper slope as Q 2 → 0 as compared to the dipole fit (see left panel of Fig. 9) the resulting values of the radius are consistent within the uncertainties. The strange induced pseudoscalar form factor G s P (Q 2 ) as a function of Q 2 is shown in the right panel of Fig. 9. As in the case of G s A (Q 2 ), G s P (Q 2 ) is clearly negative and large in magnitude especially at low Q 2 . The dipole and the z-expansion fits describe well the data. However, when we limit the fit range up to Q 2 = 0.5 GeV 2 the r.m.s and the value of the form factor at Q 2 = 0 are significantly larger. This is due to the curvature observed for small Q 2 .
We follow the same analysis described for the strange form factors to extract the charm axial form factors G c A (Q 2 ) and G c P (Q 2 ) that are shown in Fig. 10. and G c P (bottom) for Q 2 = 0.057 GeV 2 extracted using the plateau and the summation methods. The notation is the same as in Fig. 4. we can determine the same parameters as in the case of the strange form factors. The values are given in Table V. Since the slope of the z-expansion fit as Q 2 → 0 is steeper, the r.m.s. radius determined from the z-expansion tends to be larger as it was the case for the corresponding strange r.m.s. radius. It is worth mentioning that the z-expansion describes better the data as compared to the dipole Ansatz as indicated by the χ 2 /d.o.f. For the charm axial charge we find g c A = −0.0098 (17).  ) and G s P (Q 2 ) using the dipole Ansatz and the z-expansion. The notation is the same as that in Table III up to column five. The next columns are G s p (0), the value of the induced pseudoscalar form factor for Q 2 = 0, m s P the dipole mass and (r s P ) 2 the r.m.s radius.  , as a function of Q 2 . The notation is the same as that in Fig. 9.

VII. ANALYSIS OF THE FLAVOR SINGLET AND OCTET AXIAL FORM FACTORS AND THE SU(3) SYMMETRY BREAKING
The determination of isoscalar and strange form factors allows us to construct the corresponding SU(3) flavor octet and singlet form factors. We would like to highlight that these quantities are computed for the first time directly at the physical point.
In Fig. 11 we present results for the SU(3) flavor octet axial form factor G u+d−2s A (Q 2 ) and for the singlet G u+d+s A (Q 2 ).  If SU(3) was exact, the disconnected contributions would cancel in the octet combination. In practice, we find deviations from SU(3) symmetry especially at low Q 2 where the form factor is larger (see Fig. 12). This demonstrates that SU(3) flavor symmetry is violated due to the different mass between light and strange quarks. This is an important result since many phenomenological analyses assume SU(3) flavor symmetry introducing an uncontrolled systematic error. We find that there is up to 10% breaking for the axial and up to 50% for the induced pseudoscalar form factors. Due to the suppression of disconnected contributions in the octet combination, as can be seen in Fig. 12, G u+d−2s A (Q 2 ) is more precise as compared to G u+d+s A (Q 2 ) shown in Fig. 11. The data for both octet and singlet form factors are well described by our two fit Anzätze, namely the dipole form and the z-expansion. The resulting values of χ 2 /d.o.f are given in Table VI. The value of the form factors at zero momentum transfer, gives the octet and singlet axial charges g u+d−2s A and g u+d+s A , respectively. We find g u+d−2s A = 0.530 (22) and g u+d+s A = 0.384 (33). These charges have been also extracted from phenomenological analyses. In Ref. [54], the authors use polarised deep inelastic scattering data to extract g u+d−2s A = 0.46(5) and g u+d+s A = 0.36(3)(5) both in agreement with our findings but with larger uncertainties. It is worth mentioning that the analysis of Ref. [54] assumes SU(3) flavor symmetry.
In Table VI, we collect the parameters extracted from these fits. The SU(3) flavor octet axial mass m u+d−2s A tends to have a smaller value than the corresponding singlet, m u+d+s A , which translates to a bigger octet r.m.s. radius. However, statistical errors on the singlet quantities are large and the two values agree within the statistical errors. This is particularly true for the parameters extracted from the z-expansion where the statistical errors are even larger. The axial mass and radius determined from fitting the SU(3) flavor octet and singlet axial form factor G u+d−2s A (Q 2 ) and G u+d+s A (Q 2 ), respectively, using the dipole Ansatz and the z-expansion. The notation is the same as the one in Table III The Q 2 -dependence of the induced octet pseudoscalar form factors G u+d−2s P (Q 2 ) is shown in Fig. 13 with the corresponding extracted parameters provided in Table VII. It is well-known that the isovector induced pseudoscalar form factor, G u−d P (Q 2 ) has a pion pole behavior. Results on this form factor using the same ensemble were reported in Ref. [14]. By similar arguments, the SU(3) flavor octet form factor G u+d−2s P (Q 2 ) is expected to have an η pole behaviour. Since the η-meson has a much larger mass compared to the mass of the pion, the relations that hold in the chiral limit for G u−d P (Q 2 ) are expected to be significantly violated in this case.  In Fig. 13 we show also results on (m 2 η + Q 2 )G u+d−2s P (Q 2 ) that cancel the η-meson pole, as well as the dipole and z-expansion fits. This allows the extraction of the eta-nucleon coupling g ηN N in analogy to the determination of g πN N since where F 8 η is the decay constant of the η meson and m η its mass. We note that the mixing with the η has been neglected in Eq. (25). The η decay constant F 8 η can be determined directly in lattice QCD in an analogous manner to the computation of F π [55]. This will be computed for the current ensemble in a future work. Here, we use the value of F 8 η determined from phenomenology [56] to extract the coupling constant The extrapolation to Q 2 = −m 2 η is shown in the right panel of Fig. 13. As can be seen, the fact that one needs to perform a large extrapolation in the negative Q 2 region increases the statistical uncertainty as compared to the isovector case. The values extracted are in agreement with the ones extracted from phenomenological studies [56][57][58].
If one defines a Goldberger-Treiman discrepancy for the octet in a similar manner as done for the isovector combination can assess how much the Goldberger-Treiman relation is violated in this case. We find that We find a violation of about 40-50% for the octet combination of ∆ GT , which is much larger than the 2% determined for the isovector combination [14]. This is a consequence of the large η-meson mass. The flavor singlet induced pseudoscalar form factor G u+d+s P (Q 2 ) is noisy because of the disconnected contributions are large and of opposite sign to the connected partly canceling each other, as can be seen in Fig.14  FIG. 14: Results on the flavor singlet induced pseudoscalar form factor. The notation is the same as that in Fig. 6.

VIII. THE UP AND DOWN AXIAL AND INDUCED PSEUDOSCALAR FORM FACTORS
Having determined the isovector [14] and isoscalar form factors we can disentangle the up and down quark contributions to the these form factors.
In Fig. 15 we show results for the up and down quark axial form factors, G u A (Q 2 ) and G d A (Q 2 ) as a function of quark charges obtained at Q 2 = 0 are g u A = 0.859 (18) and g d A = −0.423 (17) in agreement with the values found in Ref. [53]. Since the value of the form factors at zero momentum is known, we use it to eliminate one fit parameter in the jackknife analysis. The values of the up and down quark axial masses and r.m.s. radii extracted from the dipole fit and using the z-expansion are given in Table VIII. We find that m u A ∼ m d A and (r u A ) 2 ∼ (r d A ) 2 within statistical errors. The axial mass and radius determined from fitting G u A (Q 2 ) and G d A (Q 2 ), using the dipole Ansatz and the z-expansion. The notation is the same as that in Table VI.  form and the z-expansion we eliminate the pion pole and consider instead for the fits. Note that GP (Q 2 ) has units of GeV 2 . Like in the case of G u,d A (Q 2 ), G u P (Q 2 ) is positive and G d P (Q 2 ) is negative. However, unlike G u,d A (Q 2 ), both G u P (Q 2 ) and G d P (Q 2 ) have similar magnitude. For the case of GP (Q 2 ), the Q 2 -dependence is different being more linear for the up quark as compared to down quark. The dipole fits to both G ũ P (Q 2 and G d P (Q 2 ) do not describe the curvature as well as the z-expansion fit does, producing more curvature for the former and less for the latter as compared to the lattice QCD data.
In Table IX, we provide the parameters extracted from the up and down induced pseudoscalar form factors. To relate the parameters, we utilize the relations , and r 2 P = 6 m 2 π + r 2 P .

IX. FINAL RESULTS
In this section we collect our final results extracted from the fits to the axial and induced pseudoscalar form factors for the various flavor combinations. Results are provided using the z-expansion given in Eq. (10), since in most cases it fits better the form factors as in, e.g., G u,d P (Q 2 ). Fits to the dipole Ansatz are used as a determination of the systematic error due to the choice of the fit form, by taking the difference between the z-expansion and dipole fit values. In addition, we use the two different Q 2 fit ranges, namely Q 2 0.5 and Q 2 1 GeV 2 to extract a systematic due to the fit range dependence. We quote as the parameters extracted using as upper range Q 2 1 GeV 2 in the fit and the difference between the mean values extracted using the two ranges as the systematic error. In the first column, we give the quark flavor combination considered, in the second and third columns the axial mass and r.m.s radii and in the rest three columns the value of the form factor extrapolated to Q 2 = 0, the dipole mass and r.m.s radii extracted from fitting the induced pseudoscalar form factor. The first error is purely statistical, the second is a systematic due to different fit ranges and the third is a systematic due to the two different forms used to fit the Q 2 -dependence.
Comb In Table X, results for the axial and induced pseudoscalar masses and r.m.s radii for the quark flavor combinations considered are provided. It is worth mentioning that this is the first time that these radii are determined for each quark flavor separately but also for the octet and singlet combinations providing us with detailed information on the structure properties of the nucleon. A notable finding is that the axial strange and charm r.m.s radii tend to be larger than those for the light quarks. However, the uncertainties are still large and we would need to improve the accuracy in order to draw any definite conclusion. The values of G u,d P (0) for both up and down quarks are very large compared to the rest due to the presence of the pion pole. This sharp rise of these two form factors is reflected in the extracted r.m.s radii which are significantly larger than all the rest.
From the flavor octet combination we can determine the pseudoscalar η-meson-nucleon coupling. The value is given in Eq. (33), where the first error is statistical and the second is a systematic due to the fit form used. This is the first determination at the physical point. It is, however, in agreement with a previous lattice QCD study [19] for an ensemble with pion mass of 317 MeV. The Goldberger-Treiman discrepancy ∆ GT is also determined for the octet combination and it is found to be 50% which highlights that such relations are badly broken for mesons with much larger mass than the pion. , and the N f = 2 ensemble cA2.09.48 (filled blue squares) presented in Ref. [18]. See Table I for details on the parameters of the two ensembles. The upper plots show results for the isoscalar and strange axial form factor while the lower plot for the charm axial form factor.
The form factors presented in this work were studied previously by only another lattice QCD group, namely the LHPC [19] but using an ensemble with pion mass m π = 317 MeV. Here we restrict the comparison to studies performed directly at the physical point and, therefore the only other available results are provided from our previous work [18] using the N f = 2 ensemble cA2.09.48 of Table I. In that study we didn't employ the improved noise reduction approaches for the evaluation of the quark loops that we use in this work and presented in Sec. III D. In our previous study we used volume sources without spin nor color dilution. For the light quarks we used 2250 stochastic sources while for the strange we used 1024 and for the charm 1250. For the strange and charm quark loops we also used the truncated solver method [59]. We note that hierarchical probing used in this work was not used for the analysis of the cA2.09.48. Results for the axial form factors are compared in Fig. 18. For the isoscalar combination only few Q 2 values are available in the case of the cA2.09.48 ensemble, namely up to Q 2 = 0.3 GeV 2 . For the strange axial form factor, the results using the cA2.09.48 ensemble are very noisy. This comparison provides a nice demonstration of the improvements accomplished in this work with about only twice the computational effort. The situation is similar for the case of the charm axial form factor.
In Fig. 19 we compare the results for the induced pseudoscalar form factor. We observe agreement between the results using the two ensembles with the results of the current work being significantly more precise.

XI. CONCLUSIONS
The complete flavor decomposition of the axial and induced pseudoscalar form factors of the nucleon is determined directly at the physical point using one N f = 2 + 1 + 1 ensemble of twisted mass fermions. We obtain non-zero results for the up, down, strange and charm quark form factors to increased accuracy as compared to our previous study using an N f = 2 twisted mass ensemble [18]. This is accomplished by using a combination of deflation of lower mode, hierarchical probing, spin-colour dilution and the one-end trick.
These results provide valuable input to on-going and planned parity-violating experiments. They are also crucial for the cross sections for a class of popular cold dark matter candidates [1]. Having the complete flavor decomposition allows us to check for SU(3) flavor symmetry. We find that SU(3) symmetry is broken up to about 10% for the octet axial and up to 50% for the induced pseudoscalar form factors with the breaking being larger at low Q 2 values. This is an important result since many phenomenological studies assume SU(3) flavor symmetry, and thus carry an uncontrolled systematic error.
In the future we plan to analyze two additional N f = 2 + 1 + 1 twisted mass fermion ensembles with smaller lattice spacings so that we can take the continuum limit. This will also enable us to check the PCAC relation directly in the continuum limit eliminating any cut-off effects that may cause violations.