An integral involving Airy functions $\int_0^\infty\frac{x^p}{\operatorname{Ai}^2 x + \operatorname{Bi}^2 x}\mathrm dx$

I need your help with this integral: $$\mathcal{K}(p)=\int_0^\infty\frac{x^p}{\operatorname{Ai}^2 x + \operatorname{Bi}^2 x}\mathrm dx,$$ where $\operatorname{Ai}$, $\operatorname{Bi}$ are Airy functions: $$\operatorname{Ai}\,x=\frac{1}{\pi}\int_0^\infty\cos\left(x\,z+\frac{z^3}{3}\right)\,\mathrm dz,$$ $$\operatorname{Bi}\,x=\frac{1}{\pi}\int_0^\infty\left(\sin\left(x\,z+\frac{z^3}{3}\right)+\exp\left(x\,z-\frac{z^3}{3}\right)\right)\,\mathrm dz.$$ I am not sure that $\mathcal{K}(p)$ has a general closed form, but I hope so, because approximate numerical calculations suggest these conjectured values: $$\mathcal{K}(3)\stackrel?=\frac{5\,\pi^2}{32},\ \ \mathcal{K}(6)\stackrel?=\frac{565\,\pi^2}{512}.$$

Let us start with a warm-up exercise. Introduce the functions $$g_{\pm}(x)=\operatorname{Ai}(x)\pm i\operatorname{Bi}(x).$$ Computing the Wronskian of these two solutions of the Airy equation, one can check that $$\frac{1}{\operatorname{Ai}^2(x)+\operatorname{Bi}^2(x)}=\frac{\pi}{2i}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]$$ This gives the integral $\mathcal{K}(0)$ as $$\mathcal{K}(0)=\pi\left[\arg g_+(\infty)-\arg g_+(0)\right]=\pi\left[\pi-\frac{\pi}{3}\right]=\frac{\pi^2}{6}.$$

To compute the integral $\mathcal{K}(3n)$, we will need to develop a more sophisticated approach. First note that (see here) $$g_{\pm}(x)=-2e^{\mp 2\pi i/3}\operatorname{Ai}\left(e^{\mp2\pi i/3}x\right).$$ Therefore \begin{align}\mathcal{K}(3n)&=\frac{\pi}{2i}\int_0^{\infty}x^{3n}\left[\frac{g_+'(x)}{g_+(x)}-\frac{g'_-(x)}{g_-(x)}\right]dx=\\ &=\frac{\pi}{2i}\lim_{R\rightarrow\infty}\int_{S_R}z^{3n}\frac{\operatorname{Ai}'(z)}{\operatorname{Ai}(z)}dz, \end{align} where the contour $S_R$ in the complex $z$-plane is composed of two segments: one going from $Re^{2\pi i/3}$ to $ 0$ and another one going from $0$ to $ Re^{-2\pi i/3}$.

It is a well-known fact that the Airy function $\operatorname{Ai}(z)$ has zeros (i.e. our integrand has poles) on the negative real axis only. Therefore by residue theorem our integral is equal to $$\mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{R\rightarrow \infty}\int_{C_R}z^{3n}\left[\ln\operatorname{Ai}(z)\right]'dz,\tag{1}$$ where $C_R$ is the arc of the circle of radius $R$ centered at the origin going counterclockwise from $Re^{-2\pi i/3}$ to $Re^{2\pi i/3}$.

The limit (1), on the other hand, can be computed using the asymptotics of the Airy function as $z\rightarrow\infty$ for $|\arg z|<\pi$: \begin{align} \ln\operatorname{Ai}(z)\sim -\frac23 z^{3/2}-\ln2\sqrt{\pi}-\frac14\ln z+ \ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}\left(\frac43 z^{3/2}\right)^{-k}.\tag{2} \end{align} Note that if we introduce instead of $z$ the variable $s=\frac43z^{3/2}$, then the integration will be done over the circle of radius $\Lambda=\frac43 R^{3/2}$, i.e. a closed contour in the complex $s$-plane. The corresponding integral can therefore be computed by residues by picking the necessary term in the large $s$ expansion of $\ln \operatorname{Ai}(z)$.

More precisely, we have the following formula: \begin{align} \mathcal{K}(3n)=-\frac{\pi}{2i}\lim_{\Lambda\rightarrow \infty}\oint_{|s|=\Lambda} \left(\frac{3s}{4}\right)^{2n}d\left[-\frac16\ln s+\ln\sum_{k=0}^{\infty}\frac{(-1)^k\left(\frac16\right)_k\left(\frac56\right)_k}{k!}s^{-k}\right]\tag{3} \end{align} To compute the residue, it suffices to expand the logarithm-sum up to order $2n$ in $s^{-1}$. Note that the Pochhammer symbol coefficients are in fact some rational numbers.

In the simplest case $n=0$, the residue is determined by the term $-\frac16\ln s$ and we readily reproduce the previous result $$\mathcal{K}(0)=-\frac{\pi}{2i}\cdot 2\pi i\cdot\left(-\frac16\right)=\frac{\pi^2}{6}.$$ The general formula for arbitrary $n$ would look a bit complicated (but straightforward to obtain) due to the need to expand the logarithm of the sum.

Example. The calculation of the corresponding values $M(n)=\mathcal{K}(3n)$ in Mathematica can be done using the command

\begin{align}\mathtt{\text{ M[n_] := -Pi^2 SeriesCoefficient[ Series[(3 s/4)^(2 n) D[-Log[s]/6 +}} \\ \mathtt{\text{ Log[Sum[(-1)^k Pochhammer[1/6, k] Pochhammer[5/6, k]/(k! s^k), }} \\ \mathtt{\text{{k, 0, 2 n}]], s], {s, Infinity, 1}], 1]}} \end{align}

This yields, for instance, $$M(0)=\frac{\pi^2}{6},\quad M(1)=\frac{5\pi^2}{32},\quad M(2)=\frac{565\pi^2}{512},$$ $$\ldots, M(10)=\frac{2\,660\,774\,144\,147\,177\,521\,025\,228\,125\,\pi^2}{2\,199\,023\,255\,552},\ldots$$ and so on.

Added: The large $s$ expansion (2) can also be found directly by using Airy equation. Moreover, by transforming it into an equation for $\ln \operatorname{Ai}(z)$, one can avoid reexpanding the logarithm of the sum. The price to pay will be that the expansion coefficients will be determined by a nonlinear recurrence relation instead of explicit formulas.

I can give a partial answer. Using DLMF, Chapter 9 "Airy and Related Functions", Airy Functions:

• §9.2 (ii) Initial Values, formulae 9.2.3, 9.2.5,
• §9.8 (i) Modulus and Phase, Definitions, formulae 9.8.3, 9.8.4, and
• §9.8 (ii) Modulus and Phase, Identities formula 9.8.14,

we can see that $$\frac{d}{dx}\arctan\frac{\operatorname{Ai}x}{\operatorname{Bi}x}=-\frac1\pi\frac1{\operatorname{Ai}^2x+\operatorname{Bi}^2x}$$ and $$\int_0^\infty\frac{dx}{\operatorname{Ai}^2x+\operatorname{Bi}^2x}=\frac{\pi^2}6.$$

This settles the question for the particular case $p=0$.