Constraining beyond the Standard Model nucleon isovector charges

Historically nuclear and neutron beta decays have played an important role in determining the vector-axial (V-A) structure of weak interactions and in shaping the Standard Model (SM). However, more recently, neutron and nuclear $\beta$-decays can also be used to probe the existence of beyond the Standard Model (BSM) tensor and scalar interactions. The interaction of the $W$ boson with the neutron and proton during neutron $\beta$-decay is proportional to the matrix element of flavour changing vector and axial-vector currents between the initial neutron state and final proton state, with coupling constants $g_{A}/g_{V}=1.2756(13)$ [1]. It has been identified that the potential existence of BSM tensor and scalar couplings would provide additional contributions to neutron $\beta$-decay [2]. These new BSM contributions are proportional to analogous matrix elements of flavour-changing tensor or scalar operators. To gain sensitivity to these effects the majority of previous and proposed neutron beta decay studies aim to determine one or more of the correlation coefficients included in the differential decay rate for a beam of polarised neutrons [2]:

where $\bm{\sigma}_{n}$ is the neutron spin, $p_{e}$ is the momentum of the electron and $p_{\nu}$ is the momentum of the neutrino with energies $E_{e}$ and $E_{\nu}$, respectively, and $E_{0}$ is the end-point energy of the electron. In the SM, $\xi=G_{F}^{2}V_{ud}^{2}(1+3\lambda^{2})$, where $\lambda=g_{A}/g_{V}$ is the ratio of the axial-vector and vector coupling constants and $G_{F}$ is the Fermi constant. The neutron decay observables include, $a$, the neutrino-electron correlation coefficient, $b$, the Fierz interference term, $A$, the beta asymmetry, and $B$, the neutrino asymmetry. Within the SM, the correlation coefficients $a,A$ and $B$ depend solely on the ratio of the axial-vector and vector coupling constants, $\lambda=g_{A}/g_{V}$. However the parameter, $b$, is included to account for the case of the hypothetical scalar or tensor couplings in addition to the (V-A) interaction of the SM. Many experiments are underway worldwide with the aim to improve the precision of measurements of these neutron decay observables, two importantly being the neutrino asymmetry $B$ [3], and the Fierz interference term $b$ [4, 5]. The parameter $b$ has linear sensitivity to BSM physics through [6]:

where $\epsilon_{T}$ and $\epsilon_{S}$ are the new-physics effective couplings and $g_{T}$ and $g_{S}$ are the tensor and scalar nucleon isovector charges. Here $b^{\text{BSM}}_{v}$ is a correction term to the neutrino asymmetry correlation coefficient, $B$, and $b^{\text{BSM}}$ is an addition to the Fierz interference term $b$ in Eq. 1. Data taken at the Large Hadron Collider (LHC) is currently looking at probing scalar and tensor interactions at the $\lesssim 10^{-3}$ level [7]. However to fully assess the discovery potential of experiments at the $10^{-3}$ level it is crucial to identify existing constraints on new scalar and tensor operators.

Another quantity of interest is the neutron electric dipole moment (EDM), which is a measure for CP violation. In extensions of the Standard Model quarks acquire an EDM through the interaction of the photon with the tensor current [8]. The contribution of the quark EDMs, $d_{q}$, to the EDM of the neutron, $d_{n}$, is related to the quark tensor charges, $g_{T}^{q}$, by [9, 10, 11]:

Here $d_{u},~{}d_{d},~{}d_{s},$ are the new effective couplings which contain new CP violating interactions at the TeV scale. The current experimental data gives an upper limit on the neutron EDM of $|d_{n}|<1.8\times 10^{-26}e$.cm [12]. In calculating the tensor charges and knowing a bound on $d_{n}$, we are able to constrain the couplings, $d_{q}$, and hence BSM theories.

In recent years there has been an increase in interest from lattice QCD collaborations in calculating the axial, scalar and tensor isovector charges due to their importance in interpreting the results of many experiments and phenomena mediated by weak interactions [13, 14, 15, 16, 17, 18, 19]. The QCDSF/UKQCD/CSSM collaborations have an ongoing program investigating various hadronic properties using the Feynman-Hellmann theorem [20, 21, 22, 23, 24, 25, 26, 27]. Here we extend this work to a dedicated study of the nucleon tensor, scalar and axial charges. We discuss a flavour symmetry breaking method to systematically approach the physical quark mass. We then extend this existing flavour breaking expansion to also account for lattice spacing and finite volume effects to quantify systemic uncertainties. Finally, we look at the potential impact of our results on measurements of the Fierz interference term and the neutron EDM.

