Help with $\int_{0}^{\infty} \frac{\cos^k (x)}{a^2 + x^2} \, dx$

Question

I'm having trouble evaluating the following integral: $$\int_{0}^{\infty} \frac{\cos^k(x)}{a^2+x^2} \, dx \,\text{ where } k \in \mathbb{N} \text{ and } a >0$$ I've tried to follow this question which approaches tan instead. $$I = \frac{1}{2} \int_{-\infty}^{\infty} \frac{\cos^k(x)}{a^2 + x^2} \, dx$$

$$I = -i \int_{0}^{i \infty} \frac{\cosh^k(x)}{(a-x)(a+x)} \, dx$$ From here, I'm not sure how to progress. I am aware of the following expansions of $\cos^n$: $$\cos^{2n} (z) = 2^{-2n} \binom{2n}{n} + 2^{1-2n} \sum_{k=0}^{n-1} \binom{2n}{k} \cos(2(n-k)z)$$ $$\cos^{2n} (z) = 1 + 2^{1-2n} \sum_{k=0}^{n-1} \binom{2n}{k}( \cos(2(n-k)z)-1)$$ $$\cos^{2n+1} (z) = 2^{-2n} \sum_{k=0}^{n} \binom{2n+1}{k} \cos((2n-2k+1)z)$$ $$\cos^{n} (z) = 2^{-n} \binom{n}{\frac{n}{2}}(1- \, n \text{ mod } 2) + 2^{1-n} \sum_{k=0}^{ \lfloor{\frac{n-1}{2}\rfloor}} \binom{n}{k} \cos((n-2k)z)$$ Based on @Claude Leibovici response, it is possible to evaluate the following two equations: \begin{equation} \int_{0}^{\infty} \frac{\cos^{2n} (x)}{a^2+x^2} \, dx = \frac{\sqrt{\pi } \Gamma \left(n+\frac{1}{2}\right) \left(1-2 \, _2F_1\left(1,-n;n+1;-e^{2 a}\right)\right)}{2 a \Gamma (n+1)}+\frac{\pi 4^{-n} e^{-2 a n} \left(e^{2 a}+1\right)^{2 n}}{a} \end{equation} \begin{equation} \int_{0}^{\infty} \frac{\cos^{2n+1} (x)}{a^2+x^2} \, dx = \frac \pi a e^{-a (2 n+1)} (1-\tanh (a))^{-(2 n+1)}- \frac{\sqrt{\pi } e^a \Gamma \left(n+\frac{3}{2}\right)}{a \Gamma (n+2)}{}_2F_1\left(1,-n;n+2;-e^{2 a}\right) \end{equation}

Answer 1

At least, for odd values of the exponent, the expansion given by @aaaaaaaaabbbbbbbbbcccccc in comments is very useful since $$\int_0^\infty \frac {\cos \big[ (2 n+1-2k)x\big] }{x^2+a^2}\,dx=\frac \pi{2a}\,\exp\big[-(2n+1-2k)a\big]$$ Summing from $k=0$ to $k=n$ then gives for $$I_n=\int_{0}^{\infty} \frac{\cos^{2n+1} (x)}{x^2+a^2 } \, dx$$ $$I_n=\frac \pi a e^{-a (2 n+1)} (1-\tanh (a))^{-(2 n+1)}- \frac{\sqrt{\pi } \,e^a\, \Gamma \left(n+\frac{3}{2}\right) \, }{a \,\Gamma (n+2)}\,_2F_1\left(1,-n;n+2;-e^{2 a}\right)$$ where appears the gaussian hypergeometric function.

Now, $$\,_2F_1\left(1,-n;n+2;-e^{2 a}\right)=\frac{P_n(t)}{\binom{2 n+1}{n+1}} \qquad \text{where} \qquad t=e^{2a} $$ The first terms are $$\left( \begin{array}{cc} n & P_n(t) \\ 0 & 1 \\ 1 & t+3 \\ 2 & t^2+5 t+10 \\ 3 & t^3+7 t^2+21 t+35 \\ 4 & t^4+9 t^3+36 t^2+84 t+126 \\ 5 & t^5+11 t^4+55 t^3+165 t^2+330 t+462 \\ 6 & t^6+13 t^5+78 t^4+286 t^3+715 t^2+1287 t+1716 \\ 7 & t^7+15 t^6+105 t^5+455 t^4+1365 t^3+3003 t^2+5005 t+6435 \\ 8 & t^8+17 t^7+136 t^6+680 t^5+2380 t^4+6188 t^3+12376 t^2+19448 t+24310 \end{array} \right)$$