Nucleon axial form factors using lattice QCD simulations with a physical value of the pion mass

TL;DR

This paper reports on lattice QCD simulations at the physical pion mass to compute nucleon axial form factors, including isovector, isoscalar, strange, and charm contributions, highlighting significant disconnected quark loop effects.

Contribution

It introduces a novel lattice QCD approach with deflation techniques to accurately evaluate disconnected contributions to nucleon axial form factors at the physical pion mass.

Findings

Disconnected quark loops are non-zero and large for the induced pseudo-scalar form factor.

The study provides detailed isovector, isoscalar, strange, and charm axial form factors.

Results improve understanding of nucleon structure and quark flavor contributions.

Abstract

We present results on the nucleon axial and induced pseudo-scalar form factors using an ensemble of two degenerate twisted mass clover-improved fermions generated at the physical value of the pion mass. We evaluate the isovector and the isoscalar, as well as, the strange and the charm axial form factors. The disconnected contributions are evaluated using recently developed methods that include deflation of the lower eigenstates, allowing us to extract the isoscalar, strange and charm axial form factors. We find that the disconnected quark loop contributions are non-zero and particularly large for the induced pseudo-scalar form factor.

Tables

Table 1. Table 1: Simulation parameters of the ensemble used here. The nucleon and pion mass and the lattice spacing have been determined in Ref.  Alexandrou and Kallidonis ( 2017, 1704.02647 ) .

$\beta$=2.1, $c_{\mathrm{SW}}$=1.57751, $a$=0.0938(3) fm, $r_0/a$=5.32(5)
48${}_{}{}^3\times$96, $L$=4.5 fm	$a{\mu }_l$=	0.0009
	$m_{\pi }$=	0.1304(4) GeV
	$m_{\pi }L$=	2.98(1)
	$m_N$=	0.932(4) GeV
	$m_N/m_{\pi }$=	7.15(4)

Table 2. Table 2: The statistics of our calculation. N conf subscript 𝑁 conf N_{\rm conf} is the number of gauge configurations analyzed and N src subscript 𝑁 src N_{\rm src} is the number of source positions per configuration. For the disconnected contributions, N r HP superscript subscript 𝑁 𝑟 HP N_{r}^{\rm HP} is the number of high-precision stochastic vectors produced, and N r LP superscript subscript 𝑁 𝑟 LP N_{r}^{\rm LP} is the number of low-precision vectors used when employing the TSM.

Connected			Disconnected
$t_s/a$	$N_{\mathrm{conf}}$	$N_{\mathrm{src}}$	Flavor	$N_{\mathrm{conf}}$	$N_r^{\mathrm{HP}}$	$N_r^{\mathrm{LP}}$	$N_{\mathrm{src}}$
10	579	16	light	2120	2250	-	100
12	579	16	strange	2057	63	1024	100
14	579	16	charm	2034	5	1250	100

Table 3. Table 3: We give the values extracted from the plateau method for the isovector and isoscalar axial charges and for the axial charge of the individual quarks. The first error is the statistical error determined using jackknife and the second is the systematic error due to the excited states computed as the difference in the mean value between the plateau fit and the two-state fit. The experimental values have been taken from Ref.  Patrignani et al. ( 2016 ) .

	This work	Experiment
$g_A$	1.212(33)(22)	1.2723(23)
$g_A^{u+d}$ (Conn.)	0.595(28)(1)	-
$g_A^{u+d}$ (Disc.)	-0.150(20)(19)	-
$g_A^{u+d}$	0.445(34)(19)	0.416(18)
$g_A^u$	0.827(30)(5)	0.843(12)
$g_A^d$	-0.380(15)(23)	-0.427(12)
$g_A^s$	-0.0427(100)(93)	-
$g_A^c$	-0.00338(188)(667)	-

Table 4. Table 4: Extracted values for the axial masses and corresponding axial radii using dipole fits to Eq. ( 32 ) with their associated χ 2 / d.o.f superscript 𝜒 2 d.o.f \chi^{2}/\textrm{d.o.f} . The central value and statistical error is from fits to results using the plateau method at t s = 1.31 subscript 𝑡 𝑠 1.31 t_{s}=1.31  fm. The first error is statistical while the second is the systematic due to excited states, taken as the difference between the central value and the value extracted from the two-state fit.

Form factor	$m_A$ [GeV]	$⟨r_A^2⟩$ $[$fm${}_{}{}^2]$	${\chi }^2$/d.o.f
$G_A^{u-d}$	1.322(42)(17)	0.266(17)(7)	0.35
$G_A^{u+d}$	1.736(244)(374)	0.155(43)(96)	0.64
$G_A^u$	1.439(28)(114)	0.225(28)(40)	0.61
$G_A^d$	1.243(49)(133)	0.301(24)(55)	0.42
$G_A^s$	0.921(228)(90)	0.549(272)(93)	0.30

Table 5. Table 5: Our values for the axial form factor for various values of Q 2 superscript 𝑄 2 Q^{2} . The first column gives Q 2 superscript 𝑄 2 Q^{2} in GeV 2 superscript GeV 2 \mathrm{GeV^{2}} , the second is the isovector axial form factor, the third is the connected contribution to the isoscalar axial form factor, the forth column gives the total value of the isoscalar, while the fifth and the sixth columns are the total contributions from the up and down quarks separately.

$Q^2\left[Ge{V}^2\right]$	$G_A^{u-d}$	$G_A^{u+d}$(Conn.)	$G_A^{u+d}$(Tot.)	$G_A^u$(Tot.)	$G_A^d$(Tot.)
0.0000	1.212(33)(22)	0.595(29)(1)	0.445(35)(18)	0.827(30)(5)	-0.380(15)(23)
0.0753	1.110(20)(16)	0.551(17)(11)	0.439(23)(20)	0.772(18)(2)	-0.339(12)(16)
0.1477	1.035(14)(17)	0.518(14)(10)	0.430(18)(44)	0.728(13)(12)	-0.308(10)(29)
0.2174	0.970(18)(10)	0.484(18)(13)	0.389(27)(33)	0.676(19)(8)	-0.294(13)(23)
0.2849	0.911(20)(5)	0.458(17)(1)	0.377(31)(45)	0.641(20)(29)	-0.270(17)(29)
0.3502	0.855(18)(9)	0.438(14)(1)	-	-	-
0.4135	0.802(20)(9)	0.413(14)(3)	-	-	-
0.5351	0.701(25)(21)	0.385(17)(23)	-	-	-
0.5936	0.678(23)(18)	0.364(15)(1)	-	-	-
0.6506	0.636(28)(53)	0.321(18)(17)	-	-	-
0.7064	0.588(28)(65)	0.289(24)(35)	-	-	-
0.7609	0.573(50)(45)	0.307(42)(29)	-	-	-
0.8143	0.520(30)(73)	0.265(23)(7)	-	-	-
0.8666	0.533(37)(52)	0.301(24)(22)	-	-	-
0.9683	0.096(1.665)(161)	0.056(470)(169)	-	-	-

Equations

$$⟨N(p^{\prime },s^{\prime })|A_{\mu }|N(p,s)⟩$$

$$⟨N(p^{\prime },s^{\prime })|A_{\mu }|N(p,s)⟩=i\sqrt{\frac{m_N^2}{E_N({\vec{p}}^{\prime })E_N(\vec{p})}}{\bar{u}}_N(p^{\prime },s^{\prime })\left({\gamma }_{\mu }{G}_A(Q^2)-i\frac{Q_{\mu }}{2m_N}{G}_p(Q^2)\right){\gamma }_5{u}_N(p,s),$$

$$A_{\mu }^{\mathrm{isov}}=\bar{\psi }(x){\gamma }_{\mu }{\gamma }_5{\tau }^3\psi (x),A_{\mu }^{\mathrm{isos}}=\bar{\psi }(x){\gamma }_{\mu }{\gamma }_{5}1\psi (x),A_{\mu }^s=\bar{s}(x){\gamma }_{\mu }{\gamma }_{5}s(x)\text{and}{A}_{\mu }^c=\bar{c}(x){\gamma }_{\mu }{\gamma }_{5}c(x).$$

$$S_F\left[\chi ,\overline{\chi },U\right]=a^4\sum _x\overline{\chi }(x)\left(D_W[U]+i{\mu }_l{\gamma }_5{\tau }^3-\frac{1}{4}{c}_{\mathrm{SW}}{\sigma }^{\mu \nu }{\mathcal{F}}^{\mu \nu }[U]\right)\chi (x),$$

$${\mathcal{F}}^{\mu \nu }[U]=\frac{1}{8}\left[P_{\mu ,\nu }(x)+P_{\nu ,-\mu }(x)+P_{-\mu ,-\nu }(x)+P_{-\nu ,\mu }(x)-(h.c.)\right],$$