For this work we use gauge field configurations that have been generated with $N_{f}=2+1$ flavours of dynamical fermions, using the tree-level Symanzik improved gluon action and non-perturbatively $\mathcal{O}(a)$ improved Wilson fermions [28]. In our simulations, we have kept the bare average quark mass, $\bar{m}=(m_{u}+m_{d}+m_{s})/3$, held fixed approximately at its physical value, while systematically varying the quark masses around the $SU(3)$ flavour symmetric point, $m_{u}=m_{d}=m_{s}$, to extrapolate results to the physical point [29]. We also have degenerate $u$ and $d$ quark masses, $m_{u}=m_{d}\equiv m_{l}$. The coverage of lattice spacings and pion masses is represented graphically in Fig. 1.

Further information about these ensembles is presented in Table 1 and Appendix A, Table 7. We have five lattice spacings, $a=0.082,~{}0.074,~{}0.068,~{}0.059,~{}0.052$ fm [30], enabling an extrapolation to the continuum limit as well as three lattice volumes, $32^{3}\times 64$, $48^{3}\times 96$ and $64^{3}\times 96$, allowing an extension to the flavour-breaking expansion, which describes the quark mass-dependence of the matrix elements, to also account for lattice spacing and finite volume effects. We also use a bootstrapping resampling technique to compute all statistical uncertainties in our study.

In order to compare with existing results in the literature we use the renormalisation constants given in Table 2. Table 2 summarises the renormalisation constants at each value of $\beta$ after chiral and continuum extrapolation across multiple masses with conversion from R$\text{I}^{\prime}$-MOM to $\overline{\text{MS}}$ at $\mu~{}=~{}2$ GeV. The renormalisation constants are calculated following the method in Ref. [32] and the results first appeared in Ref. [31].

The Feynman-Hellmann (FH) theorem is used to calculate hadronic matrix elements in lattice QCD through modifications to the QCD Lagrangian. The expression for the FH theorem in the context of field theory is [20]:

where $S$ is a modified action of our theory so that it depends on some parameter $\lambda$, $S\rightarrow S(\lambda)$ and $E_{H,\lambda}(\vec{p})$ is the energy of a hadron state, $H$. This result relates the derivative of the total energy to the expectation value of the derivative of the action with respect to the same parameter.

The dependency of the fits on the time ranges used is a source of systematic uncertainty. To address these issues, we use a weighted averaging method on the fit results to limit the impact of the fit window selection. The weighted averaging method we use is a simplified variation of that outlined in detail in Ref. [53] and has similarities to that proposed in Ref. [54]. We proceed by determining the energy shifts, $\Delta E=E_{\lambda}-E$, by fitting the ratio of perturbed to unperturbed correlation functions using Eq. 17 for a variety of different Euclidean time fit windows. The largest time slice employed in each fit for each ensemble and operator is fixed to be the last time slice before the signal is lost due to statistical noise. For example, in Fig. 2 this would be chosen to be $t_{\text{max}}\approx 17$. The start of the fit range, $t_{\text{min}}$, is varied between $t_{\text{min}/a}=6,~{}7,~{}8,~{}9,~{}10$ for ensembles with $\beta=5.40,~{}5.50,~{}5.65,~{}5.80,~{}5.95,$ respectively and up to the largest value of $t_{\text{min}}$ such that no less than four time slices are used in a fit. By adjusting the minimum time of the fit range, $t_{\text{min}}$, based on the lattice spacing of each ensemble, we are ensuring that each fit starts at an earlier scale. In the following we refer to the value of $\Delta E$ for a single fit, $f$, as $E^{f}$. Each fit result is then assigned a weight:

where $f$ labels the choice of fit range specified by $t_{\text{min}}$ for a fixed $t_{\text{max}}$, $p_{f}=\Gamma(N_{\text{dof}}/2,~{}\chi^{2}/2)/\Gamma(N_{\text{dof}}/2)$ is the p-value of the fit and $\delta E^{f}$ is the uncertainty in the energy shift, $E^{f}$, for fit $f$. Taking a weighted average of the $N$ fit findings, $E^{f}$, provides the final estimate of the energy shift, $\overline{E}$, and associated uncertainty $\delta\overline{E}$,

