# On the Feynman-Hellmann theorem in quantum field theory

and the calculation of matrix elements

###### Abstract

The Feynman-Hellmann theorem can be derived from the long Euclidean-time limit of correlation functions determined with functional derivatives of the partition function. Using this insight, we fully develop an improved method for computing matrix elements of external currents utilizing only two-point correlation functions. Our method applies to matrix elements of any external bilinear current, including nonzero momentum transfer, flavor-changing, and two or more current insertion matrix elements. The ability to identify and control all the systematic uncertainties in the analysis of the correlation functions stems from the unique time dependence of the ground-state matrix elements and the fact that all excited states and contact terms are Euclidean-time dependent. We demonstrate the utility of our method with a calculation of the nucleon axial charge using gradient-flowed domain-wall valence quarks on the MILC highly improved staggered quark ensemble with lattice spacing and pion mass of approximately 0.15 fm and 310 MeV respectively. We show full control over excited-state systematics with the new method and obtain a value of with a quark-mass-dependent renormalization coefficient.

## I Introduction

The Feynman-Hellmann theorem (FHT) in quantum mechanics relates matrix elements to variations in the spectrum Güttinger, P. (1932); Pauli (1933); Hellmann (1937); Feynman (1939):

(1) |

where the Hamiltonian is given by . This simple relation follows straightforwardly at first order in perturbation theory. The method is applicable beyond perturbation theory and is often used in lattice QCD (LQCD) calculations, for example, to compute the scalar quark matrix elements in the nucleon Procura *et al.* (2004, 2006); Alexandrou *et al.* (2008); Walker-Loud *et al.* (2009); Ohki *et al.* (2008); Young and Thomas (2010); Durr *et al.* (2012); Horsley *et al.* (2012); Freeman and Toussaint (2013); Bali *et al.* (2013); Ohki *et al.* (2013); Semke and Lutz (2012); Shanahan *et al.* (2013); Ren *et al.* (2012); Junnarkar and Walker-Loud (2013); Durr *et al.* (2016)

(2) |

for the light () and strange () quarks. Quantitative knowledge of these matrix elements is necessary for interpreting direct searches for dark matter which look for the elastic recoil of nuclei. In the scenario that dark matter is heavy and couples through the electroweak sector, the uncertainty on the strange and charm nucleon matrix elements is one of the largest uncertainties in spin-independent constraints upon direct dark matter detection Hill and Solon (2014). In particular, due to cancellations in the amplitude at the level of quarks and gluons, there is a particular sensitivity to the scalar charm quark matrix elements with current uncertainties allowing for several orders of magnitude variability in the cross section; see Fig. 3 of Ref. Hill and Solon (2014). A significant reduction over the current uncertainty in these matrix elements would be a welcome advancement for the field.

Recently, the FHT has been used to compute other nucleon matrix elements, such as the spin content of the nucleon Chambers *et al.* (2014, 2015).
More recently, a hybrid method using ideas from background field methods Fucito *et al.* (1982); Martinelli *et al.* (1982); Bernard *et al.* (1982); Detmold *et al.* (2006); Engelhardt (2007); Detmold *et al.* (2009, 2010) and the FHT has been introduced to compute few-nucleon electroweak matrix elements Savage *et al.* (2016).
An advantage of the FHT is that it relates a three-point correlation function to a change in a two-point correlation function induced by an external source.
Thus, one can take advantage of the simplified analyses of two-point functions.
Traditional lattice calculations of three-point functions, particularly those involving nucleons, face a number of challenging systematics beyond those present for two-point functions: the stochastic noise of three-point functions is more severe than the corresponding two-point functions and also three-point functions have systematic contamination from excited states which is constant in Euclidean time for fixed source-sink(insertion) separation with identical initial and final states at zero momentum transfer. Controlling these systematics requires a significant increase in the numerical cost.

Previous implementations of the FHT and related methods Chambers *et al.* (2014, 2015); Savage *et al.* (2016) are also costly, as the calculation must be performed for several values of the external parameter, .
In the case of the scalar quark matrix elements, the QCD action contains the operators of interest, . The FHT is then simply used by varying the values of the quark masses and determining the resulting variation of the spectrum, a routine step in present LQCD calculations. In the case of the nucleon spin, the operator is perturbatively added to the theory for varying values of and the resulting spectrum is computed such that can be approximated via finite difference.

In this work, we develop an improved implementation of the FHT and explore its connection with the partition function of quantum field theory. This new method offers several advantages including: an improved implementation, improved stochastic sampling over computations of equal computing time, a complete discussion of all systematics, and demonstrably rigorous control over all systematics associated with analysis of correlation functions. To demonstrate these claims, we present the formulation of our method, and perform a sample calculation of the nucleon axial-vector charge. We then discuss the generalizations and conclude.

## Ii The Feynman-Hellmann Theorem and a new method

### ii.1 The new method

Consider a two-point correlation function computed in the presence of some external source

(3) |

with the external source coupled through some bilinear current density

(4) |

and partition function in the presence of the source,

(5) |

Here, is a general field operator representing the various quantum fields of the theory. The state is the vacuum state in the presence of the external source. We denote the sourceless vacuum state, partition function, and two-point correlation function by

(6) | ||||

(7) | ||||

(8) |

respectively. The operator creates a tower of states with specified quantum numbers out of the vacuum at time , which are later destroyed by a conjugate operator at time .

We are interested in the partial derivative of this correlation function with respect to , at . This partial derivative can be built from an integral of uniform functional derivatives over the space-time volume or, if we wish for more general matrix elements (such as those involving momentum transfer), an integral over nonuniform values of . For now, we will focus on the simplest case of a constant source, .

The partial derivative of interest is related to the matrix elements of the current

(9) |

The first term is proportional to the vacuum matrix element of the current and vanishes unless the current has vacuum quantum numbers. The second term involves an integral over matrix elements involving the current and the creation/annihilation operators:

(10) |

where we have defined . The second term is related to the hadronic matrix of interest in the time region . In the other time regions, and , the current creates/destroys a tower of states that also couple to the states created by (in the case of quark bilinear operators in QCD, these are just the mesons coupled to the currents):

(11) |

Recall that the FHT relates matrix elements to derivatives of the spectrum. In Euclidean calculations, the effective mass is a derived quantity which asymptotes to the ground-state energy in the long-time limit,

(12) |

In analogy with the FHT, consider the linear response of the effective mass to the external current

(13) |

A first observation to make is that the term proportional to the vacuum matrix element in Eq. (II.1) exactly cancels in the difference in Eq. (13). The linear response of the effective mass is therefore given by

(14) |

where

(15) |

This expression can be analyzed with the usual spectral decomposition. The two-point correlation function in time-momentum space, with , is given by

(16) |

which can be obtained by inserting the identity operator

(17) |

The overlap factors are defined as

(18) |

(I) | (II) | (III) | (IV) |

The numerator in Eq. (15) can be similarly decomposed. It is useful to work in discrete Euclidean time relevant to LQCD calculations in which the numerator is

(19) |

There are four time regions to consider, depicted in Fig. 1:

(20) |

The matrix elements of interest occur in time region I, when the current is inserted between the source/sink creation/annihilation interpolating fields. The contributions from regions II, III and IV are systematic corrections which must be accounted for and controlled. In region II, the current creates/destroys a meson with the quantum numbers of the current outside the region of interest. Regions III and IV give rise to contact operators when the current insertion time occurs precisely at the creation or annihilation time of the hadron.

In the first region, , the numerator is

(21) |

where for the vacuum and otherwise. All thermal states will be exponentially suppressed and we can consider the contribution from just the vacuum state . In this case, we have

(22) |

The first term is independent of and is thus enhanced by an overall factor of . The dependence of the second term is contained entirely in the exponential factor

(23) |

where we have defined

(24) |

The contributions from region I are then

(25) |

where we have defined and written the expression to expose the symmetry. The only terms which do not also appear in the two-point correlation functions are the matrix elements of interest, . It is worth noting that the entire contribution from region I vanishes at .

Let us now consider the contribution from region II,

(26) |

To understand the systematics from this region, we cannot immediately drop the nonvacuum contributions from the thermal sum. Inserting a complete set of states, one finds

(27) |

For small , we can neglect all but the vacuum contribution in the sum over and for the first and second terms respectively. The region II contributions are then

(28) |

where we have defined

(29) |

Note the first term in Eq. (28) only contributes for scalar currents.

Finally, there is the contribution from the contact terms, regions III and IV, when or . These contact contributions are standard two-point functions with different interpolating operators

(30) | ||||

(31) |

In both terms, the thermal contributions are suppressed for all but the vacuum contribution. These two terms then contribute

(32) |

where we have defined

(33) |

Because the states are annihilated by the same operator as in the two-point function, the sum over states in Eq. (32) is over the same set of states as the two-point function but with modified overlap factors.

Putting all the regions together, the numerator is

(34) |

In order to make the expressions more manageable, we introduce the following substitutions

(35) | ||||

(36) | ||||

(37) |

leading to a two-point functional form

(38) |

and the resulting numerator expression is

(39) |

The simplicity of our method is recovered by substituting Eqs. (38) and (39) into Eq. (13). All the contributions from regions II, III and IV are encoded in the constants. These contributions, as well as the other excited-state contributions, are suppressed in the difference in the two terms in Eq. (13). It is straightforward to show that in the long-time limit, we recover the ground-state matrix element of interest

(40) |

### ii.2 Implementation

The numerical implementation of this method is straightforward. What is needed is the construction of the derivative correlation function, Eq. (II.1). We provide an explicit example of a nonscalar current with interpolating operators coupling to the proton. Standard proton creation and annihilation operators are given by

(41) |

where and are up and down quark field operators at .
, , and are spin projectors used to project onto the total spin of the proton. For example, working in the Dirac basis, the dominant spin-up local interpolating field can be constructed using Basak *et al.* (2005a, b) (where the spinor indices are labeled 1,2,3,4)

(42) |

and similar operators for the sink. Denoting the up and down quark propagators as

(43) |

the proton correlation function is given by

(44) |

The derivative correlation function () is trivially determined. Applying the partial derivative in Eq. (II.1) at the level of the path integral, one immediately observes that the derivative correlation function is obtained with the replacement of one of the quark propagators in the two-point correlation function with a Feynman-Hellmann (FH) propagator, summed over all possible insertions, where the FH propagator is simply a sequential propagator which is also summed over the current insertion time

(45) |

This extra sum leads to an stochastic enhancement of the resulting derivative correlation function as compared to the standard method with a sequential propagator. The idea of using this propagator in a two-point function can be traced back to Ref. Maiani *et al.* (1987), where the equivalent of Eq. (15) is approximated with its long-time limit and computed for various hadronic matrix elements.
The idea was further discussed in Ref. Bulava *et al.* (2012) and applied to couplings in Ref. Bernardoni *et al.* (2015).
This idea has been extended recently in Refs. de Divitiis *et al.* (2012); Alexandrou *et al.* (2014) with an application to derivatives with respect to of mesonic structure functions in Ref. de Divitiis *et al.* (2012).

For our present example, the proton derivative correlation function with a current insertion on the down quark is given by

(46) |

while for a current insertion on the up quark, one has

(47) |

These correlators are functions of the source-sink separation time in contrast to the fixed source-sink separation time dependence of the standard three-point correlators constructed with sequential propagators. It is trivial to generalize this construction to an arbitrary correlation function with the successive replacement of each quark propagator with its respective FH propagator.

For many matrix elements, one must also consider contributions from disconnected diagrams. While disconnected diagrams are stochastically noisier, we may be able to improve upon the general method as we have the freedom to compute the disconnected quark loop as a function of Euclidean time. Instead of summing over all time as in Eq. (II.1), we can explicitly choose to only sum over the time window , thus including only the contributions of interest [Eq. (25)]. It is worth exploring if this suggestion provides an improved determination of the disconnected contributions.

## Iii An application to the nucleon axial charge