$$\psi (x)=\frac{1}{\sqrt{2}}\left(11+i{\tau }^3{\gamma }_5\right)\chi (x),\overline{\psi }(x)=\overline{\chi }(x)\frac{1}{\sqrt{2}}\left(11+i{\tau }^3{\gamma }_5\right).$$

$$J(x)={ϵ}^{abc}\left(u^{a⊺}(x)\mathcal{C}{\gamma }_5{d}^b(x)\right)u^c(x),$$

$$G_{\mu }({\Gamma }_{\nu },\vec{q},\vec{p}{}^{\prime };t_s,t_{\mathrm{ins}},t_0)=\sum _{{\vec{x}}_{\mathrm{ins}},{\vec{x}}_s}{e}^{i({\vec{x}}_{\mathrm{ins}}-{\vec{x}}_0)\cdot \vec{q}}\mathrm{Tr}\left[{\Gamma }_{\nu }⟨J(t_s,{\vec{x}}_s)A_{\mu }(t_{\mathrm{ins}},{\vec{x}}_{\mathrm{ins}})\bar{J}(t_0,{\vec{x}}_0)⟩\right]e^{-i({\vec{x}}_s-{\vec{x}}_0)\cdot \vec{p}{}^{\prime }},$$

$$C({\Gamma }_0,\vec{p};t_s,t_0)=\sum _{{\vec{x}}_s}{e}^{-i({\vec{x}}_s-{\vec{x}}_0)\cdot \vec{p}}\mathrm{Tr}\left[{\Gamma }_0⟨J(t_s,{\vec{x}}_s)\bar{J}(t_0,{\vec{x}}_0)⟩\right],$$

$$R_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p};t_s,t_{\mathrm{ins}})=\frac{G_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p};t_s,t_{\mathrm{ins}})}{C({\Gamma }_0,\vec{p}{}^{\prime };t_s)}\sqrt{\frac{C({\Gamma }_0,\vec{p};t_s-t_{\mathrm{ins}})C({\Gamma }_0,\vec{p}{}^{\prime };t_{\mathrm{ins}})C({\Gamma }_0,\vec{p}{}^{\prime };t_s)}{C({\Gamma }_0,\vec{p}{}^{\prime };t_s-t_{\mathrm{ins}})C({\Gamma }_0,\vec{p};t_{\mathrm{ins}})C({\Gamma }_0,\vec{p};t_s)}}$$

$$R_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p};t_s,t_{\mathrm{ins}})\underset{t_s-t_{\mathrm{ins}}\gg 1}{\overset{t_{\mathrm{ins}}\gg 1}{\to }}{\Pi }_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p}).$$

$$G_{\mu }(\vec{p}{}^{\prime },\vec{p},t_s,t_{\mathrm{ins}})$$

$$C(\vec{p},t_s)=c_0(\vec{p})e^{-E_0(\vec{p})t_s}+c_1(\vec{p})e^{-E_1(\vec{p})t_s}.$$

$$\mathcal{M}=\frac{A_{00}(\vec{p}{}^{\prime },\vec{p})}{\sqrt{c_0(\vec{p}{}^{\prime })c_0(\vec{p})}}.$$

$$R_{\mu }^{\mathrm{sum}}({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p},t_s)\equiv \sum _{t_{\mathrm{ins}}=a}^{t_s-a}{R}_{\mu }({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p},t_s,t_{\mathrm{ins}})=C+t_s\mathcal{M}+\mathcal{O}(e^{-\Delta t_s})+\cdots ,$$

$$d{R}_{\mu }^{\mathrm{sum}}({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p},t_s)=\frac{R_{\mu }^{\mathrm{sum}}({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p},t_s+d{t}_s)-R_{\mu }^{\mathrm{sum}}({\Gamma }_{\nu },\vec{p}{}^{\prime },\vec{p},t_s)}{d{t}_s}$$

$$L^{(f)}(t_{\mathrm{ins}},\vec{q})=\sum _{{\vec{x}}_{\mathrm{ins}}}\mathrm{Tr}\left[M_f^{-1}(x_{\mathrm{ins}};x_{\mathrm{ins}})\mathcal{G}\right]e^{+i\vec{q}\cdot {\vec{x}}_{\mathrm{ins}}}$$

$$\sum _{{\vec{x}}_{\mathrm{ins}}}\mathrm{Tr}\left[i{\gamma }_5{M}_{{\chi }_u}^{-1}(x_{\mathrm{ins}};x_{\mathrm{ins}})-i{\gamma }_5{M}_{{\chi }_d}^{-1}(x_{\mathrm{ins}};x_{\mathrm{ins}})\right]$$

$$L^{u+d}(t_{\mathrm{ins}};\vec{q})$$

$$\frac{1}{N_r}\sum _r|{\xi }_r⟩⟨{\xi }_r|=1+\mathcal{O}\left(\frac{1}{\sqrt{N_r}}\right),$$

$$L^{u+d}(t_{\mathrm{ins}};\vec{q})$$

$${\sigma }^2\propto 2(1-r_c)+\frac{N_r^{\mathrm{HP}}}{N_r^{\mathrm{LP}}},$$

$$\frac{N_r^{\mathrm{LP}}}{N_r^{\mathrm{HP}}}=\sqrt{\frac{n_{\mathrm{HP}}}{2(1-r_c)n_{\mathrm{LP}}}},$$

$$G_{\mathcal{G}}^{ns}(p)=\frac{a^{12}}{V}\sum _{x,y,z}⟨u(x)\bar{u}(z)\mathcal{G}d(z)\bar{d}(y)⟩e^{-ip(x-y)},G_{\mathcal{G}}^s(p)=\frac{a^{12}}{V}\sum _{x,y,z}⟨u(x)\bar{u}(z)\mathcal{G}u(z)\bar{u}(y)⟩e^{-ip(x-y)}$$

$${\Lambda }_{\mathcal{G}}(p)=(S(p){)}^{-1}{G}_{\mathcal{G}}(p)(S(p){)}^{-1}$$

$$Z_{q}^{-1}Z_{\mathcal{G}}\operatorname{Tr}\left[\Lambda_{\mathcal{G}}(p)\Lambda_{\mathcal{G}}^{\text{tree}}\right]=\operatorname{Tr}\left[\Lambda_{\mathcal{G}}^{\text{tree}}\;\Lambda_{\mathcal{G}}^{\text{tree}}\right]\;\;\;\text{with}\;\;\;Z_{q}=-\frac{i}{4}\operatorname{Tr}\left[\frac{\frac{1}{a}\gamma_{\rho}\sin(ap_{\rho})}{\frac{1}{a^{2}}\sum_{\rho}\sin^{2}(ap_{\rho})}S^{-1}(p)\right]\Bigg{|}_{p_{\rho}=\mu_{\rho}}.$$

$${\Pi }_k({\Gamma }_k,\overrightarrow{p^{\prime }},\vec{p})$$

$${\Pi }_i({\Gamma }_k,\overrightarrow{p^{\prime }},\vec{p})$$

$$C=\frac{2m_N^2}{E(\vec{p})(E(\overrightarrow{p^{\prime }})+m_N)}\times \sqrt{\frac{E(\vec{p})(E(\overrightarrow{p^{\prime }})+m_N)}{E(\overrightarrow{p^{\prime }})(E(\vec{p})+m_N)}}.$$

$${\Pi }_i(k,\vec{p}{}^{\prime },\vec{p})=D_i(k,\vec{p}{}^{\prime },\vec{p})F(Q^2),$$

Full text

DESY 17-064

**Nucleon axial form factors using $N_{f}=2$ twisted mass fermions with a physical value of the pion mass **

C. Alexandrou(a,b), M. Constantinou*(c), K. Hadjiyiannakou(b), K. Jansen(d), C. Kallidonis(b), G. Koutsou(b), A. Vaquero Aviles-Casco(e)* **

(a) Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus

(b) Computation-based Science and Technology Research Center, The Cyprus Institute, 20 Kavafi Str., Nicosia 2121, Cyprus

(c) Temple University,1925 N. 12th Street, Philadelphia, PA 19122-1801, USA

(d) NIC, DESY, Platanenallee 6, D-15738 Zeuthen, Germany

(e) Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA

We present results on the nucleon axial and induced pseudo-scalar form factors using an ensemble of two degenerate twisted mass clover-improved fermions with mass yielding a pion mass of $m_{\pi}=130$ MeV. We evaluate the isovector and the isoscalar, as well as, the strange and the charm axial form factors. The disconnected contributions are evaluated using recently developed methods that include deflation of the lower eigenstates, allowing us to extract the isoscalar, strange and charm axial form factors. We find that the disconnected quark loop contributions are non-zero and particularly large for the induced pseudo-scalar form factor.

