I have a function such as
$$ g(\epsilon) = \frac{1}{2}\left[ 1 - \frac{\epsilon - \mu}{\sqrt{(\epsilon - \mu)^2 + \Delta^2}}\right]. $$
When $\Delta \rightarrow 0$, $g(\epsilon)$ reproduces very well the Fermi-Dirac distribution graphically
$$ f(\epsilon) = \frac{1}{e^{(\epsilon - \mu)/kT} + 1}$$
for $T \rightarrow 0$ (which is approximately the step function at this regime).
I have tried to expand $g(\epsilon)$ in order to arrive in $f(\epsilon)$ without success. Is there a possibility find $f(\epsilon)$ from $g(\epsilon)$ analytically?
You might well be envisioning the wrong asymptotic. W.l.o.g., take $\Delta=kT$. Both functions g and f are really functions of just one variable, $$ x\equiv \frac{\Delta}{\epsilon-\mu },\\ g= G(x) = {1\over 2} \left ( 1-{1\over \sqrt{1+x^2} } \right ), \\ f=F(x)= {1\over e^{1/x} +1} ~~ . $$
You wish to investigate the function $$ H(x)=G(x) - F(x), $$ which, as @J.G. reminded you, goes to 0 for the "majority" cases $|x|\ll 1$, (i.e. both $x>0$, $G(x)\to 0; ~~F(x)\to 0$, and also $x<0$, $G(x)\to 1; ~~F(x)\to 1$ ); and also the "transitional" ones $|x|\gg 1$, $G(x)\to 1/2; ~~F(x)\to 1/2$.
The only place where it does not vanish, e.g. $H(1)\approx -0.12$, is in the region $x=O(1)$, where no violent behavior is observed—recall the transitional region $|x|\gg 1$ is $O(x^{-1})$ safe!
All you are reassured about is the cases where you actually use the approximation, $|x|\ll 1$, which cannot fail to work, $H(\pm 0.01)\approx 0.003\%$.