We demonstrate the use of this method by computing the nucleon isovector axial charge on one of the publicly available 2+1+1 HISQ Follana *et al.* (2007) ensembles from the MILC Collaboration Bazavov *et al.* (2013), with fm, MeV and a lattice volume of .
The present work utilizes 1960 configurations with six sources per configuration.
The HISQ configurations are first gradient flowed to smear out the UV fluctuations.
Möbius domain-wall fermion (MDWF) Brower *et al.* (2012) propagators are then solved with the QUDA library Clark *et al.* (2010) with multi-GPU support Babich *et al.* (2011).
The action and efficiency of the software was described in Ref. Berkowitz *et al.* (2017a).
The valence quark mass is tuned such that the MDWF pion mass matches the taste-5 HISQ pion mass within 1%.
We then construct the isovector nucleon two-point and derivative FH correlation functions.
We further double our statistics by including the time-reversed correlation functions constructed with negative-parity nucleon operators.
The calculation requires solving both the regular and the FH propagator for each source.
The motivation and advantages of such a mixed action were described in Ref. Berkowitz *et al.* (2017a).
This action has also been used to compute the matrix element relevant to (neutrinoless double beta-decay) Nicholson *et al.* (2016).

### iii.1 Fit strategy

We extract the isovector nucleon charge by applying two analysis strategies:
we apply a standard frequentist analysis of our results as well as a Bayesian constrained fit with Gaussian priors on the two- and three-point correlation functions following the framework of Ref. Lepage *et al.* (2002). While we find consistent answers with both analysis strategies, the Bayesian analysis allows us to more rigorously and completely explore the fitting systematics, thus demonstrating the efficacy of this new method. Controlling fitting systematics for nucleons is critical for current and future LQCD calculations. Similar ideas of Bayesian constrained curve fits were explored in Ref. Chen *et al.* (2004) and recently used in Ref. Yoon *et al.* (2017).

In our present Bayesian analysis, we observe stability in all ground-state parameters under variations of the number of time slices in our data, and number of states in our fit Ansatz, demonstrating complete control over systematic uncertainty originating from the fitting procedure. We rate the quality of our fits by the -statistic, which is related to the frequentist -value as defined in Eq. (B4) of Ref. Bazavov *et al.* (2016) and consider only fits with . We select results from different fit Ansätze by comparing Bayes factors, and make distinctions between models if the Bayes factor is larger than 3. Constrained curve fits are implemented by the software package lsqfit Lepage (2016).

In the following sections, the fit procedure is first discussed for the two-point correlation function. The preferred two-point fit is then performed simultaneously with the three-point correlator leading to our final value for . We adopt the notation of placing a tilde on top of the parameter (e.g. ) to denote its prior distribution (e.g. ), and a hat (e.g. ) to denote its posterior distribution.

### iii.2 Two-point correlator analysis

#### iii.2.1 Fit strategy

The two-point fit Ansatz is given by Eq. (38) and has a sum over the infinite tower of states. We focus on the case of zero momentum insertion and drop the superscript momentum label, . In the following section, we detail how the priors are set for the infinite number of states, and quantitatively show where the series may be truncated.

We begin our two-point analysis by looking at the effective mass as defined in Eq. (12) and shown in Figs. 2 and 3. We observe a plateau indicating that the ground-state energy has a value of approximately . We set a loose prior of as indicated by the light blue band. The width of the prior is observed to be approximately 1 order of magnitude larger than the resulting posterior distribution, which is indicated by the dark blue band.

Using the estimate of , we construct the scaled two-point function in order to set prior distributions for and , which are respectively the ground-state smeared and point overlap factors. We define the scaled correlators as

(48) | ||||

(49) |

where and are the smeared-smeared and point-smeared two-point functions, respectively. From the redefinition of Eq. (38), the scaled correlators plateau to in the limit. From Figs. 4 and 5, we assign approximate values for the overlap factors, and set prior widths to half the magnitude of the expected central value leading to and . These priors are unconstraining since the widths of the posterior distributions are approximately 2 orders of magnitude smaller.

