Let $n\ge 1$ be an integer, let
$$F(x,y)=\int_0^\infty u^{n(x+y)} (K_{x-y}(u))^n du$$
for $x,y\ge 0$.
When $n=1$, this is just Mellin transform of the Bessel K function. When $n=2$, $F(x,y)$ has an explicit form in product of Gamma functions, given by the Parseval formula for Mellin transform.
For general $n$, I expect some Stirling formula type estimation for $F(x,y)$. I tried with Laplace method but didn't get anywhere.