Seeking for help to find a formula for $\int_{0}^{\pi} \frac{d x}{(a-\cos x)^{n}}$, where $a>1.$ When tackling the question, I found that for any $a>1$,
$$
I_1(a)=\int_{0}^{\pi} \frac{d x}{a-\cos x}=\frac{\pi}{\sqrt{a^{2}-1}}.
$$
Then I started to think whether there is a formula for the integral
$$
I_n(a)=\int_{0}^{\pi} \frac{d x}{(a-\cos x)^{n}}, $$
where $n\in N.$
After trying  some substitution and integration by parts, I still failed and got no idea for reducing the power n. After two days, the Leibniz Rule for high derivatives come to my mind.
Differentiating $I_1(a)$ w.r.t. $a$ by $(n-1)$ times yields
$$
\displaystyle \begin{array}{l}
\displaystyle \int_{0}^{\pi} \frac{(-1)^{n-1}(n-1) !}{(a-\cos x)^{n}} d x=\frac{d^{n-1}}{d a^{n-1}}\left(\frac{\pi}{\sqrt{a^{2}-1}}\right) \\ \displaystyle 
\int_{0}^{\pi} \frac{d x}{(a-\cos x)^{n}}=\frac{(-1)^{n-1} \pi}{(n-1) !} \frac{d^{n-1}}{d a^{n-1}}\left(\frac{1}{\sqrt{a^{2}-1}}\right) \tag{*}\label{star}
\end{array}
$$
I am glad to see that the integration problem turn to be merely a differentiation problem.
Now I am going to find the $(n-1)^{th} $ derivative by Leibniz Rule.
First of all, differentiating $I_1(a)$ w.r.t. $a$ yields $$
\left(a^{2}-1\right) \frac{d y}{d a}+a y=0 \tag{1}\label{diffeq}
$$
Differentiating \eqref{diffeq} w.r.t. $a$ by $(n-1)$ times gets $$
\begin{array}{l}
\displaystyle \left(a^{2}-1\right) \frac{d^{n} y}{d a^{n}}+\left(\begin{array}{c}
n-1 \\
1
\end{array}\right)(2 a) \frac{d^{n-1} y}{d a^{n-1}}+2\left(\begin{array}{c}
n-1 \\
2
\end{array}\right) \frac{d^{n-2} y}{d a^{n-2}}+x \frac{d^{n-1} y}{d a^{n-1}}+(n-1) \frac{d^{n-2} y}{d a^{n-2}}=0
\end{array}
$$
Simplifying, $$
\left(a^{2}-1\right) y^{(n)}+(2 n-1) ay^{(n-1)}+(n-1)^{2} y^{(n-2)}=0 \tag{2}\label{diffrec}
$$
Initially, we have $ \displaystyle y^{(0)}=\frac{1}{\sqrt{a^{2}-1}}$ and $ \displaystyle y^{(1)}=-\frac{a}{\left(a^{2}-1\right)^{\frac{3}{2}}}.$
By \eqref{diffrec}, we get $$
y^{(2)}=\frac{2 a^{2}+1}{\left(a^{2}-1\right)^{\frac{5}{2}}}
$$ and $$
\displaystyle  y^{(3)}=-\frac{3 a\left(2 a^{2}+3\right)}{\left(a^{2}-1\right)^{\frac{7}{2}}}
$$
Plugging into \eqref{star} yields $$
\begin{aligned}
\int_{0}^{\pi} \frac{d x}{(a-\cos x)^{3}} &=\frac{\pi}{2} y^{(2)}=\frac{\pi\left(2 a^{2}+1\right)}{2\left(a^{2}-1\right)^{\frac{5}{2}}} \\
\int_{0}^{\pi} \frac{d x}{(a-\cos x)^{4}} &=-\frac{\pi}{6} \cdot \frac{3 a\left(2 a^{2}+3\right)}{\left(a^{2}-1\right)^{\frac{7}{2}}} =-\frac{\pi a\left(2 a^{2}+3\right)}{2\left(a^{2}-1\right)^{\frac{7}{2}}}
\end{aligned}
$$
Theoretically, we can proceed to find $I_n(a)$ for any $n\in N$ by the recurrence relation in $(2)$ .
By Mathematical Induction, we can further prove that the formula is $$
\int_{0}^{\pi} \frac{d x}{(a-\cos x)^{n}}=\frac{\pi P(a)}{\left(a^{2}-1\right)^{\frac{2 n-1}{2}}}
$$
for some polynomial $P(a)$ of degree $n-1$.
Last but not least, how to find the formula for $P(a)$?  Would you help me?