Nucleon axial form factors, Nucleon Structure, Lattice QCD

I Introduction

Understanding the structure of the nucleon from first principles constitutes one of the key endeavors of both nuclear and particle physics. Despite the long history of experimental activity its structure is not yet fully understood. This includes the portion of its spin carried by quarks as well as the charge radius of the proton. While electromagnetic form factors have been well studied experimentally, the axial form factors are known to less accuracy. An exception is the nucleon axial charge, which has been measured from $\beta$-decays to very high precision. Two methods have been extensively used to determine the momentum dependence of the nucleon axial form factor. The most direct method is using elastic scattering of neutrinos and protons, typically $\nu_{\mu}\;+p\;\longrightarrow\;\mu^{+}\;+\;n$  Ahrens et al. (1988). The second method is based on the analysis of charged pion electro-production data Bernard et al. (1992) off the proton, which is slightly above the pion production threshold. The induced pseudo-scalar form factor $G_{p}(q^{2})$ is even harder to measure experimentally. For the case of the induced pseudo-scalar coupling $g_{p}$, a range of muon capture experiments, as proposed in Ref. Bardin et al. (1981), have been carried out for its determination (see Ref. Gorringe and Fearing (2004) for a review). The form factor $G_{p}(q^{2})$ is less well known, and has only been determined at three values of the momentum transfer from the longitudinal cross section in pion electro-production Choi et al. (1993).

Lattice QCD presents a rigorous framework for computing the axial form factors from first principles, in particular in light of the tremendous progress made in simulating the theory at near physical values of the quark masses, large enough volumes and small enough lattice spacings. Having simulations using the physical values of the light quarks eliminates chiral extrapolations, which for the baryon sector introduced a large systematic uncertainty. In addition, improved algorithms and novel computer architectures have enabled the computation of contributions due to disconnected quark loops, which previously were mostly neglected.

In this work we present results for the nucleon axial and induced pseudo-scalar form factors from an ensemble generated with two degenerate quarks with masses fixed approximately to their physical value Abdel-Rehim et al. (2017). We study both the isovector and isoscalar combinations as well as the strange and charm form factors, which receive only disconnected contributions.

The paper is organized as follows: in Section II we introduce the axial form factors and the nucleon axial matrix element and in Section III we give details of the lattice action used. In Section IV we explain our set-up, the correlation functions used and the methods employed to extract the nucleon matrix elements from the lattice data. The renormalization process is described in Section V and in Section VI we present our results. Finally, in Section VII we conclude.

II Axial form factors

To extract the axial and pseudo-scalar form factors one needs to evaluate the nucleon matrix element

[TABLE]

where $A_{\mu}(x)=\bar{\psi}(x)\gamma_{\mu}\gamma_{5}\tau^{a}\psi(x)$ is the axial-vector current with $\psi(x)=(u(x),d(x))$ the doublet of up- and down-quarks, $\tau^{a}$ a Pauli matrix acting on flavor components and $p,s$ ($p^{\prime},s^{\prime}$) are the momentum and spin of the initial (final) nucleon state, $N$. The nucleon matrix element of the axial-vector current decomposes into two form factors, $G_{A}(q^{2})$ and $G_{p}(q^{2})$, which are functions of the momentum transfer squared $q^{2}=(p^{\prime}_{\mu}-p_{\mu})^{2}=-Q^{2}$. In a lattice QCD computation one performs a Wick rotation to imaginary time. Working in Euclidean space the nucleon matrix element of the axial-vector operator can be written in the continuum as

[TABLE]

where $u_{N}$ are nucleon spinors and $m_{N}$ and $E_{N}(\vec{p})$ are the nucleon mass and energy at momentum $\vec{p}$. In this work, we consider the isovector, isoscalar as well as strange and charm combinations

[TABLE]

In the isovector case disconnected contributions cancel in the isospin limit. For the isoscalar combination both connected and disconnected contributions enter, while for the strange and charm form factors we have only disconnected contributions. In this work the disconnected contributions are computed for the first time using simulations with a physical value of the pion mass. The connected and disconnected three-point functions are represented schematically in Fig. 1.

III Lattice action

In this work we use a single gauge ensemble of two degenerate ($N_{f}=2$) up and down twisted mass quarks with mass tuned to reproduce approximately the physical pion mass Abdel-Rehim et al. (2017). The parameters of our calculation are shown in Table 1. The “Iwasaki” improved gauge action is used Iwasaki (1983); Abdel-Rehim et al. (2014a) for the gluonic part. In the fermion sector, the twisted mass fermion action for a doublet of degenerate quark flavors Frezzotti et al. (2001); Frezzotti and Rossi (2004a) is employed including in addition a clover-term Sheikholeslami and Wohlert (1985).

The fermionic action is given by

[TABLE]

where $D_{W}$ is the Wilson-Dirac operator, $\mu_{l}$ is the bare twisted light quark mass and $-\frac{1}{4}c_{\rm SW}\sigma^{\mu\nu}\mathcal{F}^{\mu\nu}[U]$ is the clover-term, with $c_{\rm SW}$ the so-called Sheikoleslami-Wohlert improvement coefficient. The field strength tensor $\mathcal{F}^{\mu\nu}[U]$ is given by Sheikholeslami and Wohlert (1985)

[TABLE]

where $P_{\mu,\nu}(x)$ is a fundamental $1\times 1$ Wilson plaquette and $\sigma^{\mu\nu}=\frac{1}{2}[\gamma_{\mu},\gamma_{\nu}]$. We take $c_{\rm SW}=1.57551$ from Ref. Aoki et al. (2006). The quark fields denoted by $\chi$ in Eq. (4) are in the so-called “twisted basis”. The fields in the “physical basis” denoted by $\psi$, are obtained at maximal twist by the transformation

[TABLE]

In this paper, unless otherwise stated, the quark fields will be understood as “physical fields”, $\psi$, in particular when we define the interpolating fields.

Twisted mass fermions (TMF) provide an attractive formulation for lattice QCD allowing for automatic ${\cal O}(a)$ improvement, infrared regularization of small eigenvalues and fast dynamical simulations Frezzotti and Rossi (2004a). However, the $\mathcal{O}(a^{2})$ lattice artifacts that the twisted mass action exhibits lead to instabilities in the numerical simulations, particularly at lower values of the quark masses and influence the phase structure of the lattice theory Aoki (1984); Sharpe and Singleton (1998); Farchioni et al. (2005). The clover-term was added in the TMF action to allow for smaller $\mathcal{O}(a^{2})$ breaking effects between the neutral and charged pions that lead to the stabilization of simulations with light quark masses close to the physical pion mass retaining at the same time the particularly significant ${\cal O}(a)$ improvement that the TMF action features. The reader interested in more details regarding the twisted mass formulation is referred to Refs. Frezzotti et al. (2001); Frezzotti and Rossi (2004a, b); Frezzotti et al. (2006); Boucaud et al. (2008) and for the simulation strategy to Refs. Abdel-Rehim et al. (2017, 2015).

IV Lattice evaluation of the nucleon matrix elements

In order to extract the nucleon matrix elements, we need an appropriately defined three-point function and the nucleon two-point function. To construct these correlation functions one creates states with the quantum numbers of the nucleon from the vacuum at some initial time (source) and annihilates them at a later time (sink). The commonly used nucleon interpolating field is given by

[TABLE]

where $\mathcal{C}$ is the charge conjugation matrix. To improve the overlap of this operator with the ground state we employ Gaussian smearing Gusken (1990); Alexandrou et al. (1994) to the quark fields at the source and the sink. In addition, we apply APE smearing Albanese et al. (1987) to the gauge links entering the hopping matrix in order to reduce unphysical ultra-violet fluctuations.

The three-point function in momentum space can be written as

[TABLE]

and the two-point function is given by

[TABLE]

where $\vec{q}=\vec{p}\,^{\prime}-\vec{p}$ is the momentum transfer. For the two-point function we use the projector \Gamma_{0}=\frac{1}{4}(\hbox{1!!1}+\gamma_{0}), whereas for the three-point function the projector is $\Gamma_{j}=i\Gamma_{0}\gamma_{5}\gamma_{j}$, $j=1,2,3$, which permits the extraction of the axial $G_{A}(Q^{2})$ and the induced pseudo-scalar $G_{p}(Q^{2})$ form factors. The matrix element can be extracted by taking appropriate combinations of three- and two-point functions. An optimal ratio Hagler et al. (2003), which cancels unknown overlap terms and time dependent exponentials is,

[TABLE]

