I have a bivariate polynomial for each $n=0,1,2...$ $$ f_n(x,y)=\sum _{k=0}^n \frac{(-1)^k}{2 k+1} \binom{n}{k} \left(x ^2-y ^2\right)^{2 n-2 k}\left([y ( x^2 -1) +x(1 -y^2 )]^{2 k+1}\\ \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad-[y ( x^2 -1) -x(1 -y^2 )]^{2 k+1}\right) $$ looking at low orders explicitly in Mathematica reveals that this polynomial contains $n+1$ factors of $(1-y^2)$. I know that one factor comes from the difference of two odd powers, but cannot see where the other $n$ factors come from. The factorized expressions appear quite complex, for example: $$ f_3=\frac{2}{35} \left(y ^2-1\right)^4 x \left(-35 y^6+35 \left(5 y ^2-1\right) x ^{10}+7 \left(15 y ^4-73 y ^2+3\right) x ^8+(y -1) (y +1) \left(5 y^4-388 y^2+5\right) x ^6-7 y ^2 \left(3 y ^4-73 y ^2+15\right) x ^4+35 y ^4 \left(y^2-5\right) x ^2+35 x ^{12}\right) $$

Assuming $f_n$ contains $(1-y^2)^{n+1}$, How can I find the formula for $f_n$ in its factorized form, for general $n$?

If it is any use, the series coefficients are related to the following integral: $$ \sum _{k=0}^n \frac{(-1)^k}{2 k+1}\binom{n}{k} z^{2k+1}=\int_0^z(1-u^2)^n\mathrm{d}u $$