Nucleon Electric Dipole Moment from the θ Term with Lattice Chiral Fermions

We calculate the nucleon electric dipole moment (EDM) from the θ term with overlap fermions on three domain wall lattices with diﬀerent sea pion masses at lattice spacing 0.11 fm. Due to the chiral symmetry conserved by the overlap fermions, we have well deﬁned topological charge and chiral limit for the EDM. Thus, the chiral extrapolation can be carried out reliably at nonzero lattice spacings. We use three to four diﬀerent partially quenched valence pion masses for each sea pion mass and ﬁnd that the EDM dependence on the valence and sea pion masses behaves oppositely, which can be described by partially quenched chiral perturbation theory. With the help of the cluster decomposition error reduction (CDER) technique, we determine the neutron and proton EDM at the physical pion mass to be d n = − 0 . 00148 (14) (31) ¯ θ e · fm and d p = 0 . 0038 (11) (8) ¯ θ e · fm. This work is a clear demonstration of the advantages of using chiral fermions in the nucleon EDM calculation and paves the road to future precise studies of the strong CP violation eﬀects.

Introduction: Symmetries and their breaking are essential topics in modern physics, among which the discrete symmetries C (charge conjugation), P (parity), and T (time reversal) are of special importance.This is partially because the violation of the combined C and P symmetries is one of the three Sakharov conditions [1] that are necessary to give rise to the baryon asymmetry of the universe (BAU).However, despite the great success of the standard model (SM), the weak baryogenesis mechanism from the CP violation (¨CP ) within the SM contributes negligibly (∼ 16 orders of magnitude smaller than the observed BAU [2][3][4][5][6]).This poses a hint that, besides the possible θ term in QCD, there could exist beyond-standard-model (BSM) sources of ¨CP and thus the study of ¨CP plays an important role in the efforts of searching for BSM physics.
The electric dipole moment of nucleons (NEDM) serves as an important observable to study ¨CP .The first experimental upper limit on the neutron EDM (nEDM) was given in 1957 [7] as ∼ 10 −20 e•cm.During the past 60 years of experiments, this upper limit has been improved by 6 orders of magnitude.The most recent experimental result of the nEDM is 0.0(1.1)(0.2) × 10 −26 e•cm [8], which is still around 5 orders of magnitude larger than the contribution that can be offered by the weak ¨CP * jianliang@scnu.edu.cnphase.Currently, several experiments are aiming at improving the limit down to 10 −28 e•cm in the next ∼10 years.This still leaves plenty of room for the study of ¨CP from BSM interactions and the QCD θ term.
As a reliable nonperturbative method for solving the strong interaction, lattice QCD provides us the possibility of studying the nucleon EDM (NEDM) from first principles and with both the statistical and systematic uncertainties under control.To be specific, lattice QCD can be used to calculate the ratio between the neutron and proton EDM induced by strong ¨CP and the parameter θ, which is the most crucial theoretical input to determine θ from experiments.
Many lattice calculations have been carried out on this topic.However, there was a watershed in 2017 when it was pointed out [9] that all the previous lattice calculations, e.g.[10][11][12][13][14], used a wrongly defined ¨CP form factor such that all of those old results need a correction.Although the fixing is numerically straight forward, none of the previous lattice calculations gives statistically significant results after the fixing, leaving a great challenge to the lattice community.Since then, several attempts [15][16][17][18] have been made to tackle the problem, but the signal-to-noise ratios of the new results are still not satisfying, and no calculation performed directly at the physical point gives nonzero results.
A possibility to bypass this difficulty is to perform the computations with several heavier pion masses and ex-trapolate to the physical point.However, only with chiral fermions can a correct chiral limit be reached at finite lattice spacings.Otherwise, extrapolating to the continuum limit for each pion mass becomes an inevitable prior step before a reliable chiral extrapolation, which complicates the calculation and potentially leads to hard-to-control systematic uncertainties.The best result, so far, of this approach, using clover fermions, obtained a 2-sigma signal [16].
In this article, we demonstrate that using chiral fermions to extrapolate to the physical point from heavier pion masses is the most efficient choice to study NEDM on the lattice at the current stage.We employ 3 gauge ensembles with different sea pion masses ranging from ∼300 to ∼600 MeV and we use 3 to 4 valence pion masses on each lattice.Therefore, we can study both the valence and sea pion mass dependence of the NEDM and better control the chiral extrapolation.The results we obtain at the physical pion mass are d n = −0.00148(14) (31) θ e•fm and d p = 0.0038 (11) ( 8) θ e•fm for neutron and proton, respectively.
Nucleon EDM and the θ term: The QCD Lagrangian in Euclidean space with the θ term reads (detailed conventions can be found in the Supplemental Materials [19]): where where θ is the original coefficient of the θ term and M is the quark mass matrix generated by the spontaneous breaking of SU (2) × U (1) in the electroweak sector.For simplicity, we will not distinguish θ and θ in the following content.A crucial point is that, if Det [M ] = 0, phase of the U A (1) transformation is arbitrary, which means one can always find a chiral rotation that lets θ = 0, leaving no net effect of ¨CP .This indicates a zero NEDM in the chiral limit [20], which poses a very strong constraint in the chiral extrapolation numerically.However, as mentioned before, for lattice fermions which violate the chiral symmetry this constraint cannot be used at finite lattice spacing.
Given that θ is small, one can expand the theta term in the action in the path integral and obtain the correlation functions and matrix elements to the leading order in θ as θ ... θ = ... + iθ ...Q t , where |0 θ denotes the vacuum with the θ term (namely, the θ vacuum), and is the topological charge of the gauge field geometrically.Based on this expansion, the ¨CP electromagnetic (EM) form factor F 3 (q 2 ) can be extracted from normal and Q t weighted nucleon matrix elements with initial momentum Tr Γ e M (2) , α where the matrix elements are with V µ being the EM current operator, Γ e = 1+γ4 2 is the unpolarized spin projector, Γ i = −iγ 5 γ i Γ e the polarized projector along the i'th direction, q 2 = (p f − p i ) 2 = −Q 2 the momentum transfer, and q i the nonzero component of the momentum transfer.The above formalism is the same for both neutron and proton.In the end, the nucleon EDM can be extracted from the ¨CP form factor F 3 (q 2 ) in the forward limit for neutron and proton respectively using An interesting fact, as seen in Eq. ( 2), is that the neutron ¨CP form factor at the zero momentum transfer limit, F 3,n (0) has no ¨CP angle α 1 dependence since G E,n (0) = 0, and thus one actually needs no information about M (2)Q in the neutron case.
Numerical setups: This study is carried out on three 2 + 1-flavor RBC/UKQCD gauge ensembles of domain wall fermions [21] with the same lattice spacing 0.1105(3) fm and lattice volume 24 3 × 64 but different sea quark masses.Using the overlap fermion action [22] on the HYP (hyper-cubic) smeared [23] gauge links, multiple partially quenched valence quark masses (as listed in Table I with other parameters) are calculated utilizing the multi-mass inversion algorithm; thus both the sea and valence pion mass dependencies of NEDM can be studied and the chiral extrapolation can be more reliable.Generally, using overlap fermions can be O(100) times more costly compared to the traditional Wilson-like discretized fermion actions.To improve the computational efficiency, 12-12-12 grid sources with Z 3 -noise and Gaussian smearing are placed at t src = 0 and t src = 32 in one inversion with randomly chosen spatial positions on different configurations, and low-mode substitution (LMS) [24] is applied to suppress the statistical contamination between different source positions.We also use the stochastic sandwich method (SSM) [25] with LMS to make the cost of using multiple nucleon sinks be additive instead of multiplicative.We use 8 sets of source noises and 16 sets of sink noises (for each of the source-sink separations 6a, 7a, and 8a) to improve the statistics.Five nonzero momentum transfers are calculated such that we can reliably do the q 2 extrapolation to get F 3 (0); the details of the q 2 extrapolation are given in the Supplemental Materials [19].
CDER improvement and results: To further suppress the statistical uncertainty of M (2)Q and M (3)Q , we take advantage of a technique called cluster decomposition error reduction (CDER) for the disconnected insertion [26].As illustrated in Fig. 1, we write the total topological charge as the summation of the local charge density q t (x) derived from the overlap operator [27,28] as q t (x) = 1  2 Tr [γ 5 D ov (x, x)], where the trace is over the color-spin indices, and convert the two-point function weighted with the total topological charge Q t into a summation of the three-point functions involving q t (x) where χ is the nucleon interpolating operator, G denotes the source grid, and x = (t f , x).We then use the cluster decomposition property to limit the sum to a range commensurate with the correlation length which reduces the variance by a volume factor [26].In Eq. ( 6), R is the 4-dimensional truncated size of the topological operator, δm is the effective mass gap between the with different mπ,v and mπ,s = 339 MeV.We can see that the value saturates at R ∼ 9a.
nucleon and its excited states, and m η is the mass of the pseudoscalar meson η.
Similarly, the three-point function with Q t can be converted into a four-point function with q t (x) where y = (t c , y), and δE( q) is the energy gap of the nucleon and its excited states with 3-momentum q at the sink.Using Eqs. ( 6) and ( 7), the ¨CP form factor F 3 can be calculated as a function of cutoff R. Due to the cluster decomposition principle, operators far enough separated have exponentially small correlation.When the distance between operators is larger than the correlation length ∼ 1/m η , the signal falls below the noise while the errors still accumulate in the disconnected insertions [26].So we bind the topological charge to the sink of the nucleon in the three-point functions or to the inserted currents in the four-point function to see if a proper cutoff R exists, such that the physics is not altered while the errors can be reduced.
Then we do the two-state fit to eliminate the excitedstate contamination of nucleon matrix elements at each value of R, and obtain F 3 (Q 2 ) as a function of R. The corresponding systematic uncertainty is estimated to be the difference between the value from the two-state fits and that from single-exponential fits using only the middle point at different separations.Taking F 3,n (Q 2 = 0.2 GeV 2 ) at m π,s = 339 MeV and different m π,v as an example (shown in Fig. 2), the central value starts to saturate at around R = 9a ∼ 2/m η as expected.Since the R dependence for different pion masses are similar, we choose R c = 9a as our optimal cutoff in the neutron case.For the proton, we use R c = 10a.The systematic uncertainty of this cutoff will be estimated by two independent ways: 1) taking the difference between the value at the cutoff R c and the constant fit result with R ≥ R c ; 2) fitting the correlation between the topological charge density and the current operator in the nucleon state to an exponential form first, and then taking the summation of the correlation in the tail R ≥ R c .Either way suggests a ∼12% systematic uncertainty.
Benefited from CDER, the data points of MeV, while there is no significant deviation from a linear shape.Thus we use a linear fit for the extrapolation to Q 2 = 0, and estimate the corresponding systematic uncertainty to be the difference between the extrapolated value and the data value with the smallest Q 2 .
After the Q 2 → 0 extrapolation, the final chiral extrapolation of the neutron EDM is shown in the upper panel of Fig. 4 with both valence and sea pion mass dependencies.We observe that the partially quenched data behave differently from those with unitary points in the lower panel.The former tend to move away from zero as the valence quark mass decreases.Using the overlap fermion allows us to fit our data with the partially quenched chiral perturbation form [29] at finite lattice spacing, where c 1,2,3,n/p are free parameters.Our lattice data are well fitted with χ 2 /d.o.f.= 1.2, and our numerical results suggest that the different valence and sea quark mass dependence is consistent with the chiral perturbation expression.It is also interesting to point out that the chiral log term is crucial to ensure that the NEDM approaches zero in the chiral limit of both the valence and sea quark masses.With the zero NEDM constraint  The chiral extrapolation of dn/θ on both the sea and valence quark masses (upper panel) and on only the unitary points (lower panel).
at the chiral limit, our interpolated result for neutron is d n = −0.00148(14), where the statistical uncertainty is less than 10%.This is quite an improvement from the 2 σ statistical error in Ref. [16].
We also carry out another chiral extrapolation using only the unitary pion mass points, as shown in the lower panel of Fig. 4. It gives d n = −0.00142(20)θ, which is consistent with the prediction using partially quenched data points but with larger statistical uncertainty.We take the difference between the extrapolated results with and without partially quenched data points as an estimation of the systematic uncertainty in the chiral extrapolation.
The proton EDM and its systematic uncertainties can be obtained with a similar procedure.More detailed discussion on the fits, systematic uncertainty estimation, and proton EDM can be found in the Supplemental Materials [19].
Summary: We calculate the nucleon electric dipole moment with overlap fermions on 3 domain wall lattices at lattice spacing 0.11 fm.Since the overlap fermion preserves chiral symmetry, we have well-defined topological charge and the chiral extrapolation is carried out reliably without the need of doing continuum extrapolations first.We have in total 3 sea pion masses and 10 partially quenched valence pion masses in the chiral fitting and find that the EDM dependence on the sea and valence pion masses behaves oppositely.
With the help of the cluster decomposition error reduction (CDER) technique, we determine the neutron and proton EDM at the physical pion mass point to be d n = −0.00148(14) (31) θ e•fm and d p = 0.0038 (11) (8) θ e•fm, respectively.The two uncertainties are the statistical uncertainty and the total systematic uncertainty from the excited-state contamination, the CDER cutoff, and the Q 2 and chiral extrapolations.By using the most recent experimental upper limit of d n , our results indicate that θ < 10 −10 .This work demonstrates the advantage of using chiral fermions in the NEDM calculation and paves the road to future precise studies of the strong ¨CP effects.

