How can we generalize factorisation of $(a+b)^n-(a^n+b^n)$ 
How can we generalize factorisation of
$$(a+b)^n-(a^n+b^n)\,?$$
where $n$ is an odd positive integer.

I found the following cases:
$$(a+b)^3-a^3-b^3=3ab(a+b)$$
$$(a+b)^5-a^5-b^5=5 a b (a + b) (a^2 + a b + b^2)$$
$$(a+b)^7-a^7-b^7=7ab(a+b)(a^2+ab+b^2)^2$$
Exapansions and factoring seem terrible. Maybe there is a good way for generalization.
The question is, can we generalize these factorisations?
$$(a+b)^n-(a^n
+b^n)$$
for all odd integer $n$?
 A: For $n=9$ we have
$$
(a+b)^9-(a^9+b^9)=3(3a^6 + 9a^5b + 19a^4b^2 + 23a^3b^3 + 19a^2b^4 + 9ab^5 + 3b^6)(a + b)ab,
$$
where the first part is irreducible.
A: There doesn't seem to be any general factored form, but on the few polynomial factors you mention in the question, we can answer exactly when they appear. From symmetry and homogeneity, it is enough to consider the simpler $$p(x)=(1+x)^n-(1+x^n)$$ as with $x = b/a$, we have $a^n p(x)$ is the expression you want to factor.

*

*Now $p(0)=0 \implies x$ is always a factor of $p(x) \implies ab$ is always a factor of the original bivariate expression.  From the Binomial theorem, it is also clear that no higher power of $x$ or $(ab)$ will be a factor.

*Noting $p(-1)= -1-(-1)^n \implies x+1$ is a factor when $n$ is odd, or equivalently $(a+b)$ is a factor for odd $n$.  From $p'(x) = n(1+x)^{n-1}-nx^{n-1}, p'(-1) \neq 0$ so the multiplicity of this factor is always $1$.

*Let $\omega$ be a non-real cube root of unity, i.e. $\omega^2+\omega+1=0$.  Then
$p(\omega) = (-1)^n(\omega)^{2n}-1-\omega^n$.
Case $n=3k$: we have $p(\omega)=(-1)^n$
Case $n=3k+1$: we have $p(\omega)= (-1)^n\omega^2+\omega^2$, which is zero iff $n$ is also odd, i.e. for $n=6k+1$.
Case $n=3k+2$: we have $p(\omega)= (-1)^n\omega+\omega$, which vanishes iff $n$ is also odd, i.e. for $n=6k+5$.

