Double integral over modified Bessel function of the second kind of zeroth order

Question

I am working on a project we I arrive at the integral
\begin{align}
\int_0^a \int_0^a K_0(|x-y|) \, \mathrm{d}x \mathrm{d}y \,,
\end{align}
where $K_0(\cdot)$ is the modified Bessel function of the second kind of zeroth order. I wonder if this integral can be evaluated analytically to arrive at an expression depending on the upper bound $a$, or if it is only possible to evaluate this integral numerically for a given value of $a$. I found an analytical expression for $a \rightarrow \infty$ on https://dlmf.nist.gov/10.43, which is
\begin{align}
\int_{0}^{\infty}K_{\nu}\left(t\right)\mathrm{d}t=\tfrac{1}{2}\pi\sec\left(%
\tfrac{1}{2}\pi\nu\right),
\end{align}
but this is only a single integral and in my case the upper bound should not be infinity. I was not able to come up with an answer to this, my guess would be that it is only possible to evaluate it numerically. But maybe someone else knows more about this, all ideas are very much appreciated.

Answer 1

As commented, the integral expresses in terms of modified Struve functions. Using $$\int_0^a\!\!\!\int_0^a f(|x-y|)\,dx\,dy=2\int_0^a(a-z)f(z)\,dz$$ (consider the double integral as iterated, substitute $x=y+z$ and exchange the integrations, keeping track of the possible values of $y$ with $z\in(-a,a)$ fixed) and these formulae, one gets $$\int_0^a\!\!\!\int_0^a K_0(|x-y|)\,dx\,dy=\pi a^2\big(K_0(a)\mathbf{L}_{-1}(a)+K_1(a)\mathbf{L}_0(a)\big)-2\big(1-aK_1(a)\big).$$