I. CONVENTIONS AND FORMALISM
In this part of Supplemental Materials, we list our notations and conventions in a very detailed manner, which we think is quite worthwhile since the final sign of EDM depends directly on the conventions used.

I.1. Gamma Matrices
First, for the gamma matrices in Minkowski space, we use where η µν = (+, −, −, −) is the corresponding metric tensor.Similarly, we have, for the Euclidean ones, with η E µν = (+, +, +, +).Our choice is to let γ E 4 = γ 0 while γ E i = −iγ i .For the momentum we have Then, with the above definitions, we come to the following convention of the spinors and we define in our notations.

I.2. QCD Lagrangian with the θ Term
The Minkowski QCD Lagrangian reads where the covariant derivative is D µ ≡ ∂ µ − igA µ with a minus sign in front of A µ .Along with this convention, we use To have the QCD Lagrangian in Euclidean space, we first notice And for the gauge fields, the conversion is the same as that of p E and p: Combining the above relations, we come to and Plugging in the conversions of the gamma matrices, we have The Minkowski field tensor satisfies where , and It is easy to check that and such that Then, we finally reach the form of the QCD Lagrangian in Euclidean space When the θ term is taken into consideration, in Minkowski space, we have L → L + L θ and where F µν = µνρσ F ρσ and q t is the topological charge density.Based on the above conversions, we have and in the end