where we measure all times relative to the time of the source, i.e. $t_{\rm ins}$ and $t_{s}$ measure the time separation of the current insertion and the sink, respectively, from the source. The ratio becomes time independent in the large time limit yielding a plateau $\Pi_{\mu}$ from where the matrix element of the ground state is extracted, defined via:

[TABLE]

In practice, the source-insertion and insertion-sink time separations cannot be chosen arbitrarily large because the gauge noise becomes dominant, thus several time separations must be tested to ensure convergence to the ground state. It is expected that different matrix elements have different sensitivity to excited states. In the case of the scalar operator, it has been shown that source-sink separations larger than $t_{s}=1.5$ fm are required in order to damp out sufficiently excited state effects Abdel-Rehim et al. (2016). For the axial-vector current excited state contamination is found to be less severe at least for pion masses larger than physical used in previous calculations Alexandrou et al. (2011, 2007). In this work we use three values of $t_{s}$ in the case of the connected contributions to assess the influence of excited states. As will be explained, for the case of the disconnected contributions all $t_{s}$ and $t_{\rm ins}$ values are available. We also employ different methods to analyze the ratio of Eq. (10) as explained below. Identifying a time-independent window in this ratio and extracting the desired matrix element by fitting to a constant is referred to as the plateau method. We seek convergence of the extracted value as we increase $t_{s}$.

Instead of using the aforementioned plateau method to extract the matrix element of the ground state, another option is to take into account the contribution of the first excited state. The three-point function can then be expressed as

[TABLE]

while the two-point function is

[TABLE]

$E_{0}(\vec{p})$ and $E_{1}(\vec{p})$ are the energies of the ground state and first excited state at momentum $\vec{p}$, respectively. For non-zero momentum transfer, fitting to the two- and three-point functions taking into account the contribution of the first excited state involves twelve fit parameters, namely $A_{00},\,A_{01},\,A_{10},\,A_{11},\,E_{0}(\vec{p}\,^{\prime}),E_{0}(\vec{p}),\,E_{1}(\vec{p}\,^{\prime}),\,E_{1}(\vec{p}),\,c_{0}(\vec{p}\,^{\prime}),\,c_{0}(\vec{p}),\,c_{1}(\vec{p}\,^{\prime}),$ and $c_{1}(\vec{p})$. We note that $A_{01}\neq A_{10}$ for non-zero momentum transfer. The desired nucleon matrix element $\mathcal{M}$ is obtained via

[TABLE]

In what we refer to as the two-state fit method a simultaneous fit is performed to the three- and two-point functions for several values of $t_{s}$ to obtain $\mathcal{M}$. For the connected three-point function we have three values of $t_{s}$, namely $t_{s}/a=$10, 12 and 14, while for the disconnected we have all values since in our approach the loops are computed for all time slices. We find practical to use a maximal time separation $t_{s}/a$=18, since beyond this separation the correlation functions have large errors and do not contribute to the fit. An alternative technique to study excited state effects is the summation method Maiani et al. (1987); Capitani et al. (2012). Summing over the insertion time $t_{\rm ins}$ of the ratio in Eq. (10) we obtain

[TABLE]

where we omit the source and sink time slices and sum over the geometric series of exponentials. The constant $C$ is independent of $t_{s}$ and $\Delta$ is the energy gap between the first excited state and the ground state, while the matrix element of interest ${\cal M}$ is extracted from a linear fit to Eq. (15) with fit parameters $C$ and ${\cal M}$. Alternatively, as described in Ref. Savage et al. (2016, 1610.04545), one can fit to the finite difference,

[TABLE]

of the summation method, which cancels $C$. We have checked that the two analyses yield consistent results and errors for ${\cal M}$. In the results we quote for the summation method here, we use a linear fit to Eq. (15).

We now briefly describe the so-called fixed-sink method employed to compute the connected contributions to the three-point functions depicted in Fig. 1. Within this method, the sink time, momentum, projector and final and initial hadron states are fixed, but any insertion time-slice and operator with any momentum transfer is allowed, making it the most appropriate method for the study of form factors. An alternative approach is to use stochastic methods to compute the all-to-all quark propagator from the current insertion to the sink, which would allow for both varying the current as well as the sink parameters. This more versatile approach, however, introduces stochastic noise that has to be controlled Alexandrou et al. (2014a). Since in this work we are interested in nucleon form factors, we use the fixed-sink approach for the connected three-point function, which allows for obtaining all insertion momenta with practically no additional computational cost. New sets of inversions are needed for each of the three sink times and each of the three projectors, while the sink momentum is fixed to zero $\vec{p}\,^{\prime}=0$. To increase further our statistics, we average over sixteen randomly selected source positions per gauge configuration.

The disconnected contribution to the three-point function requires the computation of the disconnected quark loop given by

[TABLE]

correlated with the nucleon two-point function. With $M_{f}^{-1}$ we denote the inverse of the twisted mass clover-improved Dirac matrix for the quark flavor $q_{f}$ and with $\mathcal{G}$ a general $\gamma$-structure. For the axial-vector current $\mathcal{G}=\gamma_{\mu}\gamma_{5}$. Eq.(17) requires information from the all-to-all propagator, which is prohibitively expensive to calculate by the standard inversion of the Dirac matrix. For a typical lattice size one needs $\mathcal{O}(10^{8})$ inversions to compute the disconnected quark loop exactly. The standard approach to overcome this difficulty is to use stochastic techniques to obtain an unbiased estimate for the quark loop at the expense of introducing stochastic noise Bitar et al. (1989). Stochastic techniques have been employed successfully in recent studies including our previous studies as for example in Refs. Alexandrou et al. (2012a, 2014b). For certain flavor and operator combinations, such as for example the isoscalar of a flavor doublet of the scalar operator, the twisted mass formulation has a powerful advantage. Such an operator transforms into an isovector of the pseudo-scalar operator in the twisted mass formulation at maximal twist. For the u- and d-flavor doublet we have $\bar{u}u+\bar{d}d=i\bar{\chi}_{u}\gamma_{5}\chi_{u}-i\bar{\chi}_{d}\gamma_{5}\chi_{d}$ where ${\chi}_{u}$ and ${\chi}_{d}$ are the two degenerate light quark fields in the twisted mass basis. The disconnected quark loop contribution to $\sigma_{\pi N}$ therefore becomes Michael and Urbach (2007)

[TABLE]

In other words a subtraction of propagators is replaced by a multiplication resulting in increasing the signal-to-noise ratio from $1/\sqrt{V}$ to $V/\sqrt{V^{2}}$ due to the appearance of an effective double sum over the volume. In this form, stochastic techniques can be employed to obtain the trace via the so-called one-end trick McNeile and Michael (2006) enabling the accurate computation of the quark loops at all time insertions $t_{\rm ins}$ Alexandrou et al. (2014b); Abdel-Rehim et al. (2014b). This method was applied to compute the light, strange and charm $\sigma$-terms with good accuracy Abdel-Rehim et al. (2016). In the case of the axial-vector operator the isoscalar combination does not result into a subtraction in the twisted basis. However, we can generalize the one-end trick to convert the addition of propagators appearing inside a trace into a multiplication. Namely, one can write

[TABLE]

where $D_{WC}$ is the Wilson-Clover operator with bare quark mass set to its critical value. Introducing the stochastic noise vectors $\xi_{r}$ with the properties

[TABLE]

where $N_{r}$ is the number of stochastic vectors, the solution vectors $\phi_{r}=M^{-1}_{u}\xi_{r}$, in Eq. (19) can be written as

[TABLE]

Computing the loop in this way still results in increasing the signal-to-noise ratio from $1/\sqrt{V}$ to $V/\sqrt{V^{2}}$. We refer to the specific application of the trick as in Eq. (21) as the generalized one-end trick, applicable in the case of the axial-vector current where the relative sign between u- and d-quarks does not change in the twisted mass basis. As already pointed out, the one-end trick allows for the evaluation of the quark loops for all insertion time-slices, enabling us to couple them with two-point functions for any value of $t_{s}$ and therefore study thoroughly the excited states behavior.

For computing the strange and charm axial form factors, we use Osterwalder-Seiler Osterwalder and Seiler (1978) valence strange and charm quarks with masses tuned to reproduce the $\Omega^{-}$ and $\Lambda_{c}$ mass respectively. The values we obtain are $a\mu_{s}=0.0259(3)$ and $a\mu_{c}=0.3319(15)$ following the procedure described in Ref. Alexandrou and Kallidonis (2017, 1704.02647). Since we use Osterwalder-Seiler quarks we have the choice to consider doublets with $\pm\mu$-value. We thus construct the axial-vector current as $\frac{1}{2}(\bar{f}^{+}\gamma_{\mu}\gamma_{5}f^{+}+\bar{f}^{-}\gamma_{\mu}\gamma_{5}f^{-})$, where $f=s,c$ and $f^{\pm}$ refers to $\pm\mu_{f}$, yielding the same expressions as for the light quark doublets $(u,d)$ and thus allowing us to apply the generalized one-end trick.

