Solve a polynomial equation that involves the gamma function

I would recommend using the new M12 function AsymptoticSolve for this. Your equation:

eqn = FD1((d-2)/2, ηs) + FD1((d-2)/2, ηs - vd) == 2 FD1((d-2)/2, η0);

We need to find the zero order approximation of ηs when vd it is small:

Simplify(Solve(eqn /. vd -> 0), (η0 | vd) ∈ Reals)

{{ηs -> η0}}

Now use AsymptoticSolve:

AsymptoticSolve(eqn, {ηs, η0}, {vd, 0, 5})

{{ηs ->
vd / 2 + ((-6 + d) (-2 + d) (-1 + d) vd ^ 4) / (
1536 η0 ^ 3) + ((2 – d) vd ^ 2) / (16 η0) + η0}}

Appendix

If you have an earlier version of Mathematica, so you do not have access to AsymptoticSolve, you can try using the cloud instead. For example, define:

asymptoticSolve(args__) := CloudEvaluate(System`AsymptoticSolve(args))

Then use asymptoticSolve instead of AsymptoticSolve.