The total error $\delta\overline{E}$ describes the combined statistical uncertainty on $\overline{E}$ plus the systematic uncertainty arising from the choice of fit range. The separating of this error into $\delta_{\text{stat}}\overline{E}$ and $\delta_{\text{sys}}\overline{E}$ only partially separates statistical and systematic uncertainties because $\delta_{\text{stat}}\overline{E}$ includes statistical errors plus systematic uncertainties related to fluctuations among the $\delta E^{f}$. The final estimate, $\overline{E}$, provides an estimate of the energy of the hadron with reduced systematic bias arising from choice of fit window. Fig. 3 shows the proton effective energy shift for the ratio (Eq. 17), using the standard definition of an effective mass. The final estimate of the energy shift, $\overline{E}$, when using the weighted averaging method is indicated by the red band. Fig. 3 also includes a bar graph for the weights assigned to each fit value.

The QCD interaction is flavour-blind, which means that the only distinction between quark flavours comes from the quark masses when we disregard the electromagnetic and weak interactions. The theory behind these interactions is easiest to understand when all three quark flavours share the same mass, as this allows us to use the full power of $SU(3)$ flavour symmetry. Here we have kept the bare quark mass, $\bar{m}=(m_{u}+m_{d}+m_{s})/3$, held fixed at its physical value, while systematically varying the individual quark masses around the $SU(3)$ flavour symmetric point, $m_{u}=m_{d}=m_{s}$, in order to constrain the extrapolation to the physical point. In this work we simulate with degenerate $u$ and $d$ quark masses $m_{u}=m_{d}\equiv m_{l}$, restricting ourselves to $n_{f}=2+1$.

When $SU(3)$ is unbroken all octet baryon matrix elements of a given octet operator can be expressed in terms of just two couplings $f$ and $d$. However, once $SU(3)$ is broken and we move away from the symmetric point we can construct quantities ($D_{i}$, $F_{i}$) which are equal at the symmetric point but differ in the case where the quark masses are different. The theory behind constructing these quantities is described in detail in Ref. [55] and is summarised below. The result of constructing these quantities leads to ‘fan’ plots, with slope parameters ($r_{i}$, $s_{i}$) relating them. Following the method in Ref. [55] we use the $SU(3)$ expansion to extrapolate the nucleon charges to the physical point.

In this work, we describe the quark mass dependence of the hadronic matrix elements by a perturbation in the quark masses about an $SU(3)$ symmetric point. This perturbation generates a polynomial expansion in the quark mass differences (i.e. the $SU(3)$ breaking parameter) and therefore appears distinct from a chiral formulation that generates nonanalytic behaviour (e.g. logarithms) in the vicinity of the 2- or 3-flavour chiral limits. However, it has been demonstrated in Ref. [56], that by expanding the logarithmic features about a fixed quark mass point (such as the chosen $SU(3)$ symmetric point), the infrared singularities reveal themselves in the high-order terms of the polynomial expansion — hence demonstrating that the group-theoretic expansion does encode the same physics that appears in the logarithms. A detailed numerical investigation exploring the numerical convergence from both limits goes beyond the present work. Here we assess the convergence of our expansion empirically, subject to the precision of our numerical results.

The flavour breaking expansion described in Section VI only accounts for the quark mass-dependence of the matrix elements. However, in order to quantify systematic uncertainties, here we extend this method to also account for lattice spacing, finite volume effects and second order mass terms. As we are performing a global fit over all ensembles, we are now able to place constraints on the second-order mass terms, which means that all fits will now incorporate a term of order $\mathcal{O}(\delta m_{l}^{2})$. These fits also include corrections with respect to $a$, $a^{2}$ and $m_{\pi}L$. In order to perform a global fit across all masses we substitute the quantity $\delta m_{l}$ from here on with:

where the pseudoscalar mass flavour singlet, $X_{\pi}^{2}$, is given by:

By determining $\delta m_{l}$ to now be dimensionless and given in terms of physical quantities we are now able to combine results from different lattice spacings. The fit used for the singlet quantities $X_{D}$ and $X_{F}$ are extended to [56]: