Analytical calculation of alchemical derivatives based on auxiliary density perturbation theory is described, coded, and validated. For the case where the nucleus is a hydrogen atom and the nuclear charge is changed from 1 to 0, it turns out that a good estimate of the proton binding energies can be obtained very efficiently. First-order results correspond exactly to the molecular electrostatic potential evaluated at the hydrogen nucleus location (removing self-repulsion), in agreement with previously reported extensive studies. Therefore, the second-order results reported here are refinements in accuracy that finally allow a quantitative exploration of differential acidity. Furthermore, the recently reported h function is produced in its analytical form as a byproduct and local descriptor associated with the proton binding energy values found with this approach. In an example application, proton binding energies are computed for a family of imidazolium derivatives to demonstrate the capabilities and the stability of the method with respect to changes in basis set or exchange-correlation functional.