I.3. Spinors Under the θ Vacuum
Now we have determined the Lagrangian in Euclidean space.In the following part of the Supplemental Materials, we will work in the Euclidean space and omit the superscript E unless otherwise specified.
After the θ term is plugged in, the P and CP symmetries are broken.The normal Dirac equation and spinor definition should be modified.The new Dirac equation reads the new spinors can be expressed as and such that we have Also, we define the overlapping factor where χ is the nucleon interpolating filed operator and |N is the corresponding nucleon state.Then, under the θ vacuum we define Similarly, we have

I.4. Form Factors
In Minkowski space, we use the following electromagnetic form factor decomposition where F 1 and F 2 are the Pauli and Dirac form factors respectively, q = p − p with p the momentum of the outgoing nucleon (ū (p )) and p the momentum of the incoming nucleon.For the Minkowski case, with our conventions we have and using the Dirac equation (p / − m) u = 0 we get On the other hand, with the Euclidean notation, we have −σ E µν q E ν = i γ E µ q / E − q E µ .And similarly So in order to have consistent results for both Minkowski and Euclidean space, one should use −σ E µν q E ν under our convention: For the ¨CP case, we have an additional form factor F 3 N.B., when taking the phase carried by the ¨CP spinors into consideration, this CP odd form factor should be modified as well.The relation between the correct ¨CP form factor under the θ vacuum F 3 and F 3 can be retrieved by considering the parity transformation of the normal spinors and the ¨CP ones Specifically, we have