We introduce priors for the excited-state energies as energy splittings, enforcing a hierarchy of states. This is achieved by setting the energy splitting with a lognormal prior distribution. For the nucleon on a lattice with a MeV pion mass and a box size of fm, we prior the first three excited states in the following order: the Roper resonance ( from the ground state), the two-pion excitation ( from the Roper), and the one-pion excitation ( from the two-pion excitation), with prior widths which accommodate a one-pion splitting (310 MeV) within . For the fourth excited state onward, the prior for the energy splitting is set to the one-pion splitting while accommodating the two-pion splitting within .

The excited-state overlap factors can enter with an unknown sign (point-smeared); therefore, so as to not bias our expectation, we set the prior central values of the excited-state overlap factors to zero. Due to smearing, we expect less overlap with excited states, and therefore set the width of to half the value of the central value. For the point-like overlap factor, we expect equal support from all states and therefore set the width of to the value of the central value.

#### iii.2.2 Discussion

We study the stability of from the two-point correlator fits under variation of the fit region and fit Ansatz as shown in Figs. 6 and 7.

Perfect stability under varying the number of states is demonstrated from fits in Fig. 6 with , suggesting that we have controlled all the excited-state systematic uncertainty. The fits show that more excited states are needed to achieve stability when including data with smaller due to larger excited-state contamination. We also observe that for all fit variations at a fixed the Bayes factor always prefers the stable fit result with the lowest number of excited states, as indicted by the solid symbols. We conclude that the increased uncertainty and variation of the central value from fits beyond are only due to omitting parts of the usable data set. Therefore the preferred two-point fit is a seven-state fit with . We note that all fits with do not pass the -value criterion and are not considered. Similarly, Fig. 7 shows stability under varying for the preferred fit and suggests that is adequate.

The preferred fit is plotted in Figs. 2–5 as a green band and agrees with the data in the selected fit region, which is indicated by the black data points. The ground-state parameters , and are recovered in the limit from the effective mass and scaled correlators by construction. This is demonstrated by observing that the green fit band overlaps with the dark blue posterior band asymptotically at large .

### iii.3 Three-point correlator analysis

#### iii.3.1 Fit strategy

The three-point fit Ansatz is given by Eq. (39). The three-point fit Ansatz introduces additional parameters and which are not constrained by the two-point correlator. In the following section, we detail how priors are chosen for the additional parameters, and study stability under variations of fit region and number of states.

We set the ground-state prior by constructing the derivative effective mass given by Eq. (14). By construction the derivative effective mass plateaus to in the limit. Motivated by Figs. 8 and 9, we set and observe that the prior widths are approximately 1 order of magnitude larger than the width of the posterior distributions.

The prior for is chosen by observing that

(50) |

and assuming that the contribution to the infinite sum is primarily from the state due to operator smearing. This unique feature of the correlator at is exemplified in Figs. 10 and 11. Therefore we set the central value of with a prior width accommodating within .

For the excited-state priors, we set to be the same order of magnitude as . The overlap factor is set with a central value of and with the same width as . This choice follows the same logic used to set the priors for .

#### iii.3.2 Discussion

The stability of under varying the fit region and number of excited states is shown in Figs. 12 and 13. For all simultaneous two- and three-point fits, we perform the preferred two-point fit as discussed in Sec. III.2.

For fits with , Fig. 12 shows stability under varying the number of states for fits with more than five states, while the Bayes factor prefers the five-state fit (solid yellow diamonds). For fits with , we also observe stability; however the Bayes factor prefers fits with less parameters ( four states) over fits that qualitatively look more stable ( five states), suggesting that there is not enough data to support fits with a large number of states. The fit at is stable after seven states; however Eq. (50) shows that includes no information about , the parameter of interest. Therefore for this study, we choose the five-state fit with as the preferred fit. The preferred fit is shown to be insensitive under variations of as demonstrated in Fig. 13.

The posterior distributions from the simultaneous fit are used to reconstruct the fitted curve and are shown by the green bands in Figs. 8 and 9. Similar to the effective mass and scaled correlator, we recover in the limit, as shown by the exact overlap of and the green fit curve in the large- limit. Following this analysis we quote the value of the bare axial charge