One way using complex analysis:
Since the integrand is even:
$$
I=\int_{0}^{\pi} \frac{d x}{(a-\cos x)^n} = \frac{1}{2} \int_{-\pi}^{\pi} \frac{d x}{(a-\cos x)^n}
$$
Since $a>1$, we can expand the integrand with the generalized binomial theorem:
$$ I = \frac{1}{2} \int_{-\pi}^{\pi} \frac{d x}{(a-\cos x)^n} = \frac{1}{2} \int_{-\pi}^{\pi} \sum_{j=0}^{\infty} \binom{-n}{j} a^{-n-j}(-1)^j\cos^j(x) dx $$
From Fubini-Tonelli we can interchange integral and series:
$$ I = \frac{1}{2} \sum_{j=0}^{\infty} \binom{-n}{j}a^{-n-j}(-1)^j\int_{-\pi}^{\pi}  \cos^j(x) dx $$
Do the following substitution
$$ \cos x = \frac{z+z^{-1}}{2}$$
$$ dx = \frac{dz}{zi}$$
The integral is now a contour integral round the unit complex circle:
$$ I = \frac{1}{2} \sum_{j=0}^{\infty} \binom{-n}{j}\frac{a^{-n-j}(-1)^j}{2^j}\oint_{|z|=1}  \frac{(z^2+1)^j}{z^{j+1} }dz$$
The integral has a pole at $z=0$. To find the residue, expand the integrand:
$$\frac{(z^2+1)^j}{z^{j+1} }= \sum_{k=0}^{j}\binom{j}{k} z^{2k-j-1}$$
The residue is the coefficient of $z^{2k-j-1}= z^{-1}$
Then
$$ 2k-j-1 = -1 \Longrightarrow k = \frac{j}{2}$$
Therefore, the residue exists if $j$ is $\textbf{divisible by}$ $2$:
$$\oint_{|z|=1}  \frac{(z^2+1)^j}{z^{j+1} }dx = 2\pi i \operatorname{Res}\left(\frac{(z^2+1)^j}{z^{j+1}},0\right) = 2\pi i \binom{j}{\frac{j}{2}}  $$
Hence, we have
$$ I = \pi\sum_{j=0}^{\infty} \binom{-n}{2j}\binom{2j}{j}\frac{a^{-n-2j}}{2^{2j}}  $$
Note
$$\binom{-n}{2j} = \frac{(-n-2j+1)_{2j}}{(2j!}$$
where $(x)_{n} = x(x+1)\cdots(x+n-1)$ is the rising factorial (Pochhammer polynomial)
and
$$\binom{2j}{j} = \frac{(2j)!}{j!^2}$$
Hence
$$\binom{-n}{2j}\binom{2j}{j} = \frac{(-n-2j+1)_{2j}}{j!^2}$$
From the reflection formula for the Pochhammer polynomial:
$$(-x)_{m} = (-1)^m(x-m+1)_{m}$$
we have
$$(-n-2j+1)_{2j}= (n)_{2j}$$
and from the duplication formula for the degree of the Pochhammer polynomial:
$$(x)_{2m} = 4^m\left(\frac{x}{2}\right)\left(\frac{1+x}{2}\right)$$
we have
$$(n)_{2j} = 4^{j}\left(\frac{n}{2}\right)\left(\frac{n+1}{2}\right)$$
Hence
$$ I = \pi\sum_{j=0}^{\infty} \binom{-n}{2j}\binom{2j}{j}\frac{a^{-n-2j}}{2^{2j}} = \frac{\pi}{a^n}\sum_{j=0}^{\infty}\frac{\left(\frac{n}{2}\right)\left(\frac{n+1}{2}\right)}{(1)_{j}} \frac{\left(\frac{1}{a^2}\right)^j}{j!}$$
$$\Longrightarrow I = \frac{\pi}{a^n}{}_{2}F_{1}\left({\frac{n}{2},\frac{n+1}{2}\atop 1};\frac{1}{a^2}\right)$$
where ${}_{2}F_{1}({a,b\atop c};z)$ is the Guassian function
Now using the following forumula:
$${}_{2}F_{1}\left({a,a+\tfrac{1}{2}\atop c};z\right)=2^{c-1}z^{\frac{(1-c)}{2}}%
(1-z)^{-a+\left(\frac{c-1}{2}\right)}P^{1-c}_{2a-c}\left(\frac{1}{\sqrt{1-z}}\right) $$
where $P_{n}^{(\alpha,\beta)}$ is the Jacobi polynomial
we have
$$ I =\frac{\pi}{a^n}{}_{2}F_{1}\left({\frac{n}{2},\frac{n+1}{2}\atop 1};\frac{1}{a^2}\right)=  \frac{\pi}{a^n}\left(1-\frac{1}{a^2}\right)^{-\frac{n}{2}}P_{n-1}\left(\frac{a}{\sqrt{a^2-1}}\right)$$
where $P_{n}$ is the standard Legendre polynomial which satisfies the Rodrigues' formula:
$$ P_{n}(x) = \frac{1}{2^n n!} \frac{d^n}{dx^n} (x^2-1)^n$$
Finally, we can conclude
$$\boxed{\int_{0}^{\pi} \frac{d x}{(a-\cos x)^n} = \frac{\pi}{2^{n-1}(n-1)!a^n}\left(1-\frac{1}{a^2}\right)^{-\frac{n}{2}}\left( \frac{d^{n-1}}{dx^{n-1}} (x^2-1)^{n-1} \right)_{x = \frac{a}{\sqrt{a^2-1}}}}$$
This is similar to the formula you found. However, Legendre polynomials also have an explicit formula:
$$P_{n} (x) = \sum_{j} (-1)^{\frac{j}{2}} \frac{(2n-j-1)!!}{j!!(n-j)!} x^{n-j} \quad j=0,2,4,...,n-\frac{1}{2}\pm \frac{1}{2}$$
(ends with $n$ in the case $n$ even and ends with $n-1$ in the case $n$ odd).
Hence
$$ \boxed{I =\int_{0}^{\pi} \frac{d x}{(a-\cos x)^n} =  \frac{\pi}{a^n}\left(1-\frac{1}{a^2}\right)^{-\frac{n}{2}}\sum_{j} (-1)^{\frac{j}{2}} \frac{(2n-j-3)!!}{j!!(n-1-j)!} \left(\frac{a}{\sqrt{a^2-1}}\right)^{n-1-j} \quad j=0,2,4,...,n-\frac{3}{2}\pm \frac{1}{2} \textrm{ ($n-1$ even or odd)}}   $$
It turned out that the polynomial you are looking for is this Legendre polynomial.
