Consider the following example:
{a, b}::Indices. test:=\sqrt{-\det{G_{a b}}; substitute(_,$G_{a b}->g_{a b}+\epsilon h_{a b}$);
How to expand a series to a first-order term for $\epsilon$ using the function map_sympy?
map_sympy
I might be mistaken, but I don't think that could be done. To my knowlegde sympy does not manage tensor components. However, I leave you the link to the manual (series expansions).
sympy
However, we did a sort of manual approximation to get a result, check the user notebooks, or equivalently the section 5 of our pre-print Cadabra and Python algorithms in General Relativity and Cosmology I: Generalities.
Although you won't find an answer to your question there, I hope they serve to give you alternative perspective.
Cheers.
Oh, I've seen the information you provided before, but as you said, the problem has not been solved. Still appreciate the information.