Up, down, and strange nucleon axial form factors from lattice QCD

TL;DR

This paper presents a lattice QCD calculation of the nucleon axial form factors for all three light quark flavors, providing insights into nucleon structure and related physical constants.

Contribution

It introduces a comprehensive lattice QCD approach to compute all light-flavor nucleon axial form factors, including disconnected contributions and nonperturbative renormalization.

Findings

Disconnected contributions to $G_P$ are large.

Quark contributions to nucleon spin are quantified.

Induced pseudoscalar coupling and pion-nucleon coupling are estimated.

Abstract

We report a calculation of the nucleon axial form factors $G_A^q(Q^2)$ and $G_P^q(Q^2)$ for all three light quark flavors $q\in\{u,d,s\}$ in the range $0\leq Q^2\lesssim 1.2\text{ GeV}^2$ using lattice QCD. This work was done using a single ensemble with pion mass 317 MeV and made use of the hierarchical probing technique to efficiently evaluate the required disconnected loops. We perform nonperturbative renormalization of the axial current, including a nonperturbative treatment of the mixing between light and strange currents due to the singlet-nonsinglet difference caused by the axial anomaly. The form factor shapes are fit using the model-independent $z$ expansion. From $G_A^q(Q^2)$, we determine the quark contributions to the nucleon spin and axial radii. By extrapolating the isovector $G_P^{u-d}(Q^2)$, we obtain the induced pseudoscalar coupling relevant for ordinary muon capture and the pion-nucleon coupling constant. We find that the disconnected contributions to $G_P$ form factors are large, and give an interpretation based on the dominant influence of the pseudoscalar poles in these form factors.

Tables

Table 1. Table 1: Number of Hadamard vectors used for disconnected loops of each flavor and source-operator separation τ 𝜏 \tau . Five separations were used for light quarks and three for strange. As shown in Subsec.  II.4 , sixteen Hadamard vectors is generally sufficient for the noise to saturate when using the axial current. Having 128 Hadamard vectors was particularly useful for Ref.  Green:2015wqa , which used the vector current.

$\tau /a=$	3	4	5	6	7
light	16	128	128	128	16
strange		16	128	16

Table 2. Table 2: Quark spin contributions to the nucleon spin.

$q$	$g_A^q$
$u$	$0.863\left(7\right)\left(14\right)$
$d$	$-0.345\left(6\right)\left(9\right)$
$s$	$-0.0240\left(21\right)\left(11\right)$

Table 3. Table 3: Parameters for the z 𝑧 z -expansion fits, along with their correlation matrices. The latter are symmetric, and we have omitted redundant entries.

$G_A^{u-d}\left(Q^2\right)$
k	$a_k$	Correlation matrix
0	$1.208\left(20\right)$	$1$	$-0.771$	$0.365$	$-0.101$	$0.067$	$0.058$
1	$-3.985\left(332\right)$		$1$	$-0.767$	$0.509$	$0.229$	$0.168$
2	$0.877\left(1.639\right)$			$1$	$-0.911$	$-0.471$	$-0.299$
3	$7.730\left(4.783\right)$				$1$	$0.416$	$0.201$
4	$4.324\left(3.101\right)$					$1$	$0.963$
5	$1.615\left(1.417\right)$						$1$
$G_A^{u+d}\left(Q^2\right)$
k	$a_k$	Correlation matrix
0	$0.517\left(18\right)$	$1$	$-0.712$	$0.430$	$-0.215$	$-0.271$	$-0.323$
1	$-1.582\left(274\right)$		$1$	$-0.821$	$0.427$	$0.459$	$0.467$
2	$0.947\left(1.975\right)$			$1$	$-0.822$	$-0.826$	$-0.803$
3	$-0.853\left(5.519\right)$				$1$	$0.990$	$0.965$
4	$-0.534\left(2.451\right)$					$1$	$0.991$
5	$-0.214\left(745\right)$						$1$
$G_A^s\left(Q^2\right)$
k	$a_k$	Correlation matrix
0	$-0.0240\left(24\right)$	$1$	$-0.678$	$0.478$	$0.314$	$0.177$	$0.094$
1	$0.0577\left(386\right)$		$1$	$-0.943$	$-0.802$	$-0.656$	$-0.575$
2	$0.0274\left(1445\right)$			$1$	$0.845$	$0.687$	$0.598$
3	$-0.0079\left(507\right)$				$1$	$0.963$	$0.911$
4	$-0.0049\left(156\right)$					$1$	$0.987$
5	$-0.0017\left(43\right)$						$1$
$a^2\left(Q^2+m_{\pi }^2\right)G_P^{u-d}\left(Q^2\right)$
k	$a_k$	Correlation matrix
0	$1.613\left(174\right)$	$1$	$-0.882$	$0.518$	$0.052$	$0.056$	$0.107$
1	$-3.997\left(1.583\right)$		$1$	$-0.838$	$0.086$	$0.118$	$0.078$
2	$-3.946\left(5.729\right)$			$1$	$-0.365$	$-0.423$	$-0.403$
3	$4.078\left(5.339\right)$				$1$	$0.991$	$0.974$
4	$1.825\left(2.345\right)$					$1$	$0.994$
5	$0.506\left(675\right)$						$1$
$a^2\left(Q^2+m_{\eta }^2\right)G_P^{u+d}\left(Q^2\right)$
k	$a_k$	Correlation matrix
0	$0.690\left(128\right)$	$1$	$-0.948$	$0.815$	$0.305$	$0.074$	$-0.022$
1	$-2.805\left(1.929\right)$		$1$	$-0.941$	$-0.194$	$0.071$	$0.174$
2	$3.790\left(6.441\right)$			$1$	$0.200$	$-0.083$	$-0.192$
3	$0.829\left(3.428\right)$				$1$	$0.957$	$0.909$
4	$0.045\left(1.336\right)$					$1$	$0.990$
5	$-0.031\left(360\right)$						$1$
$a^2\left(Q^2+m_{\eta }^2\right)G_P^s\left(Q^2\right)$
k	$a_k$	Correlation matrix
0	$-0.121\left(17\right)$	$1$	$-0.868$	$0.664$	$0.719$	$0.714$	$0.689$
1	$0.256\left(223\right)$		$1$	$-0.935$	$-0.909$	$-0.838$	$-0.765$
2	$0.538\left(766\right)$			$1$	$0.949$	$0.858$	$0.771$
3	$0.191\left(293\right)$				$1$	$0.976$	$0.931$
4	$0.049\left(83\right)$					$1$	$0.988$
5	$0.011\left(20\right)$						$1$

Equations

$$⟨p^{\prime },{\lambda }^{\prime }|A_{\mu }^q|p,\lambda ⟩=\bar{u}(p^{\prime },{\lambda }^{\prime })\left[{\gamma }_{\mu }{G}_A^q(Q^2)+\frac{(p^{\prime }-p{)}_{\mu }}{2m_N}{G}_P^q(Q^2)\right]{\gamma }_{5}u(p,\lambda ),$$

$$\Delta q\equiv g_A^q\equiv G_A^q(0)={\int }_0^{1}dx\left(\Delta q(x)+\Delta \bar{q}(x)\right).$$

$$G_A^q(Q^2)=g_A^q\left(1-\frac{1}{6}(r_A^2{)}^q{Q}^2+O(Q^4)\right).$$

$$C_{\text{2pt}}(\vec{p},t)=\sum _{\vec{x}}{e}^{-i\vec{p}\cdot \vec{x}}\mathrm{Tr}\left[{\Gamma }_{\text{pol}}⟨\chi (\vec{x},t)\bar{\chi }(\vec{0},0)⟩\right]$$

$$C_{\text{3pt}}^{A_{\mu }^q}(\vec{p},\vec{p}{}^{\prime },\tau ,T)=\sum _{\vec{x},\vec{y}}{e}^{-i\vec{p}{}^{\prime }\cdot \vec{x}}{e}^{i(\vec{p}{}^{\prime }-\vec{p})\cdot y}\mathrm{Tr}\left[{\Gamma }_{\text{pol}}⟨\chi (\vec{x},T)A_{\mu }^q(\vec{y},\tau )\bar{\chi }(\vec{0},0)⟩\right],$$

$$R_{\mu }^q(\vec{p},\vec{p}{}^{\prime },\tau ,T)$$

$$=M_{\mu }^q(\vec{p},\vec{p}{}^{\prime })+O(e^{-\Delta E_{10}(\vec{p})\tau })+O(e^{-\Delta E_{10}(\vec{p}{}^{\prime })(T-\tau )}),$$

$$S_{\mu }^q(\vec{p},\vec{p}{}^{\prime },T)\equiv a\sum _{\tau /a=1}^{T/a-1}{R}_{\mu }^q(\vec{p},\vec{p}{}^{\prime },\tau ,T)=c+T{M}_{\mu }^q(\vec{p},\vec{p}{}^{\prime })+O(T{e}^{-\Delta E_{\text{min}}T}).$$

$$T_{\mu }^q(\vec{k},t)\equiv -\sum _{\vec{x}}{e}^{i\vec{k}\cdot \vec{x}}\mathrm{Tr}\left[{\gamma }_{\mu }{\gamma }_5{D}_q^{-1}(x,x)\right],$$

$${\eta }_{a\alpha }^{[d,n]}(\vec{x})=\sum _{b,\beta }(P_d{)}_{a\alpha }^{b\beta }{z}_n(\vec{x}){\eta }_{b\beta }(\vec{x}).$$

$$\frac{-1}{N_{\text{hvec}}}\sum _{d,n}\sum _{\vec{x}}{e}^{i\vec{k}\cdot \vec{x}}{\eta }^{[d,n]†}(\vec{x},t){\gamma }_{\mu }{\gamma }_5{\psi }_q^{[d,n]}(\vec{x},t).$$

$$(A_{\mu }^q-A_{\mu }^{q^{\prime }}{)}^{R,I}$$

$${\left(\sum _q{A}_{\mu }^q\right)}^{R,I}$$

$$z(Q^2)=\frac{\sqrt{t_{\text{cut}}+Q^2}-\sqrt{t_{\text{cut}}}}{\sqrt{t_{\text{cut}}+Q^2}+\sqrt{t_{\text{cut}}}},$$