(51) |

on the fm lattice with MeV. Determining the axial renormalization factor requires calculations at multiple quark masses to allow for an extrapolation to the chiral limit, which we do not perform in this work. These calculations use MDWF fermions, allowing for a determination of by computing the pion decay constant with both the conserved five-dimensional axial Ward identity and with the nonconserved four-dimensional current. We find for this ensemble

(52) |

While this is a quark-mass-dependent renormalization, we can still compare our renormalized value of with that determined using clover valence quarks on similar HISQ ensembles in Ref. Bhattacharya *et al.* (2016). They found at MeV with fm, as compared to our value of

(53) |

A more comprehensive and direct comparison with the results of Ref. Bhattacharya *et al.* (2016) will be made in a forthcoming publication Berkowitz *et al.* (2017b).
When comparing on specific ensembles, it is observed that roughly an order of magnitude less inversions are required to achieve the same statistical precision with this new method as compared to the fixed source-sink separation method.

## Iv Generalizations

The generalization of this method to nonzero momentum transfer and/or flavor-changing currents is straightforward. To inject momentum with the current, instead of summing over a constant operator, the Feynman-Hellmann propagator would be constructed with a momentum phase

(54) |

With the current implementation, each choice of momentum injection would require the computation of an additional Feynman-Hellmann propagator.

More complicated space-time-dependent current insertions can also be considered. The modification of the numerator function is also straightforward. Unlike Eq. (II.1), the energy of the incoming and outgoing states would no longer be equal, for any of the states and . The exception would be if the calculation is performed in the Breit frame when . In this case, the numerator expression from region I would still be given by Eq. (25). Otherwise, the entire contribution from region I is parametrized by the second term in Eq. (25) where the state labels and now carry information about the momentum as well.

Flavor-changing interactions are just as straightforward. In this case, the Feynman-Hellmann propagator is constructed as

(55) |

for a flavor-changing interaction , with a trivial generalization to include momentum transfer as well. As with nonzero momentum transfer, the numerator expression from region I simplifies to just the second line of Eq. (25) as the energies of the incoming and outgoing states will not match.

This method can also be extended to consider two current insertions.
A feasible technique for calculating two-current insertion matrix elements is needed to permit the LQCD evaluation of, e.g., two-photon corrections to the Lamb shift in muonic helium-3 ions Carlson *et al.* (2017), which may shed light on the proton radius problem, and the - box diagram corrections to electromagnetic structure functions Rislow and Carlson (2013) needed to improve the interpretation of results from the experiment Androic *et al.* (2013).

We first generalize Eqs. (II.1)–(5) to include multiple currents with associated couplings . In the presence of multiple currents, the modified action of Eq. (4) is

(56) |

Two-current-insertion matrix elements can then be associated with the second derivative of the effective mass,

(57) | ||||

where is a generalization of Eq. (15),

(58) |

and is the ratio of the two-current insertion matrix element to the two-point function,

(59) |

Like the first derivative of the effective mass in Eq. (14), the second derivative of the effective mass (i.e. Eq. (57)) benefits from an exact cancellation of the vacuum matrix elements of , , and . The spectral decomposition of Eq. (57) precedes in analogy with that outlined in Sec. II.1, albeit with a few added complications that we intend to address in future work.

## V Conclusions and Outlook

This work presents a computationally efficient and comprehensive implementation of lattice calculations of hadronic matrix elements utilizing ideas which can be related to the Feynman-Hellmann Theorem. In particular, the derivative correlation function was analytically derived such that a background field is not explicitly needed, reducing the overall number of propagator solves required to determine the linear response of the theory to an external current.

The example calculation of the nucleon axial charge demonstrates that from the calculation of a low-statistics run, we are able to achieve a 2.1% uncertainty on , including all systematic errors from the fitting procedure. The derivation of systematic effects associated with this method, including excited-state contamination, mesonic propagating modes, and signal from contact operators were completely presented up to terms suppressed by