I have more results to share. I can accept the phase modulation artifacts, they change the NCO period by less than 7%, so it is not a big issue. I can even lie through teeth and proclaim it to be an EMI-reducing spread spectrum modulator or another continuum transfunctioner. Nonetheless, for completeness I have created a generic second order biquad structure, filled it with various coefficients and simulated in double precision to avoid the most obvious numeric instabilities. I have checked the Gaussian, Butterworth, Chebyshev and elliptic approximations and the -3dB frequencies in the range of 1..50Hz.
As it can be expected, the ripple level of these filters is superb, especially in the case of the elliptic filter for 5.9Hz, which happens to have a zero at 100Hz, BUT the loop is much harder to stabilize. The lock occurs only for a relatively narrow range of gains, while the reference exponential smoother works correctly within the range as wide as an order of magnitude. Lower gain just means increased constant phase error, while too much gain means increased distortion level. The filters are also more agile, but they always introduce a tiny phase lag, while the smoother with sufficient gain doesn't. For some reason 10Hz passband gives the best results, lower bandwidth filters converge too slowly or not at all, while the wider ones are robust, but their ripple level is comparable to that of the smoother. So since the increased complexity doesn't buy me much, I'll stick to the embarassingly primitive smoother.
Moreover, it is possible to add a frequency-dependent gain compensation based on linear interpolation, which, too, works unexpectedly well. The PLL locks in the range of 45..55Hz with a negligible phase offset.
There is something deeply mysterious about this low alpha exponential smoother, its robustness is on the verge of offensiveness, given its crude nature. I think I am onto something, but don't know what that something might be.
Best regards, Piotr