I.5. Correlation Functions
In general, path integrals under the θ vacuum can be estimated by employing the Taylor expansion in θ and keeping only the leading term where F E,µν is the total topological charge and q t is the charge density.Correlation functions can therefore be accessed by For example, the two-point functions can be expressed as where G θ 2 , G 2 , and G Q 2 are two-point functions evaluated with the θ term, normal two-point functions, and two-point functions weighted by the topological charge, respectively.Since, where Z and Z are the sink and source overlapping factors and m and E are the nucleon mass and energy, and we can get Here we are assuming t is large enough so that only the ground state survives to simplify the equations.These two-point correlation functions offer to a way of determining the ¨CP angle α 1 : where is the unpolarized projector and I 4 is the 4 by 4 identity matrix.The angle α 1 is actually the leading coefficient of the spinor dependence on θ, which, in some sense, measures the ¨CP effect of the θ term.
For the three-point function case, similarly, we have The normal three-point function is where the subscripts i and f are for the initial and final nucleons respectively.Denoting the common factor ZZ † e −E f (t f −tc) e −Eitc m 2 E f Ei = A for simplicity, we have The relation between the correlators and the form factors will be derived as follows.In general, the nucleon matrix elements in the three-point correlation functions can be decomposed into CP even and CP odd form factors W even µ and W odd µ as and Thus we have and So by doing a similar subtraction, we arrive at This is what the three-point correlator weighted by the topological charge looks like, and is what we use to extract the ¨CP form factors.

