  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
$N_f=2+1+1$ 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 $g_A = 1.213(26)$ with a quark-mass-dependent renormalization coefficient.
