Integral ${\large\int}_0^\infty\big(2J_0(2x)^2+2Y_0(2x)^2-J_0(x)^2-Y_0(x)^2\big)\,dx$

I'm interested in the following definite integral: $$\int_0^\infty\big(2J_0(2x)^2-J_0(x)^2+2Y_0(2x)^2-Y_0(x)^2\big)\,dx,\tag1$$ where $J_\nu$ and $Y_\nu$ are the Bessel functions of the first and the second kind.

Mathematica evaluates this integral symbolically to $\frac{\ln2}\pi$, but the result of a numerical integration looks more like $\frac{\ln4}\pi$, so I suspect the symbolic result is incorrect.

Moreover, it looks like both components converge and make equal contribution: $$\int_0^\infty\big(2J_0(2x)^2-J_0(x)^2\big)\,dx\stackrel?=\int_0^\infty\big(2Y_0(2x)^2-Y_0(x)^2\big)\,dx\stackrel?=\frac{\ln2}\pi,\tag2$$ but it is more difficult to check numerically, because both integrands here are oscillating (unlike $(1)$ where the integrand looks monotonic).

How can we find values of these integrals and prove them correct? Can we generalize results for values of the index $\nu$ other than $0$?

I stated an interesting thing in the comments (wisely explored by Random Variable), but it is definitely not the fastest way to go.

For the glory of Italian mathematicians, Frullani's theorem is the way (again):

$$ \int_{0}^{+\infty}\frac{x\, J_0(x)^2-2x\, J_0(2x)^2}{x}\,dx = (-\log 2)\cdot\underset{x\to+\infty}{\widetilde{\lim}} x\cdot J_0(x)^2 = \color{red}{-\frac{\log 2}{\pi}} $$ and the argument is just the same for the $Y_0$ function.

The wide tilde-limit has to be intended as a Cesàro mean: $$\underset{x\to+\infty}{\widetilde{\lim}}f(x) = \lim_{x\to +\infty}\frac{1}{x}\int_{x}^{2x}f(z)\,dz.$$

Moreover, provided that the integral converges, the previous result holds unchanged also by replacing $J_0$ with $J_\nu$ for any $\nu\geq 0$, since: $$ J_\nu(z) = \sqrt{\frac{2}{\pi z}}\left(\cos \left(z-\frac{\nu\pi}{2}-\frac{\pi}{4}\right)+O\left(\frac{1}{z}\right)\right) $$ for any large $z\in\mathbb{R}^+$.

Somewhat similar to this answer, we can use Ramanujan's master theorem to show that $$\int_{0}^{\infty}\big(2J_0(2x)^2-J_0(x)^2\big)\,dx = \frac{\ln 2}{\pi}. $$

As $x \to \infty$, $$\big(2J_0(2x)^2-J_0(x)^2\big) \sim \frac{\sin (4x)-\sin (2x)}{\pi x}. \tag{1}$$ So the integral does indeed converge (but not absolutely).

The hypergeometric representation of the square of the Bessel function of the first kind of order zero is $$J_{0}(z)^{2} = \, _1F_2\left(\frac{1}{2}; 1, 1; -z^{2} \right). $$

So for $a>0$ and $0 < s < 1$, the Mellin transform of $J_{0}(ax)^{2}$ is $$ \begin{align} \int_{0}^{\infty} x^{s-1} J_{0}(ax)^{2} \,dx &= \int_{0}^{\infty} x^{s-1} \, _1F_2\left(\frac{1}{2}; 1, 1; -(ax)^{2} \right) \, dx \\ &= \frac{1}{2a^{s}}\int_{0}^{\infty} u^{s/2-1} \, _1F_2\left(\frac{1}{2}; 1, 1; -u \right) \, du \\ &= \frac{1}{2a^{s}} \Gamma\left(\frac{s}{2} \right) \frac{\Gamma \left(\frac{1}{2}-\frac{s}{2} \right) \Gamma(1)^{2} }{\Gamma\left(\frac{1}{2}\right) \Gamma\left(1- \frac{s}{2} \right)^{2}} \\ &= \frac{1}{2a^{s}} \Gamma\left(\frac{s}{2} \right) \frac{\Gamma \left(\frac{1}{2}-\frac{s}{2} \right)}{\sqrt{\pi} \, \Gamma\left(1- \frac{s}{2} \right)^{2}} . \end{align}$$

Therefore, assuming we can bring the limit inside the integral, $$ \begin{align} \int_0^\infty\big(2J_0(2x)^2-J_0(x)^2\big)\,dx &= \frac{1}{\pi} \lim_{s \to 1^{-}} \left(\frac{1}{2^{s}}- \frac{1}{2} \right) \Gamma \left(\frac{1}{2} - \frac{s}{2} \right) \\ &= \frac{1}{\pi} \lim_{s \to 0^{-}} \left( \frac{1}{2^{s+1}} - \frac{1}{2} \right)\Gamma \left(- \frac{s}{2} \right) \\ &= \frac{1}{\pi} \lim_{s \to 0^{-}} \left(- \frac{\ln 2}{2} s + \mathcal{O}(s^{2}) \right) \left(- \frac{2}{s} + \mathcal{O}(1) \right) \\ &= \frac{\ln 2}{\pi}. \end{align}$$

The same approach seemingly also shows that $$\int_{0}^{\infty}\big(2J_v(2x)^2-J_v(x)^2\big)\,dx = \frac{\ln 2}{\pi} \, , \quad v> - \frac{1}{2}. $$

$(1)$ https://en.wikipedia.org/wiki/Bessel_function#Asymptotic_forms