How to prove that $\int_{0}^{\infty} \frac{1}{\operatorname{Ai}(x)^{2}+\operatorname{Bi}(x)^{2}} dx = \frac{\pi^{2} }{6} $?

Question

In omegadot's answer to this question, the user asserts that $$\int_{0}^{\infty} \frac{1}{\operatorname{Ai}(x)^{2}+\operatorname{Bi}(x)^{2}} dx = \frac{\pi^{2} }{6} .$$ Here, $\operatorname{Ai}(x)$ and $\operatorname{Bi}(x)$ are the Airy functions of the first and second kind, respectively.

Though numerical evidence suggests it's true, no proof is provided in the answer itself. I can't find any sources on this evaluation, and I wouldn't know where to start myself.

Question: how to prove the integral identity above?

Bonus questions:

1. Are there any articles that go more deeply into these types of integrals?
2. Is anything known about integrals of the class $I_{n} = \int_{0}^{\infty} \Big{(} \frac{1}{\operatorname{Ai}(x)^{2}+\operatorname{Bi}(x)^{2}} \Big{)}^{n} dx $ ?

Answer 1

For 1, we can note that the integral follows from the lack of a first derivative term in the differential equation, meaning the Wronskian of the differential equation is constant. Thus, for any differential equation of the form $y''(x) + b(x) y(x) = 0$ and any two linearly independent solutions $y_1$, $y_2$, we have $$ \frac{d}{dx}\left[\frac{y_1(x)}{y_2(x)}\right] = \frac{y_2(x)y_1'(x) - y_1(x)y_2'(x)}{y_2(x)^2} = \frac{C}{y_2(x)^2}, $$ and thus for any quadratic in $y_1$ and $y_2$, we have \begin{multline} \int \frac{dx}{ay_1(x)^2 + 2by_1(x)y_2(x) + c y_2(x)^2} = \int\frac{1}{a[y_1(x)/y_2(x)]^2 + 2b [y_1(x)/y_2(x)] + c}\frac{dx}{y_2(x)^2}\\= C\int\frac{du}{a u^2 + 2b u + c}\;\;,\;\; u = y_1(x)/y_2(x)\;\;,\;\; C\in\mathbb C \\ =C\tan^{-1}\left[\frac{ay_1(x)/y_2(x) + b}{\sqrt{ac-b^2}}\right] = C\ln\left[\frac{y_1(x)-r_1y_2(x)}{y_1(x)+r_2y_2(x)}\right] \end{multline} where $r_1$ and $r_2$ are the roots of the quadratic. Note that the $C$ in the above equations is not necessarily the same number each time. It represents some constant that would depend on the form of the antiderivative and the normalization of $y_1$ and $y_2$.

For 2, things aren't so easy. The trick above relied on the quadratic nature of the polynomial in the denominator. I suspect such integrals do not have a general form, though one never knows with these sorts of things.