$$G(Q^2)=\{\begin{array}{ll}{G}_A(Q^2)& \\ (Q^2+m_{\pi }^2)G_P(Q^2)& \text{isovector}\\ (Q^2+m_{\eta }^2)G_P(Q^2)& \text{isoscalar}\end{array}.$$

$$G(Q^2)=\sum _{k=0}^{k_{\text{max}}}{a}_{k}z(Q^2{)}^k.$$

$${\chi }_{\text{aug}}^2\equiv \sum _{i,j}\left(G(Q_i^2)-\sum _k{a}_{k}z(Q_i^2{)}^k\right){\Xi }_{ij}\left(G(Q_j^2)-\sum _{k^{\prime }}{a}_{k^{\prime }}z(Q_j^2{)}^{k^{\prime }}\right)+\sum _{k>1}\frac{a_k^2}{w^2}$$

$$R_{ij}\equiv \frac{C_{ij}}{\sqrt{C_{ii}{C}_{jj}}},$$

$$R_{ij}^{\star }=\{\begin{array}{ll}1& i=j\\ (1-{\lambda }_1)R_{ij}& \text{i and j have different (p′−p)2}\\ (1-{\lambda }_2)R_{ij}+{\lambda }_{2}r& \text{i and j have the same (p′−p)2}\end{array}.$$

$$S(p)=\frac{1}{V}\sum _{x,y}{e}^{-ip(x-y)}⟨q(x)\bar{q}(y)⟩,$$

$$G_{\mathcal{O}}(p,p^{\prime })=\frac{1}{V}\sum _{x,y,z}{e}^{-i{p}^{\prime }.(x-y)-ip.(y-z)}⟨q(x)\mathcal{O}(y)\bar{q}(z)⟩,$$

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

$$\begin{array}{c}\underset{m\to 0}{\lim }\frac{-i}{12p^2}\mathrm{Tr}\left[S_R^{-1}(p)/p\right]=1,\\ \underset{m\to 0}{\lim }\frac{1}{36}\mathrm{Tr}\left[{\Lambda }_{V_{\mu }}^R(p,p)\left({\gamma }_{\mu }-\frac{p_{\mu }/p}{p^2}\right)\right]=1,\\ \underset{m\to 0}{\lim }\frac{1}{36}\mathrm{Tr}\left[{\Lambda }_{A_{\mu }}^R(p,p){\gamma }_5\left({\gamma }_{\mu }-\frac{p_{\mu }/p}{p^2}\right)\right]=1.\end{array}$$

$$\begin{array}{c}\underset{m\to 0}{\lim }\frac{1}{12q^2}\mathrm{Tr}\left[q_{\mu }{\Lambda }_{V_{\mu }}^R(p,p^{\prime })/q\right]=1,\\ \underset{m\to 0}{\lim }\frac{1}{12q^2}\mathrm{Tr}\left[q_{\mu }{\Lambda }_{A_{\mu }}^R(p,p^{\prime }){\gamma }_5/q\right]=1.\end{array}$$

$$(Z_q{)}_{\text{RI′-MOM}}=\frac{Z_V}{36}\mathrm{Tr}\left[{\Lambda }_{V_{\mu }}(p,p)\left({\gamma }_{\mu }-\frac{p_{\mu }/p}{p^2}\right)\right].$$

$${\Lambda }_{\text{QCD}}\ll \mu \ll \pi /a.$$

$$Z_V\left(R(t_1)-R(t_2)\right)=1,$$

$$Z_A(p^2,p^{[4]},p^{[6]},\dots )=Z_A^0(p^2)+c_1\frac{a^2{p}^{[4]}}{p^2}+c_2{\left(\frac{a^2{p}^{[4]}}{p^2}\right)}^2+c_3\frac{a^4{p}^{[6]}}{p^2}+c_4{a}^4{p}^{[4]}.$$

$$A_{\mu }^0=\frac{1}{\sqrt{N_f}}\bar{\psi }{\gamma }_{\mu }{\gamma }_5\psi ,A_{\mu }^a=\bar{\psi }{\gamma }_{\mu }{\gamma }_5{\lambda }^a\psi ,$$

$${\mu }^2\frac{d}{d{\mu }^2}\log \left(Z_5^s{Z}_A^{\overline{\text{MS}},\text{HV}}\right)$$

$${\mu }^2\frac{d}{d{\mu }^2}\alpha$$

$$\begin{array}{c}\begin{array}{rl}{\beta }_0& =\frac{1}{4\pi }\left(\frac{11}{3}{C}_A-\frac{4}{3}{T}_F{N}_f\right)=\frac{1}{4\pi }\left(11-\frac{2}{3}{N}_f\right),\\ {\beta }_1& =\frac{1}{(4\pi {)}^2}\left(\frac{34}{3}{C}_A^2-\frac{20}{3}{C}_A{T}_F{N}_f-4C_F{T}_F{N}_f\right)=\frac{1}{(4\pi {)}^2}\left(102-\frac{38}{3}{N}_f\right),\\ {\gamma }_0& =0,\\ {\gamma }_1& =\frac{1}{(4\pi {)}^2}(6C_F{N}_f)=\frac{1}{(4\pi {)}^2}8N_f,\end{array}\end{array}$$

Full text

I Introduction

The axial and induced pseudoscalar form factors111We also denote flavor combinations using, e.g., $G_{A}^{u-d}(Q^{2})\equiv G_{A}^{u}(Q^{2})-G_{A}^{d}(Q^{2})$., $G_{A}^{q}(Q^{2})$ and $G_{P}^{q}(Q^{2})$, parameterize matrix elements of the axial current between proton states:

[TABLE]

where $Q^{2}=-(p^{\prime}-p)^{2}$ and $A_{\mu}^{q}=\bar{q}\gamma_{\mu}\gamma_{5}q$. It has been shown that $G_{A}^{q}(Q^{2})$ can be interpreted as the two-dimensional Fourier transform of the difference $q_{\uparrow}(\mathbf{b}_{\perp})+\bar{q}_{\uparrow}(\mathbf{b}_{\perp})-q_{\downarrow}(\mathbf{b}_{\perp})-\bar{q}_{\downarrow}(\mathbf{b}_{\perp})$ between transverse densities of helicity aligned and anti-aligned quarks plus antiquarks in a longitudinally polarized nucleon, in the infinite momentum frame Burkardt:2002hr .

At $Q^{2}=0$, the axial form factor gives the fractional contribution from the spin of quarks $q$ and $\bar{q}$ to the proton’s spin, which can also be obtained from a moment of polarized parton distribution functions:

[TABLE]

Understanding the constituents of the proton’s spin has been of great interest ever since the European Muon Collaboration found, by measuring the spin asymmetry in polarized deep inelastic scattering, that the total contribution from quark spin to the proton’s spin is less than half Ashman:1987hv .

Axial form factors naturally arise in the interactions of nucleons with $W$ and $Z$ bosons. Assuming isospin symmetry, the $W$ boson is sensitive to the $u-d$ flavor combination, whereas the $Z$ boson is also sensitive to strange quarks. Neutron beta decay, mediated by $W$-boson exchange, is used to determine the “axial charge” $g_{A}\equiv g_{A}^{u-d}$. Quasielastic neutrino scattering, $\nu n\to\ell^{-}p$ or $\bar{\nu}p\to\ell^{+}n$, has been used to measure the isovector axial form factor $G_{A}^{u-d}(Q^{2})$, whereas elastic neutrino scattering is also sensitive to $G_{A}^{s}(Q^{2})$. The shape of the isovector axial form factor is often assumed to be a dipole, $G_{A}^{u-d}(Q^{2})=g_{A}/(1+Q^{2}/m_{A}^{2})^{2}$; rather than assume a dipole, we will use a more general fit and characterize the shape using the squared axial radii $(r_{A}^{2})^{q}$. These are defined from the slope of the form factors at zero momentum transfer222In contrast with the strange magnetic radius $(r_{M}^{2})^{s}\equiv-6\frac{d}{dQ^{2}}G_{M}^{s}(Q^{2})|_{Q^{2}=0}$, we choose to normalize the strange axial radius relative to the value of the form factor at $Q^{2}=0$, the same as for all the axial radii. Note that this means the flavor combinations satisfy, e.g., $g_{A}^{u-d}(r_{A}^{2})^{u-d}=g_{A}^{u}(r_{A}^{2})^{u}-g_{A}^{d}(r_{A}^{2})^{d}$.:

[TABLE]

The ordinary “axial radius” is the isovector one, $r_{A}\equiv\sqrt{(r_{A}^{2})^{u-d}}$; in the dipole model, $r_{A}^{2}=12/m_{A}^{2}$. It can also be determined from pion electroproduction, using chiral perturbation theory Bernard:2001rs .

In addition to the valence up and down quarks, quantum fluctuations cause other quarks to play a role in the structure of nucleons; the strange quark is the next lightest, and is expected to be the next most important. In this paper, we report a calculation of the nucleon axial form factors using a single lattice QCD ensemble. This calculation includes both quark-connected and disconnected diagrams, which allows us to determine the up, down, and strange form factors. Using the same dataset, we previously reported a high-precision calculation of the strange nucleon electromagnetic form factors Green:2015wqa .

A lattice QCD study of the axial form factors of the nucleon is timely not least in view of experimental efforts underway using the MicroBooNE liquid Argon time-projection chamber, which, in particular, will be able to map out the strange axial form factor of the nucleon to momentum transfers as low as $Q^{2}=0.08\text{ GeV}^{2}$ Papavassiliou:2009zz . This is achieved by combining neutrino-proton neutral and charged current scattering cross section measurements with available polarized electron-proton/deuterium cross section data, and is expected to reduce the experimental uncertainty of the extrapolated value at $Q^{2}=0$, i.e., the strange quark spin contribution $\Delta s$, by an order of magnitude. Such an extraction is complementary to polarized DIS determinations that access the strange quark helicity distribution function, but suffer from lack of coverage at low and high momentum fraction $x$ when evaluating the first $x$-moment. The $Q^{2}$ range explored by the MicroBooNE experiment, between $Q^{2}=0.08\text{ GeV}^{2}$ and about $Q^{2}=1\text{ GeV}^{2}$, matches the range covered by the present lattice calculation well, enabling a future comparison of the $Q^{2}$-dependence obtained for the strange axial form factor.

This paper is organized as follows. Section II describes our methodology: the approach used to isolate the nucleon ground state and determine the form factors, the methods used to determine the numerically-challenging disconnected diagrams, the details of the lattice ensemble, and the fits to the $Q^{2}$-dependence of the form factors using the $z$ expansion. The unwanted contributions from excited states to the different observables are examined in detail, and the estimation of systematic uncertainty is described. Our nonperturbative calculation of the renormalization factors, including a nonperturbative treatment of the flavor singlet case, is presented in Sec. III. The main results are in Sec. IV: the axial and induced pseudoscalar form factors for light and strange quarks, as well as the quark contributions to the nucleon spin. Finally, we present our conclusions in Sec. V. In an appendix, we give the parameters for our fits to the form factors.

II Lattice methodology

II.1 Computation of matrix elements

To determine nucleon matrix elements, we compute two-point and three-point functions,

[TABLE]

where $\chi=\epsilon^{abc}(\tilde{u}^{T}_{a}C\gamma_{5}\frac{1+\gamma_{4}}{2}\tilde{d}_{b})\tilde{u}_{c}$ is a proton interpolating operator and $\Gamma_{\text{pol}}$ is a spin and parity projection matrix. In the interpolating operator, we use Wuppertal-smeared Gusken:1989qx quark fields $\tilde{q}=(\frac{1+\alpha H}{1+6\alpha})^{N}q$, where $H$ is the nearest-neighbor gauge-covariant hopping matrix constructed using spatially APE-smeared Albanese:1987ds gauge links.

The proton ground state can be obtained in the limit where all time separations $t$, $\tau$, and $T-\tau$ are large. In this limit, the following ratio does not depend on the time separations or on the interpolating operator:

[TABLE]

where $M_{\mu}^{q}(\vec{p},\vec{p}\,^{\prime})$ contains the desired nucleon matrix element $\langle\vec{p}\,^{\prime},\lambda^{\prime}|A_{\mu}^{q}|\vec{p},\lambda\rangle$ (with spins depending on $\Gamma_{\text{pol}}$) and some kinematic factors (see, e.g., Green:2014xba ), and $\Delta E(\vec{p})$ is the energy gap between the ground and lowest excited state with momentum $\vec{p}$.

For each source-sink separation $T$, for the ratio-plateau method, we take the average of the central two or three points $R_{\mu}^{q}(\vec{p},\vec{p}\,^{\prime},\tau,T)$ near $\tau=T/2$. This gives an estimate of $M_{\mu}^{q}(\vec{p},\vec{p}\,^{\prime})$ (and thus the nucleon matrix element) with a systematic error coming from excited-state contamination that decays exponentially as $e^{-\Delta E_{\text{min}}T/2}$, where $\Delta E_{\text{min}}=\min\{\Delta E_{10}(\vec{p}),\Delta E_{10}(\vec{p}\,^{\prime})\}$. We also use the summation method, computing the sums

[TABLE]

Fitting the slope with respect to $T$ yields an estimate of $M_{\mu}^{q}(\vec{p},\vec{p}\,^{\prime})$ that has a greater suppression of unwanted excited-state contributions Capitani:2010sg ; Bulava:2010ej , which now decay as $Te^{-\Delta E_{\text{min}}T}$.

For each $Q^{2}$, we construct a system of equations parameterizing the corresponding set of matrix elements of the axial current with $G_{A}(Q^{2})$ and $G_{P}(Q^{2})$. We combine equivalent matrix elements to improve the condition number Syritsyn:2009mx , and then solve the resulting overdetermined system of equations Hagler:2003jd . This approach makes use of all available data to minimize the statistical uncertainty. In particular, for disconnected diagrams, we are able to compute correlators for all polarizations and all equivalent momenta, maximizing the amount of averaging.

II.2 Disconnected diagrams

There are two kinds of quark contractions that contribute to $C_{\text{3pt}}$: connected and disconnected, shown in Fig. 1. We evaluate the former exactly for each source on each gauge configuration, using sequential propagators through the sink Martinelli:1988rr . For the latter, we perform a stochastic evaluation of the disconnected loop,

[TABLE]

where $D_{q}$ is the lattice Dirac operator with a fixed gauge background and $x=(\vec{x},t)$. We then obtain the disconnected contribution to $C_{\text{3pt}}$ from the correlation between this loop and the nucleon two-point function.

To evaluate the disconnected loop, we generate noise fields $\eta_{a\alpha}(x)$ that have color, spin, and space-time indices but with support only on a single timeslice333In this work we have not compared the effectiveness of placing noise on one timeslice against placing it on all timeslices., $t$. We use one $\mathbb{Z}_{2}+i\mathbb{Z}_{2}$ noise vector for each chosen timeslice and gauge configuration, i.e., the components of $\eta$ are randomly chosen from $\{\frac{1+i}{\sqrt{2}},\frac{1-i}{\sqrt{2}},\frac{-1+i}{\sqrt{2}},\frac{-1-i}{\sqrt{2}}\}$. As a result, the diagonal elements of $\eta\eta^{\dagger}$ are equal to 1, and the off-diagonal elements are random with expectation value zero. To reduce noise by replacing statistical zeros with exact zeros in targeted off-diagonal components of $\eta\eta^{\dagger}$, we use color and spin dilution Wilcox:1999ab ; Foley:2005ac , as well as hierarchical probing Stathopoulos:2013aci . The former makes use of a complete set of twelve projectors in color and spin space, $P_{d}$, such that $P_{d}\eta$ has support on only one color and one spin component. The latter makes use of $N_{\text{hvec}}$ specially-constructed spatial Hadamard vectors, $z_{n}$, that provide a scheme for progressively eliminating the spatially near-diagonal contributions to the noise. Combining these yields $12N_{\text{hvec}}$ modified noise fields,

[TABLE]

We use these as sources for quark propagators, $\psi^{[d,n]}_{q}=D^{-1}_{q}\eta^{[d,n]}$, and obtain an estimator for $T_{\mu}^{q}(\vec{k},t)$:

[TABLE]

We will separately consider the connected and disconnected contributions to nucleon matrix elements of the light quark axial current. Although the individual contributions are unphysical, they can be understood using partially quenched QCD Bernard:1993sv , by introducing a third degenerate light quark $r$ and a corresponding ghost quark to cancel its fermion determinant in the path integral. The disconnected contribution to a nucleon three-point function with current $A_{\mu}^{u}$ or $A_{\mu}^{d}$ is equal to a nucleon three-point function with $A_{\mu}^{r}$. Since it was shown in Ref. Bernard:2013kwa that partially quenched staggered fermions have a bounded transfer matrix, we expect that for our case as well we can separately isolate the ground state in the connected and disconnected contributions to three-point functions, i.e., that Eq. (6) applies to $A_{\mu}^{r}$. In Section III we will also discuss renormalization of $A_{\mu}^{r}$.

II.3 Lattice ensemble and calculation setup

We use a single lattice ensemble with a tree-level Symanzik improved gauge action ($\beta=6.1$) and 2+1 flavors of clover-improved Wilson fermions that couple to the gauge links after stout smearing (one step with $\rho=0.125$). The improvement parameters are set to their tadpole-improved tree-level values. The lattice size is $32^{3}\times 96$ and the bare quark masses are $am_{s}=-0.245$ and $am_{ud}=-0.285$.

Based on the $\Upsilon(2S)-\Upsilon(1S)$ energy splitting computed using lattice NRQCD, the lattice spacing is $a=0.11403(77)$ fm. The strange quark mass is close to its physical value: the mass of the unphysical $\eta_{s}$ meson is 672(3)(5) MeV, which is within 5% of its value determined for physical quark masses Dowdall:2011wh . The light quark mass is heavier than physical, producing a pion mass444For the pion and $\eta_{s}$ mass, the second error is from uncertainty in the lattice spacing. of 317(2)(2) MeV. The volume is quite large, such that $m_{\pi}L_{s}\approx 5.9$, and we thus expect finite-volume effects to be highly suppressed.

We performed calculations using 1028 gauge configurations, on each of which we chose six equally-spaced source timeslices. For each source timeslice $t_{0}$, we used two positions $(\vec{x}_{1},t_{0})$ and $(\vec{x}_{2},t_{0})$ as sources for three-point functions. We placed nucleon sinks in both the forward and backward directions on timeslices $t_{0}\pm T$ to double statistics and obtain a total of 24672 samples, and used five source-sink separations $T/a\in\{6,8,10,12,14\}$. We computed disconnected loops on timeslices $t_{0}+\tau$ displaced only in the forward direction from each source timeslice, yielding 6168 timeslice samples; the source-operator separations $\tau$ and number of Hadamard vectors for each flavor are listed in Tab. 1. For each source timeslice, we computed sixteen two-point functions from source positions $(\vec{x}_{i},t_{0})$, $i=1,\dots,16$, yielding 98688 samples for correlating with the disconnected loops. We imposed two constraints on our choice of momenta: $(\vec{p}\,^{\prime}-\vec{p})^{2}\leq 10(\frac{2\pi}{L_{s}})^{2}$ and $(\vec{p})^{2},(\vec{p}\,^{\prime})^{2}\leq 6(\frac{2\pi}{L_{s}})^{2}$. For the connected diagrams we used two sink momenta, $\vec{p}\,^{\prime}=\vec{0}$ and $\vec{p}\,^{\prime}=\frac{2\pi}{L_{s}}(-1,0,0)$, and all source momenta compatible with the constraints. For the disconnected diagrams we used all combinations of $\vec{p}$ and $\vec{p}\,^{\prime}$ compatible with the constraints, with the restriction that each $Q^{2}$ must match a value available from the connected diagrams.

On each set of four adjacent gauge configurations, we averaged over all spatially displaced samples of each correlator. This produced 257 blocked samples. Statistical error analysis was done using jackknife resampling.

The general form for $O(a)$ improvement of quark bilinear operators with nondegenerate quarks was given in Ref. Bhattacharya:2005rb . If we simplify the expressions by keeping only their form at one-loop order in perturbation theory, the renormalized improved operators take the form

[TABLE]

for the flavor nonsinglet and singlet cases, respectively, where $P$ is the pseudoscalar density. Matching with the improvement of the action, we take the tree-level value $c_{A}=0$. Note that in nucleon matrix elements, the term proportional to $c_{A}$ only contributes to the $G_{P}$ form factors and therefore this term is not necessary for $O(a)$ improvement of $G_{A}(Q^{2})$.555In practice lattice results for $G_{A}(Q^{2})$ could depend on $c_{A}$ indirectly due to contamination from excited states, or from a breakdown of the form factor decomposition (1) due to breaking of rotational symmetry. The latter can result from either the UV cutoff (an $O(a^{2})$ effect) or the IR cutoff (suppressed by $e^{-m_{\pi}L}$). The mass-dependent terms can effectively cause a mixing between singlet and nonsinglet axial currents; rather than determine $b_{A}$ explicitly, we absorb the mass-dependent terms into the renormalization factors, which now become a matrix. The renormalization matrix is determined nonperturbatively using the Rome-Southampton method, which we discuss in detail in Section III.

II.4 Effectiveness of hierarchical probing

On a reduced set of 366 configurations, we have data for the disconnected light-quark loops from two different methods: hierarchical probing, as used for the main calculations of this work, and “Noise only”, where the sum over $n$ in Eq. (10) is over $N$ random noise samples rather than $N_{\rm hvec}$ Hadamard vectors multiplying a single noise sample. Note that this means color and spin dilution is used in both cases. Thus, at $N=N_{\rm hvec}$ the computational cost for both methods is the same. Figure 2 shows results from both methods as a function of $N=N_{\rm hvec}$. Hierarchical probing is always guaranteed to perform at least as well as the traditional noise method. For our setup we find that the uncertainty in the disconnected light-quark $g_{A}$ saturates at $N_{\rm hvec}=16$, where it becomes dominated by gauge noise. For $g_{A}$ with $N=N_{\rm hvec}=16$, the reduction in the (combined gauge+stochastic) uncertainty is only by a modest factor of 1.4. The improvement from hierarchical probing is more significant for the disconnected electromagnetic form factors Green:2015wqa , as illustrated in Fig. 2 (right) for the disconnected light-quark contribution to $G_{M}$ at $Q^{2}\approx 0.11\>\>{\rm GeV}^{2}$. In this case, the stochastic noise dominates over the gauge noise up to a larger value of $N$ (saturation is not yet reached in the range considered), and at large $N$ the improvement from hierarchical probing is more pronounced, as expected because of the greater “coloring distance” Stathopoulos:2013aci .

II.5 Excited-state effects

It turns out that the different form factors suffer from quite different amounts of excited-state contamination. In addition, the available $(T,\tau)$ combinations are quite different between our connected-diagrams data and our disconnected-diagrams data. In particular, the former are much better suited for applying the summation method than the latter. Therefore we choose the best method for isolating the ground state separately for each form factor. We do this by examining “plateau” plots where, for each $(T,\tau)$ we determine “effective” form factors666In this subsection we show bare form factors, i.e. before renormalization. from the ratios assuming the absence of excited states. In a region where excited-state effects are negligible, these effective form factors will form a stable plateau. In addition to these plateaus from the ratio method, we also show results from the summation method, taking the sums with three adjacent points $\{T,T+2a,T+4a\}$ and fitting with a line to determine the slope.

Figure 3 (top row) shows plateau plots for the isovector axial form factor $G_{A}^{u-d}(Q^{2})$. For the axial charge $g_{A}\equiv G_{A}^{u-d}(0)$ (top left), the centers of the plateaus appear stable by $T/a=10$ and 12, which agree within uncertainty. The center of the plateau for the largest source-sink separation, $T=14a$, is shifted significantly higher, however its statistical uncertainty is quite large and the magnitude of the shift goes against expectations: in the asymptotic regime, as $T$ is increased the shift between neighboring values of $T$ is expected to decrease. Therefore we conclude that the shift at $T=14a$ is likely a statistical fluctuation777Similar behavior was previously seen in the isovector Pauli form factor computed using the same dataset Green:2013hja . and take the results from $T=12a$ as the best option using the ratio method. For the summation method, all three points are consistent within the uncertainty and we conclude that the summation method has reached a plateau already at the shortest source-sink separation, $T=6a$ (i.e., from fitting to the sums with $T/a\in\{6,8,10\}$). We take this as our primary analysis method for the isovector axial form factor $G_{A}^{u-d}(Q^{2})$. For this form factor and for any observable derived from it, we estimate systematic uncertainty due to excited-state effects as the root-mean-square (RMS) deviation between the primary result (summation with $T=6a$) and two alternatives: the ratio method with $T=12a$ and the summation method with $T=8a$. Looking at the corresponding plateau plot (top right) for the isovector axial form factor at our largest momentum transfer (about 1.1 GeV2) indicates that this approach is also reasonable at nonzero $Q^{2}$. The bottom row of the same figure shows the equivalent plots for the contribution from quark-connected diagrams to the isoscalar axial form factor $G_{A}^{u+d}(Q^{2})$. The excited-state effects appear to be slightly milder than for the isovector case, and we thus choose to apply the same analysis strategy.

Plateau plots for the contributions from quark-disconnected diagrams to axial form factors are shown in Fig. 4. Note that since these form factors were computed for several fixed source-operator separations $\tau$, we choose to use the operator-sink separation $T-\tau$ as the horizontal axis. The top row shows the light-quark case, where we computed disconnected loops for five source-operator separations, and the bottom row shows the strange-quark case where we only computed three source-operator separations. The left and right columns show $Q^{2}=0$ (i.e., the contributions to the nucleon spin) and our largest momentum transfer, respectively. In general, we do not see any significant dependence on $T-\tau$ for $T-\tau\gtrsim 5a$. Since the disconnected data were averaged over the exchange of source and sink momenta, the effective form factors are expected to be symmetric, and therefore this corresponds to a source-sink separation of $T=10a$. We use this for our primary result (averaged over the three points near $\tau=T/2$, which reduces statistical uncertainty), and use the RMS deviation with results from $T=8a$ and $T=12a$ as our estimate of systematic uncertainty due to excited states.

The isovector induced pseudoscalar form factor $G_{P}^{u-d}(Q^{2})$ at the lowest available momentum transfer (about 0.1 GeV2) is shown in Fig. 5 (left). This has very large excited-state effects (there is nearly a factor of two between the smallest and largest value on the plot), and there is no sign that a plateau has been reached using the ratio method. For the summation method, the points with $T/a=8$ and 10 are consistent, suggesting that a plateau might possibly have been reached. We take the summation method with $T=8a$ as our primary analysis method for this form factor and estimate the systematic uncertainty as the RMS deviation between the primary result and those from the ratio method with $T/a=14$ and 12. Although the latter is clearly not in the plateau regime, we nevertheless include it in order to reflect the poor control over excited-state effects that is available in our data. At larger $Q^{2}$ (right), the excited-state effects are much milder and our error estimate should be conservative.

Plateau plots for the light and strange isoscalar induced pseudoscalar form factors are shown in Fig. 6. For $G_{P}^{u+d}(Q^{2})$ at the lowest available momentum transfer (top left), we again find that the connected contributions have significant excited-state effects. On the same plot, we show the partial plateaus (limited to the available values of $\tau$) for the contributions from disconnected diagrams. Although they are a bit noisier, they also appear to contain excited-state effects, with the opposite sign. In fact, the opposite signs cause the sum of connected and disconnected diagrams to have smaller excited-state contamination. For the sum, using the ratio method with $T=10a$ appears to be a safe choice, also at the maximum momentum transfer (top right). When we examine the individual connected and disconnected contributions, we will make the same choice, with the understanding that the results include some contamination from excited states, and can only be studied qualitatively. This choice also appears safe for $G_{P}^{s}(Q^{2})$ (bottom left and right). As for the disconnected $G_{A}$ form factors, we use the RMS difference with $T/a=8$ and 12 as our estimate of systematic uncertainty due to excited states.

II.6 Form factor fits using the $z$ expansion

Having computed nucleon form factors at several discrete values of $Q^{2}$, we fit them with curves to characterize their overall shape and determine observables such as the axial radius from their slope at $Q^{2}=0$. It has been common to perform these fits using simple ansatzes, such as a dipole, which is often used to describe experimental data for the isovector $G_{A}(Q^{2})$, however these tend to be highly constrained and introduce a model dependence into the results.

Instead, we use the model-independent $z$ expansion. This was used in Refs. Bhattacharya:2011ah ; Bhattacharya:2015mpa ; Meyer:2016oeg to study axial form factors determined from quasielastic (anti)neutrino-nucleon scattering; it was found that fitting with the $z$ expansion produced a significantly larger axial radius with a larger uncertainty, compared with dipole fits. The $z$ expansion makes use of a conformal mapping from $Q^{2}$, where the given form factor is analytic on the complex plane outside a branch cut on the timelike real axis, to the variable $z$ such that the form factor is analytic for $|z|<1$. We use

[TABLE]

where we use the particle production threshold for the isovector form factors, $t_{\text{cut}}=(3m_{\pi})^{2}$. For the isoscalar form factors the actual threshold may be higher, but we use the same $t_{\text{cut}}$ everywhere for simplicity. We have chosen the mapping such that $Q^{2}=0$ maps to $z=0$.

The $G_{P}$ form factors have an isolated pole below the particle production threshold at the pseudoscalar meson mass, which we remove before fitting. We thus perform fits to

[TABLE]

Each form factor can be described by a convergent Taylor series in $z$. We truncate this series and obtain our fit form,

[TABLE]

The first two coefficients, $a_{0}$ and $a_{1}$, give the intercept and slope of the form factor at $Q^{2}=0$. Specifically, $G(0)=a_{0}$ and, for the axial form factors, $r_{A}^{2}=-3a_{1}/(2a_{0}t_{\text{cut}})$. We impose Gaussian priors on the remaining coefficients, centered at zero with width equal to $w=5\max\{|a_{0}|,|a_{1}|\}$. The series is truncated with $k_{\text{max}}=5$, but this is large enough that increasing it further has no effect in our probed range of $Q^{2}$; i.e., the priors cause $a_{k}z^{k}$ to be negligible for $k>5$.

We perform correlated fits, minimizing

[TABLE]

with respect to $\{a_{k}\}$, where $\Xi$ is an estimator for the inverse covariance matrix and the last term augments the chi-squared with the Gaussian priors. With limited statistics it can be difficult to obtain a reliable estimator, and therefore we choose to reduce statistical fluctuations by interpolating between the jackknife estimate of the covariance matrix and a simplified (less noisy but biased) estimate, and then inverting the resulting matrix. This is in the spirit of shrinkage estimators Ledoit_2004 ; Schaefer_2005 , however we do not perform an optimization step with respect to the interpolation parameter.

In order to choose the form of the target (simplified) covariance matrix, we examine the correlation matrix

[TABLE]

where $C$ is the jackknife estimate of the covariance matrix. We find that this has a quite different form between connected diagrams and disconnected diagrams. Figure 7 shows two example correlation matrices. For connected diagrams, illustrated with $G_{A}^{u-d}(Q^{2})$ (left), we find modest correlations between different values of $Q^{2}$ but no strong pattern. For disconnected diagrams, illustrated with the quark-disconnected contribution to the light-quark $G_{A}$ form factor (right), the correlation matrix is nearly block-diagonal. Each block corresponds to values of $Q^{2}$ that share the same spatial momentum transfer $(\vec{p}\,^{\prime}-\vec{p})^{2}$ and thus the same Fourier modes of the disconnected loops. There are strong correlations within each block but weak correlations between different blocks.

For connected diagrams, we set $\Xi=\left((1-\lambda)C+\lambda C_{\text{diag}}\right)^{-1}$, where $C_{\text{diag}}$ is the diagonal part of the covariance matrix. This is equivalent to multiplying the off-diagonal elements of $C$ by $1-\lambda$. We use the mild value of $\lambda=0.1$ as our main choice. For disconnected diagrams, we compute the average $r$ over all elements of $R_{ij}$ where $i$ and $j$ ($i\neq j$) correspond to the same spatial momentum transfer. We then use for $\Xi$ the inverse of the matrix $R^{\star}_{ij}\sqrt{C_{ii}C_{jj}}$, where

[TABLE]

As our main choice, we use $(\lambda_{1},\lambda_{2})=(1,\tfrac{1}{2})$.

To estimate systematic uncertainty from fitting, we perform several alternative fits. We halve the value of $w$. For connected diagrams, we perform fits with $\lambda=0$ and 1. For disconnected diagrams, we perform fits with $(\lambda_{1},\lambda_{2})=(0,0)$, $(1,0)$, and $(1,1)$. Finally, we take the RMS difference between results from all of the alternative fits as our estimate.

III Renormalization

To compare our results with phenomenology, the lattice axial current needs to be renormalized. We determine the necessary renormalization factors nonperturbatively using the Rome-Southampton approach Martinelli:1994ty . Going beyond the usual computation of the flavor nonsinglet renormalization factor, we also renormalize the flavor singlet axial current nonperturbatively. This requires disconnected quark loops but we are able to reuse the same loops that were computed for nucleon three-point functions. Since we perform these calculations on just one ensemble without taking the chiral limit, we effectively absorb the mass-dependent operator improvement terms into the renormalization (see Subsec. II.3), which requires us to determine a matrix of renormalization factors.

The singlet-nonsinglet difference in axial renormalization factors has been previously studied nonperturbatively by QCDSF Chambers:2014pea at the $SU(3)$ flavor symmetric point, using additional lattice ensembles and the Feynman-Hellmann relation to determine the contributions from disconnected quark loops. For the case of two degenerate quark flavors, nonperturbative results were presented by RQCD at the Lattice 2016 conference Bali:2017jyw , using stochastic estimation for the disconnected loops similarly to this work. The singlet-nonsinglet difference has also been studied at leading (two-loop) order in lattice perturbation theory for a variety of improved Wilson-type actions Skouroupathis:2008mf ; Constantinou:2016ieh .

This section is organized as follows: we present the Rome-Southampton method and the RI*′*-MOM and RI-SMOM schemes for the single-flavor case in Subsec. III.1, determine the light and strange vector current renormalization factors in Subsec. III.2, study discretization effects and breaking of rotational symmetry in Subsec. III.3, and discuss issues of matching to the $\overline{\text{MS}}$ scheme and running of the flavor singlet axial current in Subsec. III.4. Subsections III.5 and III.6 explain our procedure for calculating the $Z_{A}$ renormalization matrix, and finally we give the details of the calculation and its results in Subsec. III.7.

III.1 Rome-Southampton method, RI*′*-MOM, and RI-SMOM

For calculating the axial renormalization constants, we follow the Rome-Southampton approach in both RI*′*-MOM Martinelli:1994ty ; Gockeler:1998ye and RI-SMOM schemes Sturm:2009kb . In Landau gauge, we compute quark propagators

[TABLE]

Green’s functions,

[TABLE]

and amputated Green’s functions,

[TABLE]

The renormalized quantities are defined as $S_{R}(p)=Z_{q}S(p)$ and $\Lambda_{\mathcal{O}}^{R}(p,p^{\prime})=Z_{q}^{-1}Z_{\mathcal{O}}\Lambda_{\mathcal{O}}(p,p^{\prime})$. In RI*′-MOM, renormalization conditions are imposed for $p^{\prime}=p$, at scale $p^{2}=p^{\prime 2}=\mu^{2}$. For the quark field and vector and axial currents888This combination of conditions for $Z_{q}$, $Z_{V}$, and $Z_{A}$ has also been called the MOM scheme or the RI′* scheme. Note that the name RI*′*-MOM has also been used to refer to the combination of this condition for $Z_{q}$ and the original RI-MOM conditions Martinelli:1994ty for $Z_{V}$ and $Z_{A}$, even though this is not compatible with the vector and axial Ward identities.:

[TABLE]

RI-SMOM conditions are imposed at the symmetric point $p^{2}=p^{\prime 2}=q^{2}=\mu^{2}$, where $q=p^{\prime}-p$. The quark-field renormalization is the same as RI*′*-MOM, whereas for the vector and axial currents:

[TABLE]

As stated previously, in our calculations we do not take the chiral limit. We also avoid directly determining the quark-field renormalization. Instead, we impose the above renormalization conditions on the vector current, which gives $Z_{q}/Z_{V}$, and independently obtain $Z_{V}$ from three-point functions of pseudoscalar mesons. Our estimate for $Z_{q}$ in RI*′*-MOM is then obtained using

[TABLE]

In RI-SMOM, we estimate $Z_{q}$ in the same way using Eq. (22).

The renormalization scale $\mu$ should be chosen such that it is much larger than $\Lambda_{\text{QCD}}$, in order to be able to connect the nonperturbative renormalization schemes to $\overline{\text{MS}}$ using perturbation theory (in our case, this is needed for the flavor-singlet axial current), and much smaller than the inverse lattice spacing to avoid large discretization errors:

[TABLE]

As our lattice spacing is fairly coarse, we do not expect to find a stable plateau region in this window. Instead, we will perform fits to remove the leading $O(a^{2}\mu^{2})$ artifacts, and make use of the two different schemes to estimate unaccounted-for systematic uncertainties.

III.2 Vector current renormalization

We obtain the mass-dependent light and strange vector current renormalization factors from matrix elements of pseudoscalar mesons following, e.g., Ref. Durr:2010aw . For $\pi$ and $\eta_{s}$ states, we compute zero-momentum two-point functions $C_{2}(t)$ as well as three-point functions $C_{3}(t)$ with source-sink separation $T=L_{t}/2$ and an operator insertion of the time component of the local (light or strange) vector current at source-operator separation $t$. We form the ratio $R(t)=C_{3}(t)/C_{2}(T)$, so that the charge of the interpolating operator gives the renormalization condition

[TABLE]

for $0<t_{1}<T<t_{2}<L_{t}$. Taking the difference $R(t)-R(t+T)$ results in a large cancellation of correlated statistical uncertainties. Results are shown in Fig. 8. We average over the long plateau, excluding three points at each end, and obtain $Z_{V}^{l}=0.7903(2)$ and $Z_{V}^{s}=0.8337(2)$.

III.3 Study of discretization effects

We perform a dedicated study of discretization effects and breaking of rotational symmetry, for the isovector case in the RI*′*-MOM scheme. Using translation invariance to remove the sum over $y$ in Eqs. (18) and (19), we compute point-source quark propagators from a fixed point $y$, which allows us to efficiently obtain the gauge-averaged quark propagator and Green’s functions for a large set of momenta. Specifically, we save data for all momenta in the inner 1/16 of the lattice Brillouin zone, i.e., with $|p_{\mu}|\leq\frac{\pi}{2a}$. After checking that the breaking of hypercubic symmetry due to the different lattice temporal and spatial extents is negligible, we averaged the estimates for the isovector $Z_{A}$ over all hypercubic equivalent momenta.

Since the lattice breaks rotational symmetry, estimates of $Z_{A}$ will depend not only on $p^{2}$, but also the hypercubic invariants $p^{[2n]}\equiv\sum_{\mu}(p_{\mu})^{2n}$. We make use of the hypercubic fit form from Refs. Boucaud:2003dx ; Boucaud:2005rm to remove the leading terms that break rotational symmetry and collapse the data to a single function of $p^{2}$:

[TABLE]

The fit parameters are the four $c_{i}$ that control breaking of hypercubic symmetry and a separate $Z_{A}^{0}(p^{2})$ for each $p^{2}$. The data $Z_{A}(p^{2},p^{[4]},\dots)$ and the fit result $Z_{A}^{0}(p^{2})$ are shown in Fig. 9. This is effective at producing a smooth curve that depends only on $p^{2}$ and not the other hypercubic invariants. The resulting curve still contains $O(a^{2}p^{2})$ rotationally invariant lattice artifacts, so we perform a second fit in the range $a^{2}p^{2}\in[2,6]$ assuming a quadratic dependence on $a^{2}p^{2}$, and extrapolate to $a^{2}p^{2}=0$; this is also shown in Fig. 9.

An alternative approach is to pick an initial direction $p_{*}$ and restrict our analysis to points $p=\lambda p_{*}$. Then the hypercubic invariants have the form $p^{[2n]}=c_{2n}p^{2n}$ for some fixed $c_{2n}$ that depend on $p_{*}$. Thus, for this set of points along a fixed direction, the dependence on hypercubic invariants reduces to dependence only on $p^{2}$. We choose four sets of points: on-axis momenta, and momenta along 2, 3, or 4-dimensional diagonals, i.e., $p_{*}=(0,0,0,1)$, $(0,0,1,1)$, $(0,1,1,1)$, and $(1,1,1,1)$. For each set of points, we again do a fit to extrapolate $a^{2}p^{2}$ to zero. Because in this case there are fewer points available, we expand the fit range to be $a^{2}p^{2}\in[1.5,10]$. For on-axis points we use a linear fit because $a^{2}p^{2}$ does not reach very high, and for the $n$-dimensional diagonals we use a quadratic fit. The points from each set and the fit curves are shown in Fig. 9.

We find that the $Z_{A}$ determined from the hypercubic fit and from the fits along different diagonals are all consistent with one another. This indicates that we can reliably control these lattice artifacts by choosing only points along a fixed direction, which is the approach that we will use for our main results for the axial renormalization matrix.

III.4 Matching to $\overline{\text{MS}}$ and running of the singlet axial current

We consider the singlet and nonsinglet axial currents,

[TABLE]

where $\psi$ is the fermionic field and $\lambda^{a}$ is an $SU(N_{f})$ generator acting in flavor space. The nonsinglet current should be renormalized such that it satisfies the axial Ward identity associated with chiral symmetry, and the renormalized singlet current should satisfy the one-loop form of the axial anomaly. The nonsinglet axial current has no anomalous dimension and is appropriately renormalized to all orders in perturbation theory in $\overline{\text{MS}}$ (using dimensional regularization with a naive anticommuting version of $\gamma_{5}$), RI*′*-MOM and RI-SMOM schemes. Thus the matching factor between these schemes is 1, and $Z_{A}=1$ when using a chiral regulator.

For the singlet current, dimensional regularization with a naive $\gamma_{5}$ is inappropriate since the anomaly is not reproduced, and thus the ’t Hooft-Veltman prescription for $\gamma_{5}$ is necessary. Using it in $\overline{\text{MS}}$, an additional finite matching factor $Z_{5}^{s}$ is needed for the renormalized current to satisfy the one-loop form of the axial anomaly Larin:1993tq . Thus renormalized, the singlet current has an anomalous dimension, $\gamma=(\frac{\alpha}{4\pi})^{2}(-6C_{F}N_{f})+O(\alpha^{3})$ Kodaira:1979pa , where the $O(\alpha^{3})$ term is given in Ref. Larin:1993tq . Using the same dimensional regularization, it was shown in Ref. Bhattacharya:2015rsa that the conversion factor between $\overline{\text{MS}}$ (including the finite factor $Z_{5}^{s}$) and RI-SMOM is $1+O(\alpha^{2})$.

For computing the matching between RI*′-MOM and RI-SMOM, at one-loop order there should be no distinction between singlet and nonsinglet currents. Since the matching factor is $1$ for nonsinglet currents, we conclude that the conversion factor for the singlet axial current in RI′*-MOM is $1+O(\alpha^{2})$.

We remove the running of the singlet $Z_{A}$ by evolving to a fixed scale. The evolution is given by

[TABLE]

where the relevant coefficients are

[TABLE]

using $C_{A}=3$, $C_{F}=4/3$, and $T_{F}=1/2$. At two-loop order, the evolution of $\alpha$ is given by Prosperi:2006hx :

[TABLE]

where $W_{k}$ is the many-valued Lambert function defined by $W_{k}(\zeta)e^{W_{k}(\zeta)}=\zeta$. We use the PDG value, $\Lambda_{3}^{\overline{\text{MS}}}=332(19)$ MeV Olive:2016xmw . Using $\gamma_{0}=0$, the evolution of the renormalization factor at two-loop order is given by

[TABLE]

III.5 Renormalization of the axial current: $N_{f}=2+1$

Consider the flavor-diagonal axial currents, Eq. (27), with $\psi=(u\;d\;s)^{T}$. We take $a=3,8,0,$ with $\operatorname{Tr}(\lambda^{a}\lambda^{b})=\delta^{ab}$,

[TABLE]

Using $i,j$ to label quark flavors, we compute the quark propagator $S_{i}(p)$ [Eq. (18)] for quark flavor-$i$, nonamputated and amputated Green’s functions [Eq. (19), Eq. (20)] for mixed quark flavors-$i$ and -$j$, $G_{i,j}^{\mathcal{O}}(p^{\prime},p)$, and $\Lambda_{ij}^{\mathcal{O}}(p^{\prime},p)$, respectively. These renormalize as

[TABLE]

For $N_{f}=2+1$, the renormalization pattern is

[TABLE]

and for $N_{f}=3$, this reduces to two independent factors since $Z_{A}^{88}=Z_{A}^{33}$ and $Z_{A}^{80}=Z_{A}^{08}=0$.

In a RI*′*-MOM or RI-SMOM scheme, the renormalization condition for $Z_{A}$ involves tracing $\Lambda^{A_{\mu}}$ with some projector $P_{\mu}$ at kinematics corresponding to the scale $\mu^{2}$ (see Subsec. III.1). In the case of multiple flavors, this becomes

[TABLE]

so that we get

[TABLE]

Specifically, this yields for $N_{f}=2+1$

[TABLE]

where $\Sigma_{l,\text{conn}}$ is the connected contribution to the ($u$ or $d$)-quark amputated axial vertex function, traced with $P_{\mu}$. This corresponds to the usual isovector result. Writing $\Sigma_{j,\text{disc}}^{i}$ for the disconnected contribution to the amputated vertex function with the flavor-$i$ axial current and flavor-$j$ external quark states, traced with $P_{\mu}$, we get

[TABLE]

It is clear that $(Z_{A}^{-1})^{80}$ and $(Z_{A}^{-1})^{08}$ vanish when $N_{f}=3$, and the disconnected contribution to $(Z_{A}^{-1})^{88}$ is doubly suppressed by approximate $SU(3)_{f}$ symmetry.

Having evaluated an effective $Z_{A}^{-1}$ in some scheme at a scale $\mu$, we can invert the matrix and evolve to the target scale of 2 GeV:

[TABLE]

where the perturbative flavor-singlet evolution is given by Eq. (32). Finally, we fit with a polynomial in $a^{2}\mu^{2}$ to remove lattice artifacts. If we want to obtain a single-flavor axial current, such as the strange, we can write, e.g.,

[TABLE]

Similarly, we can evaluate the renormalized $u+d$ current,

[TABLE]

In order to study the disconnected light-quark current by itself, as described in Subsection II.2, we introduce a quenched third light quark $r$, degenerate with $u$ and $d$. Then the connected contribution to the matrix elements of the $u+d$ current is the same as matrix elements of the $u+d-2r$ current. Since this is a nonsinglet flavor combination formed from degenerate light quarks, it has the same renormalization factor as the isovector current. To find the disconnected light-quark contribution, we take the difference,

[TABLE]

III.6 Volume-source approach and reuse of disconnected diagrams

We evaluate our observables using quark propagators with four-dimensional volume plane-wave sources $D_{q}^{-1}(x|p)\equiv\sum_{y}D_{q}^{-1}(x,y)e^{ipy}$. For a quark-bilinear operator $\mathcal{O}=\bar{q}\Gamma q$ ($\Gamma=\gamma_{\mu}\gamma_{5}$ for the axial current), the connected contribution to the Green’s function is obtained using

[TABLE]

where $\langle\dots\rangle_{U}$ denotes the average over gauge configurations. We obtain the disconnected contribution by correlating the plane-wave-source propagators with the previously-computed disconnected loops999Recall that the loops are gauge invariant and thus do not need to be transformed to Landau gauge. $T_{\mu}^{q}(\vec{k},t)$ [Eq. (8)]:

[TABLE]

where $q$ and $q^{\prime}$ are the quark flavors of the operator and the external quark states, and we choose $p^{\prime}-p=(\vec{0},k_{4})$. Translation invariance implies that this expression is independent of $t$, and we average over all timeslices on which the disconnected loops were computed.

III.7 Results

In order to minimize cut-off effects we choose momenta on the diagonal of the Brillouin zone $p,p^{\prime}\in\frac{2\pi k}{L_{s}}(1,1,1,\pm 1)$ for $k\in\{2,3,...,8\}$. Therefore, our momenta span the range $0.6<a^{2}\mu^{2}<10$. We used for this calculation about 200 gauge configurations. This procedure involves the following steps:

Compute Landau gauge-fixed quark propagators and Green’s functions for both light and strange quarks as outlined in the previous section. Form the amputated vertex functions. 2. 2.

On the connected diagrams, impose the RI*′*-MOM or RI-SMOM vector current renormalization conditions, together with the renormalization factors from Subsec. III.2, to find estimates for $Z_{q}^{l}$ and $Z_{q}^{s}$ at each scale $\mu$. 3. 3.

Trace the axial amputated vertex functions with $P_{\mu}$ to obtain $\Sigma^{l}_{\text{conn}}$, $\Sigma^{s}_{\text{conn}}$, $\Sigma^{l}_{l,\text{disc}}$, $\Sigma^{s}_{l,\text{disc}}$, $\Sigma^{l}_{s,\text{disc}}$, and $\Sigma^{s}_{s,\text{disc}}$ at each scale $|p|$. By combining the different $\Sigma$ following Eqs. (38–42), form the matrix $Z_{A}^{-1}$. 4. 4.

Invert the matrix and evolve from scale $\mu$ to 2 GeV. 5. 5.

Optionally, convert the $Z_{A}$ matrix from the basis $\{A_{\mu}^{3},A_{\mu}^{8},A_{\mu}^{0}\}$ to $\{A_{\mu}^{u-d},A_{\mu}^{u+d},A_{\mu}^{s}\}$, using Eqs. (44) and (45). 6. 6.

Extrapolate $\mu$ to zero to remove $O(a^{2}\mu^{2})$ lattice artifacts.

For estimating the statistical and systematic errors in removing the $O(a^{2}\mu^{2})$ artifacts, we apply linear and quadratic fits for each matrix element, $Z_{A}^{ij}=c_{0}^{ij}+c_{1}^{ij}(a\mu)^{2}$ and $Z_{A}^{ij}=c_{0}^{ij}+c_{1}^{ij}(a\mu)^{2}+c_{2}^{ij}(a\mu)^{4}$. We apply these fits in different ranges of $a^{2}\mu^{2}$, all of which lie within the range $[2.5,10]$, i.e., always excluding the first two points. This fit procedure is applied to results from both RI*′-MOM and RI-SMOM schemes. We take then three best fits in each scheme (yielding six values), average all of them to get the central value and statistical uncertainty, and use the root-mean-square difference between the six values and the average to get the systematic uncertainty. Figures 10 and 11 show illustrative fits for obtaining the matrix elements in the different bases from both RI′*-MOM and RI-SMOM schemes. We obtain the following $Z_{A}$ matrices:

[TABLE]

Note that these two different matrices were obtained from independent fits to remove $O(a^{2}\mu^{2})$ artifacts, and thus they are not related exactly by Eqs. (44) and (45). For renormalizing our nucleon form factor data, we use the latter matrix. Finally, the contribution from the bare connected light axial current to the renormalized disconnected light axial current depends on the difference $Z_{A}^{u+d,u+d}-Z_{A}^{33}$, as shown in Eq. (46). In order to reduce uncertainties, we computed this difference by itself using the above procedures, and found $Z_{A}^{u+d,u+d}-Z_{A}^{33}=0.0061(18)(10)$.

From Eq. (11) and the full mass-dependent $O(a)$ improvement in Ref. Bhattacharya:2005rb , $Z_{A}^{u+d,s}$ and $Z_{A}^{s,u+d}$ first appear at two-loop order in lattice perturbation theory; since the mass-dependent part is further suppressed by $am_{s}$, it follows that these are largely sensitive to the singlet-nonsinglet difference101010This is in contrast with, e.g., $Z_{A}^{0,0}-Z_{A}^{8,8}$, which has a contribution at tree level proportional to $ab_{A}(m_{s}-m_{ud})$.. These elements are less than one percent of the diagonal ones, indicating a small difference, which is consistent with previous studies. For example, Ref. Chambers:2014pea found a singlet-nonsinglet difference $\bar{Z}_{A}-Z_{A}=0.020(3)$, using a similar lattice action. In the $SU(3)$ flavor limit, this corresponds to $\bar{Z}_{A}-Z_{A}=3Z_{A}^{s,u+d}=\frac{3}{2}Z_{A}^{u+d,s}$, so that those mixing factors are about twice as large as ours.

IV Axial form factors

IV.1 $G_{A}$ form factors

The isovector axial form factor is shown in Fig. 12 (left). From the fit, we find $g_{A}=1.208(6)(16)(1)(10)$ and $r_{A}^{2}=0.213(6)(13)(3)(0)\text{ fm}^{2}$, where the uncertainties are due to statistics, excited states, fitting, and renormalization, respectively. The dominant uncertainty is excited-state effects. The fitted value of $g_{A}$ is quite compatible with the value taken from the form factor at $Q^{2}=0$, $1.206(7)(19)(0)(10)$, with slightly smaller uncertainties. The axial charge was recently determined in a mostly independent calculation using the same ensemble Yoon:2016jzj , with somewhat higher statistics and different methodology. If we examine the bare quantity to avoid differences in renormalization factors, we get $g_{A}^{\text{bare}}=1.401(7)(18)(2)$, which differs from the result in Ref. Yoon:2016jzj , $g_{A}^{\text{bare}}=1.431(15)$, by slightly more than one standard deviation. We can compare the axial radius with the recent reanalysis of neutrino-deuteron scattering data Meyer:2016oeg that found $r_{A}^{2}=0.46(22)\text{ fm}^{2}$. Our result is slightly more than one standard deviation smaller.

Figure 12 (right) shows the light-quark isoscalar form factor $G_{A}^{u+d}(Q^{2})$. The fit yields $g_{A}^{u+d}=0.517(11)(14)(1)(3)$ and $(r_{A}^{2})^{u+d}=0.197(21)(21)(4)(0)\text{ fm}^{2}$. The statistical errors are relatively much larger than for the isovector case, and the dominant source of these errors is the connected diagrams. The uncertainty due to renormalization in $g_{A}^{u+d}$ is mostly due to the diagonal element of the renormalization matrix; the effect of mixing with strange quarks is very small.

In Fig. 13 we show the strange and light disconnected axial form factors. The strange axial form factor $G_{A}^{s}(Q^{2})$ is the most important case for mixing between light and strange axial currents, since it is small and it mixes under renormalization with $G_{A}^{u+d}(Q^{2})$, which has a contribution from connected diagrams and is much larger. The effect of this mixing is shown in the left plot: it reduces the magnitude of the form factor by up to 10%, although this effect is smaller than the total statistical uncertainty. In these plots the block-correlated nature of the statistical uncertainties is clearly visible, particularly at low $Q^{2}$: the data that are strongly correlated form clusters of nearby points, but there are large fluctuations between different clusters. This effect was previously seen in the disconnected electromagnetic form factors computed using the same dataset Green:2015wqa . Fits using the $z$ expansion to the strange and light disconnected form factors are shown in the right plot. From these fits we obtain $g_{A}^{s}=-0.0240(21)(8)(2)(7)$ and $g_{A}^{l,\text{disc}}=-0.0430(28)(46)(6)(8)$. The fit has the effect of averaging over several uncorrelated clusters of data, and produces a considerably smaller uncertainty than the value taken directly from the form factor at $Q^{2}=0$. The leading uncertainties are statistical and (for the light-quark case) excited-state effects. The uncertainty due to renormalization is dominated by uncertainty in the off-diagonal part of the renormalization matrix. We also obtain the radii $(r_{A}^{2})^{s}=0.155(73)(57)(7)(2)\text{ fm}^{2}$ and $(r_{A}^{2})^{l,\text{disc}}=0.248(57)(28)(18)(0)\text{ fm}^{2}$. Within their uncertainties, all of the squared axial radii are compatible with 0.2 fm2.

IV.2 Quark spin contributions

The axial form factors at zero momentum transfer, $g_{A}^{q}\equiv G_{A}^{q}(0)$, determined in the previous subsection, give the contribution from the spin of quarks $q$ to the proton spin. We can compare against standard experimental inputs used for phenomenological determinations of these quark spin contributions. Using isospin symmetry, the $u-d$ combination is determined from the axial charge in neutron beta decay, $g_{A}^{u-d}=1.2723(23)$ Olive:2016xmw . Our result is about 5% lower, which could be attributed to our heavier-than-physical pion mass.

The flavor nonsinglet combination $u+d-2s$ is typically obtained from semileptonic decays of octet baryons, assuming SU(3) flavor symmetry. Although there have been efforts to improve this determination using chiral perturbation theory (dating back to the original paper on the heavy baryon approach Jenkins:1990jv ), it was shown in Ref. Ledwig:2014rfa that at full next-to-leading order, there is a new low-energy constant that contributes to $g_{A}^{u+d-2s}$ but not to the octet baryon decays. Thus, in the absence of additional input, this combination cannot be predicted at NLO. The leading-order fit to octet baryon decay data Ledwig:2014rfa yields $g_{A}^{u+d-2s}=3F-D=0.608(30)$. It is therefore useful to have a lattice QCD calculation of this quantity, even for a heavy pion mass, since it will enable full NLO chiral perturbation theory analyses to be done. Our result is $g_{A}^{u+d-2s}=0.565(11)(13)$.

We find the total contribution from quark spin to the nucleon spin at $\mu=2$ GeV is $g_{A}^{u+d+s}=0.494(11)(15)$, about half. The other half must come from gluons and from quark orbital angular momentum. This is somewhat larger than results from phenomenological determinations of polarized parton distribution functions: recent analyses deFlorian:2009vb ; Nocera:2014gqa ; Sato:2016tuz give values from 0.18 to 0.28, with an uncertainty ranging from 0.04 to 0.21. There are a few possible sources for this discrepancy. First, that this is caused by our heavier-than-physical pion mass. This would require that the flavor singlet axial case be more sensitive than the isovector one to the pion mass. Second, that the unaccounted-for systematic uncertainties at this pion mass are large. These include effects due to finite lattice spacing and $O(\alpha^{2})$ corrections to the matching of the flavor singlet axial current to $\overline{\text{MS}}$. In particular, the latter does not affect the flavor nonsinglet combinations, which are in better agreement with phenomenology. A third possibility is that the phenomenological values are incorrect. The behavior at small momentum fraction $x$ is poorly constrained, and a recent estimate Kovchegov:2016weo in the large-$N_{c}$ limit of the small-$x$ asymptotics suggests that improved results at small $x$ would lead to higher values of $g_{A}^{u+d+s}$.

The individual quark contributions are summarized in Tab. 2. Our result for $g_{A}^{s}$ is compared with other lattice QCD results in Fig. 14. The results are all mutually consistent, and ours is the most precise. Our improved precision is due to much higher statistics than most previous calculations, as well as the use of a large volume and the additional constraints from data at nonzero $Q^{2}$ in the $z$-expansion fits. We also note the calculation at the physical pion mass by ETMC that was presented at Lattice 2016 Alexandrou:2016tuo , which found $g_{A}^{s}=-0.042(10)$. This differs from our result by almost two standard deviations, suggesting that the strange spin contribution to the nucleon spin becomes larger (more negative) as the light quark mass is decreased.

IV.3 $G_{P}$ form factors

Figure 15 shows the isovector induced pseudoscalar form factor $G_{P}^{u-d}(Q^{2})$. As discussed in Subsection II.6, we remove the pion pole that is present in this form factor before fitting using the $z$ expansion. With the pion pole removed, the dependence on $Q^{2}$ is much weaker. At low $Q^{2}$, there is a large systematic uncertainty from excited-state contributions. For comparison with experiment, we consider ordinary muon capture of muonic hydrogen, which (assuming isospin symmetry) is sensitive to $g_{P}^{*}\equiv\frac{m_{\mu}}{2m_{N}}G_{P}^{u-d}(Q_{*}^{2})$, where $Q_{*}^{2}=0.88m_{\mu}^{2}$. To remove the strong dependence on the pion mass arising from the pion pole, we consider Aoki:2010xg

[TABLE]

Using a modest extrapolation of our fit, we find $g_{P}^{\text{norm}}=8.47(21)(87)(2)(7)$, which is consistent with the measurement by the MuCap experiment Andreev:2012fj , $g_{P}^{*}=8.06(55)$. We can also determine the residue of the pion pole: this is related to the pion decay constant $F_{\pi}$ and the pion-nucleon coupling constant $g_{\pi NN}$ Goldberger:1958vp ,

[TABLE]

The required extrapolation in $Q^{2}$ is about twice as far as was required for $g_{P}^{*}$, but is still small compared with our probed range of $Q^{2}$. Using $F_{\pi}=106$ MeV computed on this ensemble, we obtain $g_{\pi NN}=11.5(4)(1.4)(1)(0)$. This is slightly more than one standard deviation below the recent result Baru:2010xn determined using pion-nucleon scattering lengths from measurements of pionic atoms: $g_{\pi NN}^{2}/(4\pi)=13.69(20)$, or $g_{\pi NN}=13.12(10)$. In the chiral limit, the pion-nucleon coupling constant is related to the axial charge via the Goldberger-Treiman relation, $g_{\pi NN}=g_{A}^{u-d}m_{N}/F_{\pi}$; on our ensemble the right hand side equals 12.1, and thus our precision is insufficient to resolve a nonzero Goldberger-Treiman discrepancy.

The isoscalar induced pseudoscalar form factors are shown in Fig. 16. As these contain an eta pole, we again remove the pole before fitting with the $z$ expansion. The eta mass is estimated using the leading-order relation from partially quenched chiral perturbation theory, $m_{\eta}^{2}=(m_{\pi}^{2}+2m_{\eta_{s}}^{2})/3$, yielding $m_{\eta}\approx 578$ MeV. Relative to the connected diagrams, the contributions from disconnected diagrams are not small, which is in contrast with what we saw for the $G_{A}$ form factors. This can be understood by considering the partially quenched theory, under which the connected contributions to $G_{P}^{u+d}(Q^{2})$ are equal to $G_{P}^{u+d-2r}(Q^{2})$, where $r$ is a third valence light quark, degenerate with $u$ and $d$. We would expect that this form factor has a pseudoscalar pole from the $\pi_{8}$ meson111111The presence of this pole was already argued in Ref. Liu:1991ni . (which is part of the octet of pseudo-Goldstone bosons under the exact $SU(3)$ symmetry of the valence $u$, $d$, and $r$ quarks) at $Q^{2}=-m_{\pi}^{2}$. The sharp rise of this form factor at low $Q^{2}$ is consistent with this expectation. Since the physical isoscalar form factor does not contain a pole at $Q^{2}=-m_{\pi}^{2}$, the pole must be canceled by the disconnected diagrams, which explains why the disconnected contribution to $G_{P}^{u+d}$ must also rise sharply (with opposite sign) at low $Q^{2}$. Similarly, the expectation that the octet axial current $A_{\mu}^{8}$ couples more strongly than the singlet current $A_{\mu}^{0}$ to the eta meson suggests that $G_{P}^{s}$ and $G_{P}^{u+d}$ should have opposite sign, as seen in the data.

We can attempt to quantify the couplings to the eta meson by studying the generalization of Eq. (52):

[TABLE]

where the eta decay constants are defined121212Note that using this definition for the pion decay constant would yield $f_{\pi}^{3}=\sqrt{2}F_{\pi}$, where the physical value is $f_{\pi}^{3}\approx 130$ MeV. by $\langle 0|A_{\mu}^{a}|\eta(p)\rangle=f_{\eta}^{a}p_{\mu}$ Feldmann:1998vh . As Fig. 16 shows, the extrapolation to the eta pole is rather difficult and the results have a large uncertainty. Since we have not separately computed the eta decay constants on this ensemble, we cannot determine the eta-nucleon coupling constant in this way. However, we can take the singlet-octet ratio $f_{\eta}^{0}/f_{\eta}^{8}$, which we find to be $0.96(16)(21)(4)(1)$. This is larger than expected, and three standard deviations above the value obtained from the phenomenological parameters in Ref. Feldmann:1998vh , $f_{\eta}^{0}/f_{\eta}^{8}=0.16(3)$. In particular, since our pion mass is heavier than physical, we would expect the reduced breaking of flavor $SU(3)$ symmetry to yield a value closer to zero. This unexpected behavior is likely caused by the difficulty in such a large extrapolation in $Q^{2}$; direct calculations of these decay constants such as in Ref. Michael:2013vba are much more reliable since they do not require a kinematical extrapolation. If we ignore this issue, and assume the $SU(3)$ relation $f_{\eta}^{8}=f_{\pi}^{3}$, then from $G_{P}^{8}\equiv(G_{P}^{u+d}-2G_{P}^{s})/\sqrt{6}$ we obtain an estimate for the eta-nucleon coupling constant, $g_{\eta NN}=5.2(1.0)(1.0)(0.2)(0)$.

Assuming flavor $SU(3)$ symmetry, the eta-nucleon coupling constant can also be obtained from the connected contribution to $G_{P}^{u+d}$. Provided that the considerations from the partially quenched theory are valid, the residue of the pion pole is proportional to $F_{\pi}g_{\pi_{8}NN}$, where the $\pi_{8}$-nucleon coupling constant is equal (up to $SU(3)$ breaking corrections) to $g_{\eta NN}$. Alone, the connected contribution does not benefit from the cancellation of excited-state effects with the disconnected contribution that we have seen. Therefore, to better control these effects, we determine this form factor using the summation method in the same way as $G_{P}^{u-d}$; this is shown in Fig. 17. We obtain $g_{\pi_{8}NN}=3.29(35)(45)(3)(0)$. The eta-nucleon coupling constant is not so well known phenomenologically, but both of these estimates are compatible with the value obtained using a generalized Goldberger-Treiman relation, $g_{\eta NN}=3.4(5)$ Feldmann:1999uf .

V Discussion and conclusions

As with our previous study of electromagnetic form factors Green:2015wqa , our approach of using hierarchical probing for disconnected loops and high statistics for nucleon two-point functions is effective at producing a good signal for disconnected nucleon axial form factors. In contrast with the previous study, however, we find that the gauge noise is dominant over the noise from stochastic estimation of the loops, so that further improvements in the latter would be of limited value.

A useful feature of disconnected loops is that they can be reused for calculating many different observables. We did this for computing the axial renormalization factors nonperturbatively, and we were again able to obtain a reasonable signal. At the scale $\mu=2$ GeV, the effect of mixing between light and strange axial currents is small: $G_{A}^{s}(Q^{2})$, which is most affected, is reduced in magnitude by up to 10%. The accuracy of our renormalization is limited by the unknown $O(\alpha^{2})$ term in the matching of the flavor singlet axial current to the $\overline{\text{MS}}$ scheme. Our use of two different intermediate schemes may provide some estimate of this term, but it is possible that the effect in converting between the two intermediate schemes is smaller than in converting to $\overline{\text{MS}}$. A smaller flavor-singlet renormalization factor would make both $g_{A}^{u+d+s}$ and $f_{\eta}^{0}/f_{\eta}^{8}$ more consistent with expectations. This highlights the need for higher-order conversion factors. In the flavor-nonsinglet case, these factors have been computed up to three-loop order for some operators Gracey:2003yr ; Gracey:2006zr . As lattice calculations of disconnected diagrams have made great progress, there is now a need for similar matching calculations in the flavor-singlet sector.

Since this work was performed using only one lattice ensemble, we do not provide an estimate of systematic uncertainties due to the heavier-than-physical pion mass or due to discretization effects. The former have been investigated in many lattice calculations of the isovector axial charge, where generally only modest effects have been seen. Generalizing this, we don’t expect large dependence on the pion mass for $G_{A}^{q}(Q^{2})$. On the other hand, the $G_{P}$ form factors — especially the isovector one — will have a significant dependence on light quark masses due to the presence of pseudoscalar poles. Discretization effects for this lattice ensemble have been studied in Ref. Yoon:2016jzj , where it is compared with another ensemble with similar pion mass and smaller lattice spacing. The isovector axial charge computed on the two ensembles is consistent within one standard deviation, or about 3%, which gives a rough estimate of uncertainty due to finite lattice spacing. We expect that these effects are of similar size for other nucleon matrix elements involving the axial current.

We found that the statistical correlations between the values of a form factor at different $Q^{2}$ behave differently for connected and disconnected diagrams. In the latter case, data with different spatial momentum transfers are nearly uncorrelated. This has the result of better constraining fits to the form factors; using these fits, we were able to obtain a precise value for the strange axial charge on our ensemble, $g_{A}^{s}=-0.0240(21)(11)$, which is consistent with previous lattice calculations.

For $G_{A}(Q^{2})$, the disconnected diagrams are small compared with the connected ones. For instance, $g_{A}^{u+d,\text{disc}}/g_{A}^{u+d}=-0.17$, and the strange disconnected diagrams are about half as large as the light ones. However, this is somewhat larger than we saw for the electromagnetic form factors Green:2015wqa , where the disconnected light magnetic moment, $\mu^{u+d,\text{disc}}\approx 0.11$, is about 4% of the full experimental value $\mu^{u+d}=3(\mu^{p}+\mu^{n})\approx 2.6$, and the disconnected $G_{E}(Q^{2})$ is even smaller relative to the full experimental form factor. This may change closer to the physical pion mass, since the disconnected light-quark matrix elements are expected to grow as the quark mass is decreased.

For $G_{P}(Q^{2})$, the situation is different, with disconnected diagrams not nearly as suppressed. This can be understood from the dominant influence of the pseudoscalar poles in these form factors, which leads to a significant cancellation between the connected and disconnected contributions to $G_{P}^{u+d}(Q^{2})$. As the pion mass is decreased toward the physical point, we expect that $G_{P}^{u+d}(Q^{2})$ will vary only mildly, but at low $Q^{2}$ the individual connected and disconnected contributions will become much larger since the location of the pion pole will approach $Q^{2}=0$. This growing cancellation may make it difficult to obtain a good signal for the full form factor at the physical pion mass and low $Q^{2}$.

Acknowledgements.

Computations for this work were carried out on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy, on facilities provided by XSEDE, funded by National Science Foundation grant #ACI-1053575, and at Forschungszentrum Jülich. During this research several of us were supported in part by the U.S. Department of Energy Office of Nuclear Physics under grants #DE–FG02–94ER40818 (JG, SM, JN, and AP), #DE–SC–0011090 (JN), #DE–FG02–96ER40965 (ME), #DE–FC02–12ER41890 (JL), #DE–FG02–04ER41302 (KO), #DE–AC02–05HC11231 (SS), and #DE–AC05–06OR23177 under which JSA operates the Thomas Jefferson National Accelerator Facility (KO). Support was also received from National Science Foundation grants #CCF–121834 (JL) and #PHY–1520996 (SM), the RIKEN Foreign Postdoctoral Researcher program (SS), the RHIC Physics Fellow Program of the RIKEN BNL Research Center (SM), Deutsche Forschungsgemeinschaft grant SFB–TRR 55 (SK), and the PRISMA Cluster of Excellence at the University of Mainz (JG). Calculations were performed with the Chroma software suite Edwards:2004sx , using QUDA Clark:2009wm with multi-GPU support Babich:2011:SLQ:2063384.2063478 , as well as the Qlua software suite Qlua .

Appendix A Form factor fit parameters

In this appendix we give parameters for the form factor fits and the estimated total uncertainty. Recall that we performed fits of the form

[TABLE]

where $z(Q^{2})$ is given in Eq. (12), $t_{\text{cut}}=(3m_{\pi})^{2}$, and we used the central value $am_{\pi}=0.1833$. The parameters are given in Tab. 3, where for each fit we have also given the correlation matrix. For the $G_{P}$ form factors, we give parameters for fits to $G(Q^{2})=a^{2}(Q^{2}+m^{2})G_{P}(Q^{2})$, where $m$ is either $m_{\pi}$ or $m_{\eta}$; we used the value $am_{\eta}=0.3342$. The fit curves and outer error bands for the physical form factors shown in Sec. IV (i.e., excluding the individual connected and disconnected parts) can be reproduced using the data in this table.