Solving a polynomial equation involving gamma function

I am trying to solve an equation having degenerate limit of Fermi-Dirac integral.

My code basically is

ClearAll["Global*"]
FD1[d_, η_] := η^(d + 1)/
Gamma[d + 2] ;  (* Defining the Fermi-Dirac integrals in degenerate *)

Solve[FD1[(d-2)/2, ηs] + FD1[(d-2)/2, ηs - vd] ==
2 FD1[(d-2)/2, η0], {ηs}]

I want to solve for $$\eta_S$$ in terms of $$\eta_0$$ and $$v_d$$. Is it possible to solve it using this method? (And maybe expand the solution in terms of $$vd$$ as $$vd$$ is small and we can only pick the first few terms)

I would recommend using the new in 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 zeroth order approximation of ηs when vd 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}}