As the pion mass approaches its physical value, the condition number of the Dirac operator increases, hence the conjugate gradient (CG) algorithm requires a larger number of iterations to converge. One can speedup the solver by calculating the lowest eigenvectors of the Dirac operator and then use them to precondition the conjugate gradient (CG) algorithm, by deflating the Dirac operator. In our calculations, we use the implicitly restarted Lanczos algorithm to calculate the eigenvectors. We found that deflating 600 eigenvectors results in a factor of about $20\times$ speedup for the light quark masses as shown in Fig. 2. For the light fermion loops we calculate 2250 stochastic noise vectors per configuration to high precision (HP), i.e. to a solver precision of 10*-9*. Note that no dilution has been employed, therefore one inversion per noise source is performed.

For the strange and the charm quarks the condition number of the Dirac operator is significantly smaller, thus there is no need for deflation. Instead, we employ the truncated solver method (TSM) Bali et al. (2010) where a large number of low precision (LP) noise vectors is used to reduce the stochastic variance and the bias is corrected by a small number of HP inversions. The number of iterations for a LP solve ($n_{\rm LP}$), as well as, the number of low ($N^{\rm LP}_{r}$) and high ($N^{\rm HP}_{r}$) precision inversions needs to be tuned in order to produce an unbiased estimate of the disconnected quark loop at optimal computational cost. Namely, the variance $\sigma^{2}$ can be approximated by Blum et al. (2013):

[TABLE]

where $r_{c}$ is the correlation between the targeted observable computed to high and low precision. A compromise is necessary that keeps the ratio $N^{\rm HP}_{r}/N^{\rm LP}_{r}$ small, while having $r_{c}\simeq 1$. For the strange and charm loops used in this work, we take $n_{\rm LP}$ such that we obtain $r_{c}\simeq 0.99$. We investigate the dependence of $r_{c}$ on $n_{\rm LP}$ in the left panel of Fig. 3, for various values of the twisted mass parameter. One can see that as $\mu$ decreases a larger number of iterations is needed to obtain the same value for $r_{c}$. For $a\mu=0.001$, which is very close to our value of $a\mu_{l}=0.0009$, the number of iterations needed to reach a good correlation is very large, indicating that the TSM is not efficient for light quark masses. In the right panel of Fig. 3 we show the number of iterations needed to have $r_{c}\simeq 0.99$ for a given bare quark mass. This figure shows that about $\approx 100$ iterations are sufficient for the case of the strange quark mass $a\mu_{s}\simeq 0.03$, resulting in a solver precision of $\simeq 10^{-3}$. With the values of $r_{c}$ and $n_{\textrm{LP}}$ at hand, we use:

[TABLE]

from Ref. Bali et al. (2010) to determine the ratio $N^{\rm LP}_{r}/N^{\rm HP}_{r}$. Eq. (23) is obtained by requiring minimization of the stochastic variance for equal cost. For the strange quark we test that there is no bias by increasing the resulting $N^{\rm HP}_{r}$ and observing whether the central value of our observable changes. For the charm quark, the inverter reaches very quickly our target $r_{c}$ after $\mathcal{O}(10)$ iterations. We therefore increase $n_{\rm LP}$ to yield $r_{c}\simeq 0.999$ since this increases minimally the total computational cost. This allows us to use a larger value for the $N^{\rm LP}_{r}/N^{\rm HP}_{r}$ ratio for the charm loops.

The statistics used and the parameters used for the TSM for the strange and the charm quark loops are listed in Table 2, with the disconnected fermion loops calculated for all time-slices. For the connected three-point functions three source-sink time separations have been analyzed for 16 source-positions per gauge configuration, while for the two-point functions 100 source-positions per gauge configuration have been produced in order to accumulate enough statistics for the disconnected three-point function.

V Renormalization

In order to make a comparison of form factors calculated from lattice QCD with experimental and phenomenological results, one must renormalize the lattice results. The renormalization functions can be calculated perturbatively as well as non-perturbatively. In this work, we use the non-perturbatively calculated renormalization functions Alexandrou et al. (2012b) where lattice artifacts are computed perturbatively Alexandrou et al. (2017) and subtracted from the non-perturbative results before taking the continuum limit. The Rome-Southampton method Martinelli et al. (1995), also known as the ${\rm RI^{\prime}_{MOM}}$ scheme, is used for the calculation of the renormalization functions. Note that the renormalization function $Z_{A}$ for the axial current is scheme and scale independent in the chiral limit.

In the case of flavor non-singlet operators such as the isovector axial operator, the renormalization functions can be calculated accurately with a relatively low cost, whereas the isoscalar combination receives contributions from a disconnected diagram, leading to significant increase in the computational effort. In order to calculate the renormalization functions non-perturbatively, we consider the bare vertex functions Gockeler et al. (1999)

[TABLE]

where $G^{ns}_{\mathcal{G}}$ and $G^{s}_{\mathcal{G}}$ are the non-singlet and singlet cases, respectively, $V$ is the lattice volume and, in our case, $\mathcal{G}=\gamma_{\mu}\gamma_{5}$. We employ the momentum source method which offers a high statistical accuracy. In particular, statistical errors are of the order of ‰ with ${\cal O}(10)$ measurements. The amputated vertex function can be derived from the vertex function as

[TABLE]

where $S(p)$ is the propagator in momentum space. For the singlet vertex function the disconnected contribution is amputated using one inverse propagator as the closed quark loop does not have an open leg.

In the ${\rm RI^{\prime}_{MOM}}$ scheme the renormalization functions are computed by imposing that the amputated vertex function $\Lambda_{\mathcal{G}}(p)$ at large Euclidean scale $p^{2}=\mu^{2}$, is equal to its tree-level value in the chiral limit. The renormalization condition is given by

[TABLE]

The non-singlet renormalization functions for the ensemble used in this work can be found in Ref. Alexandrou et al. (2017). In Fig. 4 we show our results for the axial singlet renormalization function for three pion masses and for several initial momenta. As can been seen, the dependence on the light quark mass is very mild. In Fig. 5 we show an example of a chiral extrapolation we perform at $(ap)^{2}=2$, which corroborates that the pion mass dependence is very weak. In Fig. 4 we also show the difference between the singlet and the non-singlet case for different pion masses and $(ap)^{2}$. We observe a small but non-zero difference. The chirally extrapolated values are shown in Fig. 5 and they are used to perform the continuum limit. In general, the momentum source method leads to small statistical errors and thus a careful investigation of systematic uncertainties is required. We eliminate the systematic effect that comes from the asymmetry of our lattices, such as due to the larger time extent and the antiperiodic boundary conditions in time, by averaging over the different components corresponding to the same renormalization function. Furthermore, remaining lattice artifacts are partially removed by the subtraction of the ${\cal O}(g^{2}\,a^{\infty})$ terms as was done in Refs. Constantinou et al. (2009); Alexandrou et al. (2017). However, the largest systematic error comes from the choice of the momentum range to use for the extrapolation to $(ap)^{2}\to 0$. To address this effect we use different intervals for the $(ap)^{2}\to 0$ fit and obtain the systematic error, shown by the black error bar in Fig. 5, by taking the largest difference in the values of the renormalization function extracted from different fit ranges. We find for the non-singlet operator that $Z_{A}^{ns}=0.7910(4)(5)$, as was originally reported in Ref. Alexandrou et al. (2017), while for the singlet $Z_{A}^{s}=0.7968(25)(91)$. Due to the large systematic error $Z_{A}^{ns}$ and $Z_{A}^{s}$ are compatible.

VI Results

VI.1 Axial charge

We first examine the extraction of the axial charge of the nucleon, which is given by $g_{A}\equiv G_{A}^{u-d}(0)$. In order to assess the effect of the excited states we study the ratio of Eq. (10) for various source-sink time separations. In Figs. 6 and 7 we show the ratio from which we extract the nucleon isovector axial charges $g_{A}$ and the isoscalar $g_{A}^{u+d}$ including the disconnected contribution. We also show the corresponding ratios from where $g_{A}^{s}$ and $g_{A}^{c}$ are determined.