If $p(\omega)=0$ we have $1+x+x^2$ as a factor, or equivalently $a^2+ab+b^2$ is a factor iff $n \equiv \pm1 \pmod 6$. Let us also check multiplicity for these cases using $p'(x) = n(1+x)^{n-1}-nx^{n-1} \implies p'(\omega)=n\left( (-\omega^2)^{n-1} - \omega^{n-1}\right)$:
Case $n=6k+1$: we always have $p'(\omega) = 0$, so always $(a^2+ab+b^2)^2$ is always a factor.  Also $p''(\omega) = n(n-1) \neq 0$, so the multiplicity is exactly $2$.
Case $n=6k+5$: we have $p'(\omega) = n\left(\omega^2-\omega \right)\neq 0$ so here the multiplicity is always $1$.
Summarising all of that, for all odd $n$, we have $ab(a+b)$ as a factor, with no higher power for any of the terms.  In case $n=1 \pmod 6$, we have additionally $(a^2+ab+b^2)^2$ as a factor (with no higher power possible) and in case $n=5 \pmod 6$ we have just $a^2+ab+b^2$ as a factor.  In no other case will $a^2+ab+b^2$ be a factor.
The above few factors with limited multiplicity are clearly not enough to provide required degree for higher order polynomials, so there would necessarily be other polynomial factors.
A: For any case obviously you have $a$ and $b$ in each term since you've removed both the $a$-only and the $b$-only terms from $(a+b)^n$. Also $a=-b$ is always a root, so you'll always get an $(a+b)$ term. So $ab(a+b)$ is a given for any odd $n$.
The roots of $a^2+ab+b^2$ are $a={1\over2}(-b\pm{\sqrt{b^2-4b^2}})=-b\omega,-b\bar\omega$ where $\omega$ is a complex number such that $\omega^3=1$, so if $n=5$ then if $a=-\omega b$ then $$a^5+b^5=b^5(1-\omega^5)=b^5(1-\omega^2)$$ and $$(a+b)^5=b^5(1-\omega)^5=b^5{\bar\omega}^5=b^5\bar\omega^2$$ since $$1-\omega=\bar\omega,\tag{1}$$
and
$$\omega^2=-\bar\omega.\tag{2}$$
You'll get a similar result if $n$ is coprime with $3$, but you don't get that root for $9$. So for example if $n=11$ you get
$$
11ab(a+b)(a^2+ab+b^2)(a^6+3a^5b+7a^4b^2+9a^3b^3+7a^2b^4+3ab^5+b^6)$$
and if $n=13$ you get
$$
13ab(a+b)(a^2+ab+b^2)^2(a^6+3a^5b+8a^4b^2+11a^3b^3+8a^2b^4+3ab^5+b^6)
$$
(obtained using Wolfram Alpha).
In general, using (1),
$$(b-\omega b)^n=b^n(1-\omega)^n=b^n(\bar\omega)^r
$$
where $n=3q+r$, $0<r<3$, and $q,r\in\mathbb{N}$, and
$$
b^n+(-\omega b)^n=b^n(1-\omega^r)=b^n(\bar\omega)^r
$$
if $r=1,2$, so
$$
(b-\omega b)^n-(b^n+(-\omega b)^n)=0
$$
and $a=-\omega b$ is a solution if $r=1,2$ and hence $a^2+ab+b^2$ is a factor.
Since the complex conjugate of a complex root of a real equation is always a root too, then obviously $a=-\bar\omega b$ is also a root.
A: 
A generalisation for non-negative integer $n$ is provided in this answer. The expansion there can be written as
\begin{align*}
&\,\,\color{blue}{(a+b)^n-\left(a^n+b^n\right)}\\
&\quad\,\,\color{blue}{=ab(a+b)\sum_{k=1}^{\lfloor n/2\rfloor}\left(\binom{n-k}{k}+\binom{n-k-1}{k-1}\right)(-ab)^{k-1}(a+b)^{n-2k-1}}\tag{1}\\
\end{align*}
where we use the binomial identity
$
2\binom{n-k}{k}-\binom{n-k-1}{k}=\binom{n-k}{k}+\binom{n-k-1}{k-1}
$.

We check formula (1) by taking $n=7$. We obtain
\begin{align*}
&\color{blue}{(a+b)^7-\left(a^7+b^7\right)}\\
&\qquad\quad=ab(a+b)\sum_{k=1}^{3}\left(\binom{7-k}{k}+\binom{6-k}{k-1}\right)(-ab)^{k-1}(a+b)^{6-2k}\\
&\qquad\quad=ab(a+b)\left(\left(\binom{6}{1}+\binom{5}{0}\right)(a+b)^4
-\left(\binom{5}{2}+\binom{4}{1}\right)ab(a+b)^2\right.\\
&\qquad\qquad\qquad\qquad\quad\left.+\left(\binom{4}{3}+\binom{3}{2}\right)(ab)^2\right)\\
&\qquad\quad=ab(a+b)\left(7(a+b)^4-14(ab)(a+b)^2+7(ab)^2\right)\\
&\qquad\quad=7ab(a+b)\left((a+b)^4-2(ab)(a+b)^2+(ab)^2\right)\\
&\qquad\quad=7ab(a+b)\left((a+b)^2-ab\right)^2\\
&\qquad\quad\,\,\color{blue}{=7ab(a+b)\left(a^2+ab+b^2\right)^2}
\end{align*}
in accordance with OPs derivation.
