With hindsight the identity is not that magical: the Bessel functions play only a secondary role, in the sense that there is a more general identity for arbitrary differentiable functions $f : [0,1] \to \mathbb{R}$, namely
\begin{align}
&12 \int_0^{1/2}\mathrm{d}x\,\cot \pi x \int_0^x \mathrm{d}y\, \cot \pi y \, f(y)\big[f(x)\,f(1-x-y)-f(1-x)\,f(x-y)\big] \\
&= -f(1)\,f(0)^2+ 2 \int_0^1\mathrm{d}x\int_0^{1-x} \mathrm{d}y\, f(x)f(y)f(1-x-y).\tag{1}
\end{align}
Substituting $f(x) = J_0(t\,x)$ and using the Laplace transform $\mathcal{L}\{J_0(t\, x)\}(\xi) = (t^2 + \xi^2)^{-1/2}$, we find that the double convolution integral satisfies
$$\int_0^1\mathrm{d}x\int_0^{1-x} \mathrm{d}y\, J_0(t\,x)J_0(t\,y)J_0(t(1-x-y)) = -\frac{1}{t}J_0'(t).$$
Together with $J_2(t) = -J_0(t)-2 J_0'(t)/t$ and $J_0(0)=1$, this reproduces the claimed identity
$$J_2(t) = 12 \int_0^{1/2}\mathrm{d}x\,\cot \pi x \int_0^x \mathrm{d}y\, \cot \pi y \, J_0(ty)\big[J_0(tx)\,J_0(t(1-x-y))-J_0(t(1-x))\,J_0(t(x-y))\big].$$
Let's see why (1) holds. For this it is convenient to introduce the integral
$$I(\epsilon):=2\int_\epsilon^{1-2\epsilon}\mathrm{d}x\int_\epsilon^{1-x-\epsilon} \mathrm{d}y\, f(x)f(y)f(1-x-y)$$
with $0<\epsilon <1/4$.
The meat of the argument is Hermite's cotangent identity,
$$\cot(\pi x)\, \cot(\pi y) + \cot(\pi y)\,\cot(\pi(1-x-y)) + \cot(\pi(1-x-y))\cot(\pi x)= 1,$$
which after insertion in the integrand and using the symmetry under permutation of $x$, $y$ and $1-x-y$, implies
$$I(\epsilon) = 6\int_\epsilon^{1-2\epsilon}\mathrm{d}x\int_\epsilon^{1-x-\epsilon} \mathrm{d}y\,\cot(\pi x)\cot(\pi y)\, f(x)f(y)f(1-x-y).$$
Since the integrand is symmetric in $x$ and $y$ we can further reduce the domain to
$$I(\epsilon) = 12\int_\epsilon^{1-2\epsilon}\mathrm{d}x\int_\epsilon^{\min(1-x-\epsilon,x)}\!\!\!\! \mathrm{d}y\,\cot(\pi x)\cot(\pi y)\, f(x)f(y)f(1-x-y).$$
Folding the region $1/2< x < 1$ onto $0 < x < 1/2$ via the substitution $x\to 1-x$ gives
\begin{align}
I(\epsilon) =& 12\int_{2\epsilon}^{1/2}\mathrm{d}x\int_\epsilon^{x-\epsilon}\! \mathrm{d}y\,\cot(\pi x)\cot(\pi y)\, f(y)\big[f(x)f(1-x-y)-f(1-x)f(x-y)\big]\\
&+ 12 \int_\epsilon^{1/2}\mathrm{d}x\int_{\max(x-\epsilon,\epsilon)}^{\min(1-x-\epsilon,x)}\! \mathrm{d}y\,\cot(\pi x)\cot(\pi y)\, f(x)f(y)f(1-x-y).
\end{align}
Both $I(\epsilon)$ and the first term on the right-hand side have limits given by the corresponding terms in (1) as $\epsilon \to 0$.
It remains to examine the last integral on the right-hand side. By the change of variables $x = \epsilon\, \hat{x}$ and $y = \epsilon\,\hat{y}$ we see that it has a limit as $\epsilon\to 0$ given by
\begin{align}
&\frac{12 f(1)f(0)^2}{\pi^2}\int_1^\infty \mathrm{d}x\int_{\max(x-1,1)}^x \mathrm{d}y \frac{1}{x\,y} \\
&=\frac{12 f(1)f(0)^2}{\pi^2}\int_1^\infty \mathrm{d}y\int_{y}^{y+1} \mathrm{d}x \frac{1}{x\,y} \\
&= \frac{12 f(1)f(0)^2}{\pi^2} \int_{1}^\infty \frac{1}{y^2}\log\left(\frac{y+1}{y}\right) = f(1)f(0)^2.
\end{align}
This finishes the proof of (1).