In the case of zero momentum transfer the square root of Eq. (10) reduces to unity, and the matrix element of Eq. (2) directly yields the axial charge. In Fig. 6 we show the ratio of Eq. (10) for various values of $t_{s}$ as a function of the insertion time. The values extracted from the plateau, summation and two-state fits are collected in the right panel of the figure. As can be seen, as $t_{s}$ increases the plateau value converges to a constant indicating that excited states become very small. When the plateau value is in agreement with the value extracted from the two-state fit we consider that contributions from excited states are sufficiently suppressed. We take the plateau value for the smallest $t_{s}$ where agreement with the two-state fit is observed as our final value for the matrix element. This value is always consistent with the result from the summation method since the statistical error of the latter is usually larger as compared to the two-state fit. As a systematic error due to excited states we take the difference between the plateau value that demonstrates convergence with $t_{s}$ and that extracted from the two-state fit.

As can be clearly seen from Fig. 7, the disconnected contributions are non-zero and negative. The value of $g_{A}^{s}$ is smaller as compared to the disconnected contribution to $g_{A}^{u+d}$. $g_{A}^{c}$, although still negative, has a large error and a small value, namely $|g_{A}^{c}|<0.005$. We note here that the value of the disconnected contribution to $g_{A}^{u+d}$ extracted from our previous study Abdel-Rehim et al. (2014b) using a TMF ensemble simulated at a pion mass of $m_{\pi}=370$ MeV is about twice smaller, namely -0.07(1), compared to the physical point value obtained here. Lattice artifacts for nucleon observables such as the ones calculated here are expected to be small. A comparison of results for the axial charge from various lattice actions including $N_{f}=2$, $N_{f}=2+1$ and $N_{f}=2+1+1$ flavors of quarks, as well as various lattice spacings and volumes shows that volume, cut-off and strange quark quenching effects are smaller than current statistical errors Alexandrou et al. (2017, 1706.02973). In Fig. 8, we show a comparison of lattice results for $g^{A}_{s}$. In particlar, we compare results using $N_{f}=2$ clover fermions at a pion mass of about $300$ MeV clover fermions from Ref. Bali et al. (2012) with results using domain wall valence fermions on $N_{f}=2+1$ asqtad gauge configurations (hybrid action) from Ref. Engelhardt (2012). Both $N_{f}=2$ and $N_{f}=2+1$ results are compatible indicating that strange sea quark effects and lattice artifacts are small compared with the statistical errors. The $N_{f}=2+1$ hybrid action result at about 370 MeV is also in agreement with the $N_{f}=2+1+1$ twisted mass fermion result, which was a high accuracy computation. Since we do expect charm quark effects to be negligible, this agreement corroborates between calculations with different actions corroborates that lattice artifacts are indeed smaller that the current statistical errors.

Our values for the nucleon axial charges are tabulated in Table 3. In the case of $g_{A}$ our result is compatible with recent results from the lattice Berkowitz et al. (2017, 1704.01114); Bhattacharya et al. (2016); Green et al. (2014); Capitani et al. (2017); Bali et al. (2015); Horsley et al. (2014) and slightly underestimates the experimental value of 1.2723(23) Patrignani et al. (2016). In the case of $g_{A}^{u+d}$ there is a good agreement with the experimental value of 0.416(18) Patrignani et al. (2016) within the current statistics.

VI.2 Axial and induced pseudo-scalar form factor

For non-zero momentum transfer, both $G_{A}$ and $G_{p}$ enter in Eq. (2), namely the large time limit of Eq. (10) is related to the form factors via

[TABLE]

in the case where $i=k$, and

[TABLE]

for $i\neq k$ with

[TABLE]

Since the form factors depend only on the momentum transfer squared ($Q^{2}$), while the plateau of Eq. (25) depends on $\vec{p}\,^{\prime}$ and $\vec{p}$, the extraction of the form factors is over-constrained. In practice, we form the system:

[TABLE]

where $D$ is an array of kinematic coefficients according to Eqs. (27, 28) and $F$ is the vector: $F^{\intercal}=(G_{A},G_{p})$. The system is solved for $F$ by taking the Singular Value Decomposition (SVD) of $D$ in order to minimize:

[TABLE]

for each $Q^{2}$, where $w$ is the statistical error of $\Pi$.

All results quoted in this paper are computed by first fitting the ratio $R_{i}(\Gamma_{k},\vec{p}\,^{\prime},\vec{p};t_{s},t_{\rm ins})$ with either the plateau, two-state or summation method to obtain $\Pi_{i}(k,\vec{p}\,^{\prime},\vec{p})$ and subsequently minimize Eq. (31) to obtain $F(Q^{2})$ without a time dependence. In order to demonstrate these plateaus we carry out an analysis in a different order. Namely we apply the SVD and minimization of Eq. (31) by inserting the ratio $R_{i}(\Gamma_{k},\vec{p}\,^{\prime},\vec{p};t_{s},t_{\rm ins})$ instead of $\Pi_{i}(k,\vec{p}\,^{\prime},\vec{p})$. In Figs. 9 and 10 we show representative examples of our obtained plateaus for a small momentum transfer, namely for $Q^{2}=0.0753\;\mathrm{GeV^{2}}$ and for a higher momentum transfer, namely $Q^{2}$=0.2848 GeV2. The corresponding results for the form factors are shown in Figs. 11 and 12 for the same momentum transfers as for Figs. 9 and 10. We observe a similar behavior with respect to the excited states as that for $Q^{2}=0$ shown in Figs. 6 and 7. We thus take the plateau value at $t_{s}=1.31$ fm for both $G_{A}^{u-d}$ and $G_{A}^{u+d}$ as our final values since they are in good agreement with the two-state and summation fits. For $G_{A}^{s}$, $t_{s}=1.31$ fm is still a reasonable choice, but for $G_{A}^{c}$ due to the large statistical uncertainty, $t_{s}=0.94$ fm is enough.

For $G_{p}^{u-d}$ and the connected part of $G_{p}^{u+d}$ we observe excited state contributions for the two smaller values of $t_{s}$. For $t_{s}=1.31$ fm the plateau value is in agreement with the two-state fit, however, indicating partial convergence. For the disconnected contribution we find better convergence and we take the value also at $t_{s}=1.31$ fm. What is particularly notable are the large disconnected contributions to the isoscalar induced pseudo-scalar form factor that are comparable in magnitude to the connected part, but with opposite sign. This has already been observed in Ref. Green et al. (2017), which used an ensemble simulated with a pion mass $m_{\pi}=317$ MeV. The explanation of such large disconnected contributions is that they are needed to cancel the pion pole of the connected isoscalar form factor in order to yield the expected $\eta$-meson pole mass dependence. Since the connected isoscalar shows a sharp rise consistent with a pion pole the disconnected contributions must also be large at small $Q^{2}$ to cancel it. This would be analogous to the case of the $\eta$-meson mass extraction on the lattice, where disconnected contributions are important since the connected contribution alone of the two-point correlation function has the mass of the pion as ground state Alexandrou et al. (2012a). From Fig. 12 where results are shown for a relatively high $Q^{2}$ value, the overall observation is that excited state contributions tend to be less severe but non-negligible. This trend continues as we increase $Q^{2}$ at least for the connected contributions where statistical uncertainties are small enough for such an investigation.

In Fig. 13 we show the isovector form factors up to $Q^{2}=1\;\mathrm{GeV^{2}}$ extracted from the plateau at the three values of $t_{s}$ considered, from the two-state and summation methods Alexandrou et al. (2016). As already noted, for $G_{p}^{u-d}(Q^{2})$, excited states contributions are notably more severe for small values of $Q^{2}$, which tend to decrease its value. Nevertheless, the values extracted from the plateau at $t_{s}=1.31$ fm are in agreement with the value extracted from the two-state fit for all $Q^{2}$-values. We thus take the plateau value at $t_{s}=1.31$ fm as our final value for the form factors with a systematic error the difference between the mean value from fitting the plateau at $t_{s}=1.31$ fm and that extracted from the two-state fit. This systematic error may be underestimated for $G_{p}(Q^{2})$ at low $Q^{2}\leqslant 0.2$ GeV2 where a larger time separation may be needed to ensure convergence.

Having a determination of the axial form factors we proceed to examine their $Q^{2}$-dependence. As customarily done in experiment we fit the axial form factor $G_{A}(Q^{2})$ to a dipole form given by

[TABLE]

where $m_{A}$ is the so-called axial mass and the axial radius, $\langle r_{A}^{2}\rangle$, is related to $m_{A}$ by

[TABLE]

We note that experimentally, one of the determinations of $m_{A}$ is obtained by fitting the axial form factor $G_{A}(Q^{2})$ extracted from pion electroproduction data, yielding a value of $m_{A}=1.077(39)$ GeV Liesenfeld et al. (1999). Recent results from charged-current muon-neutrino scattering events produced from the MiniBooNE experiment report a value of $m_{A}=1.350(170)$ GeV using a similar fit Aguilar-Arevalo et al. (2010), which is significantly higher than the historical world average. Recent results from neutrino-nucleus cross sections using deuterium target data report a smaller value of $m_{A}=1.010(240)$ GeV Meyer et al. (2016).

