Suppose $f,g : R^N \to R$ and we have an integral of the form
$$ I = \int d^N \mathbf{x} \, g(\mathbf{x}) \delta(f(\mathbf{x})), $$ with the Dirac delta function restricting the domain of integration to the submanifold on which $f$ vanishes (assume this manifold is nice and smooth, etc.). Suppose additionally that we can write $$f(\mathbf{x}) = f_0(\mathbf{x}) + \epsilon f_1(\mathbf{x}),$$ where $\epsilon$ is small in some sense and $f_0$ is not a constant.
Under what conditions, if any, can we expand the delta function as a formal Taylor series $$\delta(f(\mathbf{x})) \overset{?}{=} f_0(\mathbf{x}) + \epsilon f_1(\mathbf{x}) \delta'(f_0(\mathbf{x})) + \frac{1}{2} \epsilon^2 f_1(\mathbf{x})^2 \delta''(f_0(\mathbf{x})) +\dots,$$ where $\delta'$ is the distributional derivative, and then evaluate $I$ term by term using integration by parts? If this equation doesn't make sense, can it be modified somehow so that it does?