II. COMPARISON OF DIFFERENT TOPOLOGICAL CHARGE DEFINITIONS
In this study, we use overlap fermions as valence quarks.The overlap Dirac operator D ov satisfies the Ginsparg-Wilson relation, which ensures the lattice version of chiral symmetry at finite lattice spacing a.Moreover, since the modified quark field ψ = (1 − 1/2D ov )ψ is used for the chirally regulated current operators and interpolating fields, the effective quark propagator is then 1/(D c + m q ), where m q is the current quark mass and D c = D ov /(1 − 1/2D ov ) anticommutes with γ 5 , i.e., {D c , γ 5 } = 0 [31].This is the same form as in the continuum and the eigenvalues of D c are purely imaginary.Actually, it has been shown that all the current algebra is satisfied with overlap fermions at finite a.In particular, the anomalous Ward identity (AWI) has been proven by Peter Hasenfratz [23] for D ov with chiral axial vector current.And we have also shown numerically [32] that the normalization factor Z A obtained from the axial Ward identity in the isovector case is the same (within error) as the one from the AWI in the singlet case.
Geometrically, the θ term is related to the topological charge of the gauge field Usually, the F F definition of the topological charge with unsmeared gauge fields suffers from large UV effects and cannot give integer total topological charge values on the lattice (a review on this topic can be found in [33]).One way to solve the problem is to use the gradient flow to smooth the gauge fields and to get renormalized topological charges [34][35][36].Since we are using a lattice chiral fermion, we have an alternative way to obtain the topological charge.According to the Atiyah-Singer index theorem, the topological charge equals the numerical difference between the left-handed zero-modes of D ov and the right-handed zero-modes, that is, Q t = n − −n + , which ensures integer topological charge on each configuration with no additional renormalization.This definition is theoretically the same as the definition from the overlap Dirac operator where the trace over all color, spin and space-time indices of D ov can be estimated through noise sources.And this D ov definition can also be used to define the topological charge density q t (x).The topological charge term is essential in the nEDM calculation and the overlap definition reduces the subtleties in the evaluation of the topological charges, which is another benefit of using chiral fermions.In the left panel, the distribution with label Nν 0 corresponds to the topological charges from counting the zero modes, which should be the same as the one with label Dov.The nuanced difference between them comes from the fact that Dov is estimated by noise and has statistical fluctuations.The distribution with label F F corresponds to that using the gluonic definition with t f = 4a 2 .The brown color is the overlay of orange and blue.In the right panel, the topological susceptibility from the F F definition is plotted as a function of the flow time, while the topological susceptibility from the overlap definition is shown as a band.
It is interesting to note the difference between topological charges from the overlap definition and those from the gluonic definition with long enough gradient flow until integer topological charge values are reached.We find that, as shown in the left panel of Fig. 5, the total topological charge on individual configurations with the gluonic definition is not necessarily the same as the one with the overlap definition.This is actually natural as they involve different regulations.However, the topological charge distributions over different gauge configurations in a given ensemble are similar.All the distributions are approximately symmetric with central value at around zero, and it seems the gluonic definition gives more zero charges.Now, a further question is whether they will lead to consistent physical results at finite lattice spacing.
For the purpose of checking physical results, we calculate the topological susceptibility on the same lattice The right panel of Fig. 5 shows that at large flow time t f , the value of the topological susceptibility from the gluonic definition tends to approach that from the overlap definition.However, it is found that, even at t f = 6a 2 , the χ t value from the gluonic definition is still around 10% higher than that from the overlap definition although there is a gentle trend that the central values will be closer as the flow time t f is larger still.For the study at only one lattice spacing, it is hard to justify a precise choice of t f that is large enough.On top of this, there is O(a 2 ) error.Accordingly, in order to avoid such unnecessary systematic uncertainties, we use the overlap definition of the topological charge in our calculation.Another conclusion that can be drawn here is that the specific topological charge value on each single configuration has not much effect on the physical correlations; only the distribution matters.

