Theoretical predictions of the proton--neutron mass difference and measurements of the proton's charge radius require inputs from the Compton amplitude subtraction function. Model-dependent and non-relativistic calculations of this subtraction function vary significantly, and hence it contributes significant uncertainties to the aforementioned physical quantities. We report on the use of Feynman-Hellmann methods in lattice QCD to calculate the subtraction function from first principles. In particular, our initial results show anomalous high-energy behaviour that is at odds with the prediction from the operator product expansion (OPE). Therefore, we investigate the possibility that this unexpected behaviour is due to lattice artifacts, by varying the lattice spacing and volume, and comparing different discretisations of the vector current. Finally, we explore a Feynman-Hellmann implementation that is less sensitive to short-distance contributions and show that the subtraction function's anomalous behaviour can be attributed to these short-distance contributions. As such, this work represents the first steps in achieving a complete understanding of the Compton amplitude subtraction function.