I can't replicate the procedure on page 9 and 10 of this paper (Hill, 1878):
[...] If the second radical in this expresion is expanded in a series of proceeding according to descending powers of $R$ and the fist term omitted $\frac{1}{R}$, [...] the following expression is obtained:
Without the constants, the expression is:
$$\Omega = \frac{1}{\sqrt{(x+R)^2+y^2+z^2}}$$
And the expansion Hill made is:
$$\Omega=\frac{1}{R^3}\left [ x^2 - \frac{1}{2}(y^2+z^2) \right ] + \frac{1}{R^4}\left [ x^3 - \frac{3}{2}x(y^2+z^2) \right ]+\frac{1}{R^5}\left [ x^4 +3x^2(y^2+z^2) + \frac{3}{8}(y^2+z^2)^2 \right ]+\frac{1}{R^6}\left [ x^5 +5x^3(y^2+z^2) + \frac{15}{8}x(y^2+z^2)^2 \right ] +\dots$$
I tried Fourier series with all variable changes that occurred to me . Also try generalized Newton Binomial, but did not get anything.
I put the expression in Mathetmatica, whit the Series function, with $1/R$ as the variable, and I get the expression of Hill, but I have to do by myself. I trying Laurent Series, but I realy do not how to do for this expression.
Can you help me, please?
Starting from$$\Omega = \frac{1}{\sqrt{(x+R)^2+y^2+z^2}}= \frac{1}{\sqrt{R^2+2xR+x^2+y^2+z^2}}$$ let $a=2x$ and $b=x^2+y^2+z^2$ which make $$\Omega = \frac{1}{\sqrt{R^2+aR+b}}=\frac 1 R\frac{1}{\sqrt{1+\frac aR+\frac b{R^2}}}$$ Set $t=\frac aR+\frac b{R^2}$ and use Taylor series $$\frac{1}{\sqrt{1+t}}=1-\frac{t}{2}+\frac{3 t^2}{8}-\frac{5 t^3}{16}+\frac{35 t^4}{128}-\frac{63 t^5}{256}++\frac{231 t^6}{1024}+O\left(t^7\right)$$ Replace $t$ by $\frac aR+\frac b{R^2}$ and use the binomial theorem truncating at $O\left(\frac{1}{R^7}\right)$ $$t=\frac aR+\frac b{R^2}$$ $$t^2=\frac{a^2}{R^2}+\frac{2 a b}{R^3}+\frac{b^2}{R^4}$$ $$t^3=\frac{a^3}{R^3}+\frac{3 a^2 b}{R^4}+\frac{3 a b^2}{R^5}+\frac{b^3}{R^6}$$ $$t^4=\frac{a^4}{R^4}+\frac{4 a^3 b}{R^5}+\frac{6 a^2 b^2}{R^6}+O\left(\frac{1}{R^7}\right)$$ $$t^5=\frac{a^5}{R^5}+\frac{5 a^4 b}{R^6}+O\left(\frac{1}{R^7}\right)$$ $$t^6=\frac{a^6}{R^6}+O\left(\frac{1}{R^7}\right)$$ All of this makes $$\Omega=\frac{1}{R}-\frac{a}{2 R^2}+\frac{3 a^2-4 b}{8 R^3}+\frac{12 a b-5 a^3}{16 R^4}+\frac{35 a^4-120 a^2 b+48 b^2}{128 R^5}-\frac{63 a^5-280 a^3 b+240 a b^2}{256 R^6}+O\left(\frac{1}{R^7}\right)$$ Now, replacing $a$ by $2x$ and $b$ by $x^2+y^2+z^2$ leads to $$\Omega=\frac{1}{R}-\frac{x}{R^2}+\frac{2 x^2-y^2-z^2}{2 R^3}+\frac{-2 x^3+3 x y^2+3 x z^2}{2 R^4}+\frac{8 x^4-24 x^2 y^2-24 x^2 z^2+3 y^4+6 y^2 z^2+3 z^4}{8 R^5}+\frac{-8 x^5+40 x^3 y^2+40 x^3 z^2-15 x y^4-30 x y^2 z^2-15 x z^4}{8 R^6}+O\left(\frac{1}{R^7}\right)$$ which is not exactly the same as what you wrote : look at the coefficient of $R^{-4}$; the sign is changed. I also noticed a few other differences in the higher order terms.
Edit
WE could make things slightly easier defining $c=y^2+z^2$ which makes $b=\left(\frac a 2\right)^2+c$. Doing it, the expansion becomes $$\Omega=\frac{1}{R}-\frac{a}{2 R^2}+\frac{a^2-2 c}{4 R^3}-\frac{a \left(a^2-6 c\right)}{8 R^4}+\frac{a^4-12 a^2 c+6 c^2}{16 R^5}-\frac{a \left(a^4-20 a^2 c+30 c^2\right)}{32 R^6}+O\left(\frac{1}{R^7}\right)$$ Replace $a$ by $2x$ and $c$ by $(y^2 +z^2)$.