Fitting the momentum dependence of our results for $G_{A}^{u-d}(Q^{2})$ using Eq. (32) we obtain a value of $m_{A}=1.322(42)(17)$ GeV, which is consistent with the larger value extracted from $\nu_{\mu}$-interactions Aguilar-Arevalo et al. (2010). The fit is performed by fixing the value for $g_{A}$ directly from our lattice result for $G_{A}^{u-d}(0)$. We have checked that allowing $g_{A}$ to vary as a fit parameter yields consistent results. We also extract a consistent value for $m_{A}$ using the results from the two-state fit. We quote the difference in the mean value of $m_{A}$ extracted from fitting the $t_{s}=1.31$ fm plateau results for the form factors and that extracted from the results of the two-state fits as the systematic error due to excited states. In the left panel of Fig. 14 we show a comparison of the fits to our lattice QCD results and the experimental ones. The spread in the mean values is an indication of remaining excited state contributions, which are small and which we quote as our systematic error. The bands produced using the values from Refs. Liesenfeld et al. (1999); Meyer et al. (2016) are lower and have a steeper slope than our results, albeit with large errors.

In Fig. 14 we show our lattice QCD results for $G_{p}^{u-d}(Q^{2})$. As expected from pion pole dominance, this form factor has a much stronger $Q^{2}$ dependence as compared to $G_{A}^{u-d}(Q^{2})$. Using the partially conserved axial current relation (PCAC) and pion pole dominance one can relate the induced pseudo-scalar form factor $G_{p}^{u-d}(Q^{2})$ to $G_{A}^{u-d}(Q^{2})$ by

[TABLE]

where $C=4m_{N}^{2}$ and $m_{p}=m_{\pi}$. This relation is used to extract the induced pseudo-scalar form factor using the experimental determination of $G_{A}^{u-d}(Q^{2})$. We perform the same analysis for our lattice QCD results. Namely, in Fig. 14 we include results for $G_{p}^{u-d}(Q^{2})$ obtained by applying the pion-pole dominance hypothesis to the lattice results on $G_{A}^{u-d}$ using the lattice pion mass of $m_{\pi}=130$ MeV in Eq. (34). At low $Q^{2}$-value we observe a much steeper rise as compared to the direct lattice computation of $G_{p}(Q^{2})$, and agreement both with the experimentally determined bands taken by applying the pion-pole assumption as well as with the directly determined values of $G_{p}^{u-d}(Q^{2})$ from Ref. Choi et al. (1993). As noted above, for $Q^{2}<0.2$ GeV2 where the discrepancy is largest, excited states tend to produce smaller values. In addition, similar discrepancy at low $Q^{2}$ has been observed in previous lattice studies at heavier pion masses between multiple volumes Alexandrou et al. (2011, 2007), indicating that volume effects may also need to be investigated to resolve this tension. We plan to investigate both these systematics using a larger volume of $64^{3}\times 128$ in a future work. For the current analysis we will discard $G_{p}^{u-d}(Q^{2})$ at the two lowest values of $Q^{2}$.

In addition to the dipole form, we fit our results for the axial form factor using the so-called z-expansion Hill and Paz (2010), given by

[TABLE]

where

[TABLE]

and $t_{cut}=9m_{\pi}^{2}$. In Fig. 15 we compare the dipole fit with the z-expansion fit. For the z-expansion we used $k_{max}=3$, fixing $a_{0}=g_{A}$ and imposing Gaussian priors for the coefficients $a_{k}$ for $k>1$ with width $w=5\mathrm{max}(|a_{0}|,|a_{1}|)$. Both fit Ansätze describe the data very well, producing consistent values for the radius, namely $\langle r_{A}^{2}\rangle=0.266(17)$ fm2 in the case of the dipole fit and $\langle r_{A}^{2}\rangle=0.265(76)$ fm2 from the z-expansion. A fit using the z-expansion is more suitable when precise data are available at a large number of $Q^{2}$ values. Given the statistical errors and relatively few momenta available from our lattice calculation, the z-expansion therefore yields larger errors than a dipole fit. Given the consistency between the two fits we thus opt to use the dipole form that yields smaller errors for all the fits that follow.

PCAC relates the residue of the pion pole to the pion decay constant $f_{\pi}$, the nucleon mass $m_{N}$ and the pion-nucleon coupling constant $g_{\pi NN}$ as follows Goldberger and Treiman (1958)

[TABLE]

The relation holds when including the leading correction as obtained within the chiral perturbative framework used in Ref. Bernard et al. (1995). Using Eq. (34) we can relate $g_{\pi NN}$ to the axial form factor as

[TABLE]

where for this ensemble $f_{\pi}=89.80$ MeV Abdel-Rehim et al. (2017) and $m_{N}$=0.932(4) GeV Alexandrou and Kallidonis (2017, 1704.02647). Using $G_{A}^{u-d}(-m_{\pi}^{2})$=1.234(35)(20) obtained from our dipole fit and $C=4m_{N}^{2}$, we find $g_{\pi NN}=12.81(37)(21)$, which is consistent with the experimental value $g_{\pi NN}=13.12(10)$ measured from pion-nucleon scattering lengths Baru et al. (2011). Were we to fit directly the lattice data for $G_{p}^{u-d}(Q^{2})$ to the form

[TABLE]

taking $m_{\pi}=130$ MeV and omitting the first two $Q^{2}$ values from the fit we obtain the solid line in Fig. 14, for which $m_{p}=1.441(115)(648)$ GeV consistent with the axial mass from fitting to the axial form factor and $G^{u-d}_{p}(0)=165.62(9.82)(18.46)$ which is smaller than $4(m_{N}^{2}/m_{\pi}^{2})g_{A}$. If we then were to use Eq. (37) we would determine $g_{\pi NN}$=8.50(51)(82). This is smaller than the value determined using pion-pole dominance and our lattice results for $G_{A}^{u-d}$. Additionally one can compute also the induced pseudoscalar charge, $g_{p}^{*}$, defined as

[TABLE]

where $m_{\mu}$ is the muon mass. We find $g_{p}^{*}=7.47(30)(80)$ using our lattice results for $G_{A}^{u-d}$ and pion-pole dominance. Using the lattice results for $G_{p}^{u-d}$ the value of $g_{p}^{*}$ is lower and calls for a further study of excited states and volume effects on the lattice determination of $G_{p}^{u-d}(Q^{2})$.

In order to compute the individual light quark axial form factors one needs, besides the isovector form factors, the isoscalar combination. In Fig. 16 we illustrate our results for the connected contributions to $G_{A}^{u+d}(Q^{2})$ and $G_{p}^{u+d}(Q^{2})$ using the same analysis as for the isovector. Once more, excited states are clearly more severe for $G_{p}^{u+d}(Q^{2})$ at low $Q^{2}$ where the pion pole dominates and tends to decrease its value leading to a milder $Q^{2}$-dependence.

In Fig. 17 we show the disconnected contributions to $G_{A}^{u+d}(Q^{2})$, which are clearly non-zero and negative. The form factors for the disconnected contributions are obtained combining final nucleon states with $\vec{p}^{\prime}=\vec{0}$, the same as in the case of the connected contributions, and in addition all sink momenta which satisfy $\vec{p}^{\prime 2}=(2\pi/L)^{2}$. Since more $Q^{2}$ values are available, we plot, in Fig. 17, the sink-source separation $t_{s}=1.31$ fm and two-state fit methods alone for better clarity. The disconnected contributions reduce the value of $G_{A}^{u+d}(Q^{2})$ and for zero momentum transfer result in a value compatible with the experimental one. As already mentioned, the disconnected contributions to $G_{p}^{u+d}(Q^{2})$ are particularly large and reduce its value especially at low values of $Q^{2}$. Adding the connected and disconnected contributions obtained using $\vec{p}\,^{\prime}=\vec{0}$ for which common $Q^{2}$-values are available, yields the result shown in Fig. 18. We note that, due to the fact that the disconnected part is computed with much higher statistics as compared to the connected, the error in the total quantity is computed by adding the individual errors in quadrature. In Fig. 19 we show the resulting dipole fits to the isoscalar form factor $G_{A}^{u+d}(Q^{2})$ using Eq. (32) for the connected, disconnected and total value. The parameters extracted are collected in Table 4. The axial mass extracted by fitting $G_{A}^{u+d}(Q^{2})$ is $m_{A}^{u+d}=1.736(244)(374)$ GeV. Although the central value is larger, within the large statistical and systematic errors it is in agreement with the one extracted for the isovector case. In Table 4 we also list the corresponding axial radii, obtained from the dipole masses via Eq.(33).