III.1. Extracting Form Factors
To calculate the ¨CP form factor, we need to the make three-point function to two-point function ratios and where Γ i is the polarized projector and J µ stands for the current insertion.If we write down the explicit form of the correlators, we have, e.g., in the CP even case, Here again we assume t is large enough to simplify the equations.Details of dealing with the excited-states contamination are discussed in the systematic uncertainty section.The additional overlapping and kinematic factors in Eqs.(63, 64) are cancelled with proper combination of two-point correlation functions.Please note that in our numerical setup we always set the initial momentum p i = 0 in three-point functions.With proper selection of the momentum p f , polarization Γ i and current insertion J µ , the ratio gives the desired nucleon matrix element for particular form factors (or combinations of form factors).The relation between the corresponding form factors and the setup of the ratios are derived as follows.
For the normal EM case, we choose unpolarized projection and vector current γ 4 , which gives (in our momentum setup) where G E ≡ F 1 − q 2 4m 2 F 2 is the electric form factor.The last step used the fact that the momentum transfer and Therefore we have and We can also choose polarized projection (Γ i ≡ −i 1+γ4 2 γ 5 γ i ) and the γ j current: where G M ≡ F 1 + F 2 is the magnetic form factor, or unpolarized projection and γ i : These ratios can be used to extract the CP conserved form factors.For the ¨CP case, we can choose the polarized projection and γ 4 , which turns out to be An important fact about this ratio is that the neutron form factor F 3,n (0) has no α 1 dependence since G E,n (0) = 0.This means that one needs no information about α 1 or the other CP-even form factors if one focuses only on the neutron case.Similarly, we can also use which prefers giving the combination of α 1 F 2 + F 3 rather than (Γ i , γ 4 ) in our calculation.

