The derivative of an expression involving Bessel functions of the first kind An integral that comes up when discussing the orthogonality relations of the Bessel function of the first kind is $$ \int x J_p(\alpha x ) J_p(\beta x) \, \mathrm dx = \frac{x\left(\alpha J'_p(ax)J_p(\beta x) - \beta J_p(ax) J'_p(\beta x) \right)}{\beta^2 - \alpha^2} + C $$
The proof, which uses the defining differential equation of Bessel function of the first kind,  is quite simple.
I would like to confirm that the derivative of the right side of the equation is indeed $$x J_p(\alpha x ) J_p(\beta x)$$
I tried using identities like $$J_{p-1}(x) - J_{p+1}(x) = 2 J'_p(x)$$ and  $$J_{p-1}(x) + J_{p+1}(x) = \frac{2p}{x} J_p(x), $$ but I feel that I'm just going in circles.
 A: This s not actually anything to do with the relations between Bessel functions with different $p$: it arises purely from the fact that Bessel functions solve Bessel's equation. This can be written in Sturm–Liouville form as
$$ -(xu')' + \frac{p^2}{x}u=x u. $$
Differentiating the right-hand side of the given expression using the product rule judiciously gives
\begin{align} 
&(\alpha x J_p'(\alpha x)J_p(\beta x)-\beta x J_p'(\beta x)J_p(\alpha x))' \\
&= \alpha( x J_p'(\alpha x))'J_p(\beta x) + \alpha x J_p'(\alpha x) (J_p(\beta x))' - \beta (x J_p'(\beta x))' J_p(\alpha x)-\beta x J_p'(\beta x)(J_p(\alpha x))' \\
&= \alpha (x J_p'(\alpha x))'J_p(\beta x) - \beta (x J_p'(\beta x))' J_p(\alpha x) + \alpha\beta x J_p'(\alpha x) J_p'(\beta x) -\alpha\beta x J_p'(\beta x)J_p'(\alpha x) \\
&= ( x J_p(\alpha x)')'J_p(\beta x) - (x J_p(\beta x)')' J_p(\alpha x)
\end{align}
Now, if we change variables in the Sturm–Liouville equation to $x=ay$, then $ u' = \frac{1}{a} \frac{du}{dy}  $, so $u(ay)$ satisfies
$$ -\frac{1}{a}\frac{d}{dy}\left(y\frac{du}{dy}\right) + \frac{p^2}{ay}u =ay u, $$
or
$$ -\frac{d}{dy}\left(y\frac{du}{dy}\right) + \frac{p^2}{y}u=a^2y u, \tag{*} $$
which is also a Sturm–Liouville equation; in particular, $J_p(ay)$ satisfies this. Substituting in, we find
$$ ( x J_p(\alpha x)')'J_p(\beta x) - (x J_p(\beta x)')' J_p(\alpha x) = \left( \frac{p^2}{x}-\alpha^2 x \right) J_p'(\alpha x))'J_p(\beta x) - \left( \frac{p^2}{x}-\beta^2 x \right) J_p(\alpha x) = (\beta^2-\alpha^2)x J_p(\alpha x)J_p(\beta x), $$
as required.

Exactly the same approach works for any eigenfunctions of a Sturm–Liouville equation: if $u_1$ and $u_2$ satisfy $-(pu_1')'+qu_1 = \lambda_1 u_1$ and $-(pu_2')'+qu_2 = \lambda_2 u_2$, then
$$ \int wu_1u_2 = \frac{p (u_2'u_1-u_1'u_2)}{\lambda_1-\lambda_2} + C, \tag{$\dagger$} $$
since the derivative of the right-hand side's numerator is
$$ (pu_2')'u_1+pu_2'u_1'-(pu_1')'u_2 - pu_1'u_2' = (q-\lambda_2 w)u_2u_1 - (q-\lambda_1 w)u_1u_2 = (\lambda_1-\lambda_2)wu_1u_2.  $$

The slightly awkward thing here is that Bessel functions can be treated as solutions to the Sturm–Liouville equation (*) with two different eigenvalues: either $1/x$ or $x$ can be treated as the weight, with $p^2$ or $a^2$ being the eigenvalue. The latter is solved by ${\scr I}_p(ax)$ for fixed $p$ (where ${\scr I}$ denotes any solution to (*)), with $a$ taking values in the discrete set chosen so that the boundary conditions work), while the former is solved by varying $p$ and keeping $a$ fixed. The latter is also unusual in that $a$ is accounted for by simply scaling the argument, which is why the question's expression looks a bit different from the general version ($\dagger$) given above.