For $G_{p}(Q^{2})$ we fit using Eq. (34) for the connected and disconnected separately, allowing $C$ and $m_{p}$ to vary. We obtain the curves shown in Fig. 19 and consistent pole masses, namely $m_{p}^{\textrm{u+d,conn}}=0.324(22)(12)$ GeV for the connected and $m_{p}^{\textrm{u+d,disc}}=0.331(81)(36)$ GeV for the disconnected. We show the total isoscalar $G_{p}(Q^{2})$ in Fig. 19. As can be seen the errors are large, especially in the small $Q^{2}$ region, and do not allow us to reliably quote a value for the pole mass of $G_{p}(Q^{2})$.

In Fig. 20 we show the strange and the charm form factors, which only take disconnected contributions. For the strange quark contributions, we use sink momenta $\vec{p}\,^{\prime}=\vec{0}$ and $\vec{p}\,^{\prime 2}=(2\pi/L)^{2}$, while for the charm quark where errors are large only the $\vec{p}\,^{\prime}=\vec{0}$ case yields reasonable results. We observe a very good signal for $G_{A}^{s}(Q^{2})$ up to momentum transfer $Q^{2}=0.5$ GeV2. As already noted our results from the plateau method with $t_{s}=1.31$ fm are consistent with the two-state fit. $G_{A}^{s}(Q^{2})$ can be well fitted to a dipole form and we obtain $m_{A}^{s}=0.921(228)(90)$ GeV. The results for $G_{A}^{s}(Q^{2})$ are compatible with the experimental values measured for $Q^{2}>0.45$ GeV2, which however carry large errors Pate et al. (2008). $G_{A}^{c}(Q^{2})$ is more noisy in particular for the larger time separations. For the smallest source-sink separation of $t_{s}=0.94$ fm we obtain a non-zero negative value for the whole range of $Q^{2}$. However, for larger values of $t_{s}$ the results become noisy forbidding us to reach a conclusion on excited states contributions. For $G_{p}^{s}(Q^{2})$ we obtain a non-zero negative contribution, which is about six times smaller in magnitude compared to the disconnected $G_{p}^{u+d}(Q^{2})$. In the case of $G_{p}^{c}(Q^{2})$ results are compatible with zero even for the smallest source-sink time separation. We do not display the results produced with the summation method since these are very noisy.

In Fig. 21 we show fits to $G_{A}^{s}(Q^{2})$ using the dipole form of Eq. (32) and $G_{p}^{s}(Q^{2})$ using Eq. (34). For $G_{p}^{c}(Q^{2})$ the results are compatible with zero and no fit is attempted. Axial masses and corresponding radii extracted from the dipole fits are tabulated in Table 4.

VI.3 Comparison with other studies

The axial and induced pseudo-scalar form factors have been studied by several lattice QCD groups using recent dynamical simulations. Preliminary lattice QCD results using an ensemble with close to physical pion mass has been presented by the PNDME collaboration Jang et al. (2016). They use a mixed action approach of $N_{f}=2+1+1$ HISQ staggered fermions and clover-improved Wilson valence fermions. This action has ${\cal O}(a)$ lattice artifacts, which are shown to be sizeable for $a=0.09$ fm as compared to their results at $a=0.06$ fm. Their preliminary results on $G_{A}^{u-d}$ using an ensemble at pion mass $m_{\pi}=130$ MeV and $a=0.06$ fm are in agreement with ours. This shows that lattice artifacts for our ${\cal O}(a)$ improved action computed with $a=0.0938$ fm are small. On the other hand, their results on $G_{p}^{u-d}(Q^{2})$ for the same ensemble are larger at low $Q^{2}$-values than ours. Given that their spatial box length is $L\sim 5.76$ fm as compared to $L\sim 4.51$ fm of our lattice, these preliminary results may indicate that $G_{p}^{u-d}(Q^{2})$ suffers from sizeable finite volume effects. Additional lattice QCD results on the isovector axial form factors at higher than physical pion mass have been computed recently by two groups: LHPC has obtained results on the isovector axial form factors using $N_{f}=2+1$ clover-improved Wilson fermions with $m_{\pi}=317$ MeV Green et al. (2017), which includes the isoscalar form factors and using a mixed action for $m_{\pi}=356$ MeV Bratt et al. (2010). CLS has presented preliminary results using an ensemble of $N_{f}=2$ clover fermions at a pion mass of $m_{\pi}\sim 340$ MeV von Hippel et al. (2016). In what follows we restrict ourselves to showing published results only.

In Fig. 22 we compare our results for the isovector axial form factors to the published LHPC results, which have been produced using $N_{f}=2+1$ Asqtad staggered sea quarks on a $28^{3}\times 64$ lattice and Domain-Wall valence fermions for $m_{\pi}=356$ MeV Bratt et al. (2010). Our results for $G_{A}^{u-d}$ at the physical point show a steeper $Q^{2}$-dependence leading to a larger value of $g_{A}$. For $G_{p}^{u-d}(Q^{2})$, the LHPC results tend to be larger in particular at the smallest $Q^{2}<0.2$ GeV2. The length of their lattice is $L=3.36$ fm yielding $Lm_{\pi}\sim 6$ as compared to $Lm_{\pi}\sim 3.3$ fm for our lattice. This may again point to finite volume effects that need to be investigated.

VII Conclusions

Results on the nucleon axial form factors are presented for one ensemble of two degenerate twisted mass clover-improved fermions tuned to reproduce approximately the physical value of the pion mass. Using improved techniques we evaluate both connected and disconnected contributions to both axial and induced pseudo-scalar form factors. Our study includes an investigation of excited state effects by computing the nucleon three-point functions at several sink-source time separations. Lattice matrix elements are non-perturbatively renormalized by computing both the singlet and non-singlet renormalization functions.

We find that the isovector axial form factor $G_{A}^{u-d}(Q^{2})$ is described well by a dipole form with an axial mass $m_{A}=1.322(42)(17)$ GeV, which is larger than the historical world average but is in agreement with a recent value produced from the MiniBooNE experiment Aguilar-Arevalo et al. (2010). We can relate, via PCAC, the axial form factor to the $\pi N$ coupling constant. We find that $g_{\pi NN}=12.81(37)(21)$ consistent with the experimental value of $g_{\pi NN}=13.12(10)$ Baru et al. (2011). Similarly, we can deduce $G^{u-d}_{p}(Q^{2})$ from our results on $G_{A}^{u-d}(Q^{2})$ assuming pion pole dominance yielding agreement with experiment. However, a direct extraction of the isovector induced pseudo-scalar form factor has a weaker $Q^{2}$-dependence as compared to what is expected from pion pole dominance. Thus, although one can describe well the data using a pion pole behavior for its $Q^{2}$-dependence, one extracts a pole mass larger than the ensemble value of $m_{\pi}=130$ MeV. $G_{p}(Q^{2})$ at low $Q^{2}$ is shown to have more severe excited states effects, which tend to lower its value. Comparison to preliminary lattice results obtained on a larger volume Jang et al. (2016) indicate that volume effects may also increase its value at low $Q^{2}$. We plan to check for such volume effects in a future analysis using a larger lattice.

An important conclusion of this work is that disconnected contributions to both isoscalar and strange form factors are non-negligible. For the isoscalar $g_{A}^{u+d}$ these contributions need to be taken into account to bring agreement with the experimental value. For $G_{p}^{u+d}(Q^{2})$ the disconnected contributions are particularly large and of the same order as the connected part but with the opposite sign leading to a weaker $Q^{2}$-dependence for the isoscalar pseudo-scalar form factor. Both strange from factors $G_{A}^{s}(Q^{2})$ and $G_{p}^{s}(Q^{2})$ are found to be negative and non-zero, with the magnitude of $G_{A}^{s}(Q^{2})$ of the same order as that for the light disconnected contributions. Both charm from factors tend to be negative but given the large errors they remain compatible with zero.

VIII ACKNOWLEDGMENTS

We would like to thank the members of the ETMC for a most enjoyable collaboration. We acknowledge funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 642069. This work was partly supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project IDs s540 and s625 on the Piz Daint system, by a Gauss allocation on SuperMUC with ID 44060 and in addition used computational resources from the John von Neumann-Institute for Computing on the Jureca and the BlueGene/Q Juqueen systems at the research center in Jülich. We also acknowledge PRACE for awarding us access to the Tier-0 computing resources Curie, Fermi and SuperMUC based in CEA, France, Cineca, Italy and LRZ, Germany, respectively. We thank the staff members at all sites for their kind and sustained support. K. H. and Ch. K. acknowledge support from the Cyprus Research Promotion Foundation under contract T$\Pi$E/$\Pi\Lambda$HPO/0311(BIE)/09.

Appendix A: Table of Results