III.2. Summary on the systematic uncertainties
In this study, the main sources of systematic uncertainties are the two-state fits of the three-point (four-point) function to two-point function ratios, the momentum extrapolation, the use of the CDER technique, the final chiral extrapolation, and the finite lattice spacing effect.
1) Two-state fit: The systematic uncertainty from the two-state fit is estimated by the difference between the twostate fitted values and the results from single-exponential fits using only the middle point at different separations.Usually, one compares the two-state fits results and the values of the middle data point at the largest separations to estimate the systematic uncertainty.In our case, since we are using relatively small source-sink separations, we fit the middle points to a simplified form C 0 + C 1 e −mt f to account for the excited-state effect on different separations t f .Then, we consider the distribution of the difference between the two-state results and C 0 's (as shown in the right panel of Fig. 6), and take the 1 σ width (68% probability) to be the the final systematic uncertainty, which is determined to be 13%.
2) Momentum extrapolation: Considering the systematic uncertainty from the momentum extrapolation, although we have 5 momentum transfers, the data points show no significant deviation from a linear shape due to the large uncertainties, so we use a linear fit for the extrapolation and estimate the corresponding systematic uncertainty to be the difference between the extrapolated value and the data value with the smallest momentum transfer.An example plot can be found in Fig. 7. Similar to the two-state fit case, the systematic uncertainty is estimated to be 10% by taking the 1 σ width of the error distribution shown in the right panel of Fig. 7.A fit with an additional Q 4 term results in no significant difference.
3) CDER technique: The systematic uncertainty due to the use of the CDER technique is a crucial one.The key of d n normalized by the number of equivalent R's which is in fact the correlation in terms of the 4-D distance r between the topological charge operator and the current operator, since where N |q(x + r)J µ (x)|N denotes the nucleon matrix element that encodes the correlation.This panel demonstrates that the correlation decays exponentially and there is indeed a finite correlation length.The optimal cutoff is chosen to be R 0 = 9a.The systematic error can obtained by two ways.One is to take the difference between the value at R 0 = 9a and the constant fitted value after that cutoff.From data such as that in the left panel the systematic uncertainty is estimated to be ∼10-15% in this way.The other way is to fit the correlation to an exponential form first, and then put the fitted correlation in the summation d n (R) ∼ |r|>R0 q(0 + r)J µ (0) to estimate the contribution from the truncated tail.In this way, with the correlation data such as that in the right panel, the corresponding systematic uncertainty is estimated to be ∼10%.So the two methods give consistent systematic uncertainties and we choose ∼12% to be our final estimation.4) Chiral extrapolation: For the systematic uncertainty from the chiral extrapolation, we take the difference of the extrapolations with and without partially quenched data points to be our estimation.As shown in Fig. 9 and Fig. 10 (the chiral fits for proton), the difference is around 3%.The small systematic uncertainty of chiral interpolation is understandable since the chiral limit provides a very strong constraint to the interpolation.
The total systematic uncertainty is found to be 21%, which is simply calculated by quadrature from all the systematic uncertainties.

Figure 1 .
Figure 1.Illustration of the CDER technique used when computing the correlation functions with the local topological charge summed inside the sphere with radius R.

Figure 3 .
Figure 3.The Q 2 dependence of F3,n with mπ,v ∼ mπ,s = 339 MeV.The green band shows a linear fit in Q 2 while the red band shows the fit with an additional Q 4 term.

Figure 4 .
Figure 4.The chiral extrapolation of dn/θ on both the sea and valence quark masses (upper panel) and on only the unitary points (lower panel).

Figure 5 .
Figure 5. Topological charge distributions over gauge configurations with different definitions (left panel) and the topological susceptibility (right panel).In the left panel, the distribution with label Nν 0 corresponds to the topological charges from counting the zero modes, which should be the same as the one with label Dov.The nuanced difference between them comes from the fact that Dov is estimated by noise and has statistical fluctuations.The distribution with label F F corresponds to that using the gluonic definition with t f = 4a 2 .The brown color is the overlay of orange and blue.In the right panel, the topological susceptibility from the F F definition is plotted as a function of the flow time, while the topological susceptibility from the overlap definition is shown as a band.

Figure 7 .
Figure 7.An example of momentum transfer extrapolation of F3,n with R = 9a and mπ,v ∼ mπ,s = 339 MeV (left panel) and the systematic uncertainty distribution over CDER cutoffs and pion masses (right panel).In the left panel, blue points are lattice data and The green band shows a linear fit in Q 2 while the red band shows the fit with an additional Q 4 term.

Table I .
Parameters of the RBC/UKQCD ensembles: label, sea and valence pion masses, and the number of configurations.label mπ,s (MeV) mπ,v (MeV) N cfg Figure 6.An example of a two-state fit of F3,n with Q 2 = 0.2 GeV 2 , R = 9a and mπ,v ∼ mπ,s = 339 MeV (left panel) and the systematic uncertainty distribution over different momentum transfers, CDER cutoffs and pion masses (right panel).