The paper Ballistic Coefficient Estimation for Reentry Prediction of Rocket Bodies in Eccentric Orbits Based on TLE Data gives the formula (3) to calculate the time derivative of the semi-major axis due to atmospheric drag only:
$$\frac{da}{dt}|_{drag} = 2 \frac{a^2}{\sqrt{\mu p}}\left[ f_{r_{drag}} e \sin \theta + f_{t_{drag}} \frac{p}{r} \right], \tag{3}$$ where $p$ is the semilatus rectum, $\theta$ the true anomaly, and $f_{r_{drag}}$ and $f_{t_{drag}}$ the acceleration due to drag in radial and transverse direction, respectively.
Please, could someone explain how the formula is obtained?
EDIT: could someone confirm that "transverse direction" means orbit normal?
EDIT #2: I mistakenly thought that "transverse direction" meant "orbit normal", but it's the in-track component instead (normal to the radius vector in the orbital plane, in the direction of motion). Now the formula works like a charm.