Analytical evaluation of double integral
HINTS:
Exploiting the even symmetry and making use of
$$\int \sqrt{x^2+y^2}\,dx=\frac x2\sqrt{x^2+y^2}+\frac{x^2}{2} \log(x+\sqrt{x^2+y^2})+C$$
evaluation of the inner integral yields
$$\begin{align} \int_{-a}^a\int_{-a}^a \frac{\sqrt{x^2+y^2}}{4a^2}\,dx\,dy &= \frac1{a^2} \int_0^a \int_0^a \sqrt{x^2+y^2}\,dx\,dy \\\\ &=\frac1{a^2}\int_0^a\left(\frac a2\sqrt{a^2+y^2}+\frac{a^2}{2} \log(a+\sqrt{a^2+y^2})-\frac{a^2}{2}\log(y)\right) \,dy \tag 1 \end{align}$$
Use integration by parts for the middle term in the integral on the right-hand side of $(1)$
EDIT:
The OP has requested to see a presentation on evaluation of $\displaystyle \int_0^a \log(a+\sqrt{y^2+a^2})\,dy$ by integration of parts. To that end, we proceed.
Let $u=\log(a+\sqrt{y^2+a^2})$ and $v=y$. Then, we can write
$$\begin{align} \int_0^a \log(a+\sqrt{y^2+a^2})\,dy&=a\log(a(1+\sqrt{2}))-\int_0^a \frac{y^2}{a+\sqrt{y^2+a^2}}\frac{1}{\sqrt{y^2+a^2}}\,dy\\\\ &=a\log(a(1+\sqrt{2}))-\int_0^a \frac{\sqrt{y^2+a^2}-a}{\sqrt{y^2+a^2}}\,dy\\\\ &=a\log(a(1+\sqrt{2}))-a+a\int_0^a \frac{1}{\sqrt{y^2+a^2}}\,dy\\\\ &=a\log(a(1+\sqrt{2}))-a+a\log(a(1+\sqrt{2}))-a\log(a)\\\\ &=a(\log(a)-1+2\log(1+\sqrt{2})) \end{align}$$
Use polar coordinates (I obtain it by dividing each quarter of the square along the diagonal starting at the origin and noticing that the integral of the considered function is the same when calculated on each of these two domains) : $$I=\int^a_0\int^a_0 \frac{\sqrt{x^2+y^2}}{4a^2} \,dx \, dy=\frac{2}{4a^2} \int^{\frac\pi4}_0 \int^{a\sec \theta}_0 r^2 \, dr \, d\theta$$ So by symmetry : $$\int^a_{-a} \int^a_{-a}\frac{\sqrt{x^2+y^2}}{4a^2} \, dx \, dy = 4I = \frac2{a^2} \int^{\frac\pi4}_0\int^{a\sec \theta}_0 r^2 \, dr \, d\theta \\ $$$$= \frac2{a^2}\int^{\frac\pi4}_0 \frac{a^3}3\sec^3\theta \, d\theta = \frac {2a} 3 \int^{\frac\pi4}_0 \sec^3\theta \, d\theta\\ $$$$=\frac {2a}3\left[\frac{\sec(\theta)\tan(\theta)+\ln(\sec(\theta)+\tan(\theta)}{2}\right]_0^{\frac\pi4}\\=\frac{a}{3} (\sqrt{2} + \sinh^{-1}(1)).$$
$\newcommand{\bbx}[1]{\,\bbox[8px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
Note that $\ds{\left.\int_{-a}^{a}\int_{-a}^{a}{\root{x^{2} + y^{2}} \over 4a^{2}} \,\dd x\,\dd y\,\right\vert_{\ a\ >\ 0} = a\int_{0}^{1}\int_{0}^{1}\root{x^{2} + y^{2}}\,\dd x\,\dd y}$
With $\ds{\vec{r} = x\,\hat{x} + y\,\hat{y}}$, note that in $\ds{2}$D: \begin{align} r & = {1 \over 3}\nabla\cdot\pars{r\vec{r}\,} = {1 \over 3} \bracks{\partiald{\pars{rx}}{x} - \partiald{\pars{-ry}}{y}} = {1 \over 3}\,\bracks{\nabla \times \pars{-ry\,\hat{x} + rx\,\hat{y}\,}}_{\ z} \end{align}
With the Stokes Theorem:
\begin{align} &\int_{0}^{1}\int_{0}^{1}\root{x^{2} + y^{2}}\,\dd x\,\dd y = {1 \over 3}\oint_{\mrm{square}}\pars{-ry\,\hat{x} + rx\,\hat{y}\,}\cdot\dd\vec{r} \\[5mm] = &\ {1 \over 3}\bracks{ \int_{0}^{1}\root{1 + y^{2}} \times 1\,\dd y + \int_{1}^{0}\pars{-\root{x^{2} + 1} \times 1}\,\dd x} = {2 \over 3}\int_{0}^{1}\root{1 + x^{2}}\,\dd x \\[5mm] = &\ \bbx{\ds{\root{2} + \,\mrm{arcsinh}\pars{1} \over 3}} \end{align}
The integral is evaluated with the change of variables $\ds{x = \sinh\pars{\theta}}$:
\begin{align} \int_{0}^{1}\root{1 + x^{2}}\,\dd x & = \int_{0}^{\mrm{arcsinh\pars{1}}}\cosh^{2}\pars{\theta}\,\dd\theta = \int_{0}^{\mrm{arcsinh\pars{1}}}{1 + \cosh\pars{2\theta} \over 2}\,\dd\theta \\[5mm] & = {1 \over 2}\,\mrm{arcsinh}\pars{1} + \bracks{{1 \over 2}\,\sinh\pars{\theta} \root{1 + \sinh^{2}\pars{\theta}}}_{\ 0}^{\ \,\mrm{arcsin}\pars{1}} \\[5mm] & = {\,\mrm{arcsinh}\pars{1} + \root{2} \over 2} \end{align}