Let $n\geq 0$ be an integer and let $J_n=J_n(r)$ denote the usual Bessel function (of the first kind) of order $n$ i.e. one of the solutions to Bessel's differential equation

$$r^2\frac{d^2y}{dr^2}+r\frac{dy}{dr}+(r^2-n^2)y=0.$$

I am interested in the quantities

$$\beta_n:=\int_0^\infty J_n(r)^2 J_0(r)^2 dr$$
which arise naturally in the study of certain classes of oscillatory integrals. In particular, I would like to explore the convexity of the sequence $\{\beta_n\}$. *My question is the following:* is it true that

$$\beta_n<\frac{\beta_{n-1}+\beta_{n+1}}{2}$$ for every $n\geq 1$? This inequality seems to be numerically verifiable (Mathematica) for small values of $n$. To prove it, I was trying to use the so-called ''Neumann's formula'' (see e.g. Watson's treatise (1966), p.32)

$$J_n(r)^2=\frac{1}{\pi}\int_0^\pi J_{2n}(2r\sin\theta)d\theta,$$ together with the well-known recurrence

$$J_{n-1}(r)-J_{n+1}(r)=2 \frac{d J_n}{dr}(r)$$ and some integration by parts. Unfortunately I wasn't able to make this approach work.

Ideas or/and references are welcome, thank you very much in advance.