In 1D, the orthogonality relation for zeroth-order Bessel functions of the first kind is
$$\int_0^1xJ_0(xu_{0n})J_0(xu_{0m})dx=0.5J_1^2(u_{0m})\delta_{mn}$$
where $u_{0n}$ is the $n$-th root of $J_0(x)$.
Consider now an azimuthally symmetric Bessel function, $J_0(\rho u_{0n})$, where $\rho=\sqrt{x^2+y^2}$. What is the equivalent orthogonality relation? I tried something like
$$\int_0^1\int_0^1\sqrt{x^2+y^2}J_0\left(\sqrt{x^2+y^2}u_{0n}\right)J_0\left(\sqrt{x^2+y^2}u_{0m}\right)dxdy=0.5J_1^2(u_{0m})\delta_{mn}$$
but this doesn't work. P.s. for various reasons, I do not wish to replace $\rho=\sqrt{x^2+y^2}$ in the above integral (later I need to replace one of the two Bessel functions with an arbitrary function to compute an overlap integral, but this is beyond the scope of this question).
When changing from polar to Cartesian coordinates, the integral changes like $$\int_0^{2\pi}d\theta\int_0^1\rho J_0(\rho u_{0n})J_0(\rho u_{0n})d\rho=\int_0^1\int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}J_0(\sqrt{x^2+y^2} u_{0n})J_0(\sqrt{x^2+y^2} u_{0n})dydx.$$ The key here is that you have to make sure you are integrating over the same region in both coordinate systems.
We then have $$\int_0^1\rho J_0(\rho u_{0n})J_0(\rho u_{0n})d\rho=\frac{1}{2\pi}\int_0^1\int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}}J_0(\sqrt{x^2+y^2} u_{0n})J_0(\sqrt{x^2+y^2} u_{0n})dydx.$$