The solid angle can be expressed in terms of complete elliptic integrals of first and third kind.
I believe this is first derived by F. Paxton in 1959.
THE REVIEW OF SCIENTIFIC INSTRUMENTS. VOLUME 30, NUMBER 4 APRIL, 1959.
Solid Angle Calculation for a Circular Disk. F. PAXTON
A copy of the article can be found
here.
Below is my attempt to derive the same result in an alternate manner.
The integral is indeed trickier than I thought.
To avoid confusion with the usage of $x,y,z$ as coordinates, let us change the problem and assume the circular disk $C$ we are dealing with lies in the plane $z = b$, centered at $(a,0,b)$ with radius $s$. We will assume $a, b > 0$ and $a \ne s$.
In terms of ordinary spherical polar coordinates $(r,\theta,\phi)$, points on the plane $z = b$ can be parametrized as:
$$(x,y,z) = (b\tan\theta \cos\phi,\,b\tan\theta\sin\phi,\, b)$$
Let $\rho = \sqrt{x^2+y^2} = b\tan\theta$, the "element" for integrating the solid angle is given by:
$$d\Omega
= \sin\theta d\theta \wedge d\phi
=\frac{\tan\theta d\tan\theta}{\sqrt{1 + \tan\theta^2}^3} \wedge d\phi
= \frac{b \rho d\rho \wedge d\phi}{\sqrt{\rho^2+b^2}^3}
= \frac{b\;dx \wedge dy}{\sqrt{\rho^2+b^2}^3}
$$
Introduce complex coorindates $\eta = x + iy$ and $\bar{\eta} = x - iy$, the solid angle
extended by $C$ can be rewritten as:
$$\begin{align}
\Omega_{C}
= & \int_{C} d\Omega
= \frac{b}{2i}\int_{C} \frac{d\bar{\eta} \wedge d\eta}{\sqrt{\eta\bar{\eta}+b^2}^3}
=\color{blue}{^{[1]}} \frac{b}{i} \int_C d\left(\frac{1}{\eta}\left( \frac{1}{b} - \frac{1}{\sqrt{\eta\bar{\eta}+b^2}}\right)\right) \wedge d\eta\\
= & \frac{b}{i} \int_{\partial C} \left(\frac{1}{b} - \frac{1}{\sqrt{\eta\bar{\eta} + b^2}}\right) \frac{d\eta}{\eta}
= 2\pi \delta_{C} + ib \int_{\partial C} \frac{d\eta}{\eta \sqrt{\eta\bar{\eta} + b^2}}
\end{align}$$
where $\delta_{C} = 1 \text{ or } 0$ depends on whether $s > a$ or $< a$. In other words, whether $\partial C$ contains $0$ or not.
On $\partial C$, we can parametrize $\eta$ as $s e^{it} + a$ and $\bar{\eta}$ as $s e^{-it} + a$. Substitute this into above expression, we obtain:
$$\begin{align}
&\Omega_C - 2\pi \delta_{C}\\
= &ib \int_{-\pi}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}\frac{ is e^{it}}{s e^{it} + a }\\
= & -b
\int_{0}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}
\left( \frac{s e^{it}}{s e^{it} + a } + \frac{s e^{-it}}{s e^{-it} + a } \right)\\
= & -b
\int_{0}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}
\frac{2s(s + a\cos t)}{s^2 + a^2 + 2sa\cos t}\\
= & -b
\int_{0}^{\pi} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos t}}
\left( 1 + \frac{s^2 - a^2}{s^2 + a^2 + 2sa\cos t} \right)\\
= & -2b
\int_{0}^{\pi/2} \frac{ dt }{\sqrt{ s^2 + a^2 + b^2 + 2sa\cos(2t)}}
\left( 1 + \frac{s^2 - a^2}{s^2 + a^2 + 2sa\cos(2t)} \right)\\
= & -2b
\int_{0}^{\pi/2} \frac{ dt }{\sqrt{ (s+a)^2 + b^2 - 4sa\sin^2(t)}}
\left( 1 + \frac{s^2 - a^2}{(s+a)^2 - 4sa\sin^2(t)} \right)
\end{align}$$
The last integral can be expressed in terms of the
complete elliptic integral
of first and third kind:
$$
K(k) = \int_{0}^{\frac{\pi}{2}} \frac{dt}{\sqrt{1-k^2\sin^{2}t}}
\quad\text{ and }\quad
\Pi(n,k) = \int_{0}^{\frac{\pi}{2}} \frac{dt}{(1-n\sin^{2}t)\sqrt{1-k^2\sin^{2}t}}
$$
The final results are:
$$\Omega_{C} = 2\pi\delta_C -\frac{2b}{\sqrt{(s+a)^2+b^2}}
\left(\;K(k) + \left(\frac{s-a}{s+a}\right) \Pi(n,k)\;\right)
$$
where $\displaystyle \;\;k = \sqrt{\frac{4sa}{(s+a)^2+b^2}}
\quad\text{ and }\quad
n = \frac{4sa}{(s+a)^2}
$.
Notes
$\color{blue}{[1]}$ When $C$ contains $(0,0,b)$, we need the $\frac{1}{b}$ term. It regularize the 1-form at $\eta = 0$ and make Stroke's theorem continue to work.
As per Biot-Savart's law, magnetic field due to a small element is:
$$dB=\frac{\mu_0}{4\pi}\frac{i\vec{dl}\times \hat{r}}{r^2}$$
Notice that in polar coordinates, $\vec{dl}=dr\hat{r}+r\,d\theta \hat{\theta}$ i.e
$$\vec{dl}\times \hat{r}=-r\,d\theta \hat{z}$$
where $\hat{z}$ is the direction pointing out of the plane.
$$\Rightarrow \vec{B}=\frac{\mu_0 i}{4\pi} \int_{-\pi/2}^{\pi/2} \frac{-r\,d\theta}{r^2}\,\hat{z}$$
$$\begin{aligned}
\Rightarrow \left|\vec{B}\right| &=\frac{\mu_0 i}{4\pi} \int_{-\pi/2}^{\pi/2} \frac{d\theta}{r}=\frac{\mu_0 i}{4a\pi} \int_{-\pi/2}^{\pi/2} \cos\theta(1+\sin^2\theta)\,d\theta\\
&=\frac{\mu_0 i}{4a\pi}\int_{-1}^1 (1+t^2)\,dt\,\,\,\,(\sin\theta=t) \\
&=\boxed{\dfrac{2}{3}\dfrac{\mu_0 i}{a\pi}} \\
\end{aligned}$$
Best Answer
Start by rescaling the length dimensions by $l$, i.e. let $x=l\,\tilde x$ and $y=l\,\tilde y$. The fact that this transformation simplifies the integral we have to evaluate would be reason enough to go ahead and do this if you were just interested in this integral from a purely mathematical standpoint, but given the physical interpretation of the integral at hand the result of this simplification becomes a very interesting result in its ownright:
$$I(l) = \int_{0}^{l}\int_{0}^{l}\frac{k\,q\,l}{(x^2+y^2+l^2)^{\frac{3}{2}}}\,dxdy\\ = \int_{0}^{1}\int_{0}^{1}\frac{k\,q}{(\tilde x^2+\tilde y^2+1)^{\frac{3}{2}}}\, d\tilde x d\tilde y\\ = I(1).$$
In other words, the size of the cube's side length $l$ doesn't matter at all; the integral will have the same final value no matter what the value of $l$ happens to be, so we might as well set it equal to unity: $l=1$.
Now we evalauate $I=k\,q\int_{0}^{1}\int_{0}^{1}\frac{1}{(x^2+y^2+1)^{\frac{3}{2}}}\, dx dy$. To evaluate the integral $\int_{0}^{1}\frac{1}{(x^2+y^2+1)^{\frac{3}{2}}}\,dy$, use the trigonometric substitution,
$$y=\sqrt{x^2+1}\tan{u}.$$
Then,
$$\int_{0}^{1}\frac{1}{(x^2+y^2+1)^{\frac{3}{2}}}\,dy = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{x^2+1}}}}\frac{\sqrt{x^2+1}\sec^2{u}}{(x^2+1+(x^2+1)\tan^2{u})^{\frac{3}{2}}}\,du\\ = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{x^2+1}}}}\frac{\sqrt{x^2+1}\sec^2{u}}{(x^2+1)^{\frac32}(1+\tan^2{u})^{\frac{3}{2}}}\,du\\ = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{x^2+1}}}}\frac{\sec^2{u}}{(x^2+1)\sec^3{u}}\,du\\ = \frac{1}{x^2+1}\int_{0}^{\tan^{-1}{\frac{1}{\sqrt{x^2+1}}}}\cos{u}\,du\\ = \frac{1}{x^2+1}\frac{1}{\sqrt{x^2+2}}.$$
Next we need to integrate $\frac{I}{kq} = \int_{0}^{1}\frac{1}{(x^2+1)\sqrt{x^2+2}}\,dx$. Trigonometric substitution again does the trick. Letting $x=\sqrt{2}\tan{u}$,
$$\int_{0}^{1}\frac{1}{(x^2+1)\sqrt{x^2+2}}\,dx = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{2}}}}\frac{\sqrt{2}\sec^2{u}}{(2\tan^2{u}+1)\sqrt{2\tan^2{u}+2}}\,du\\ = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{2}}}}\frac{\sec^2{u}}{(2\tan^2{u}+1)\sqrt{\tan^2{u}+1}}\,du\\ = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{2}}}}\frac{\sec{u}}{2\tan^2{u}+1}\,du\\ = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{2}}}}\frac{\cos{u}}{2\sin^2{u}+\cos^2{u}}\,du\\ = \int_{0}^{\tan^{-1}{\frac{1}{\sqrt{2}}}}\frac{\cos{u}}{\sin^2{u}+1}\,du\\ = \int_{0}^{1/\sqrt{3}}\frac{dw}{w^2+1}=\tan^{-1}{\frac{1}{\sqrt{3}}} = \frac{\pi}{6},$$
where in the last line we have substituted $w=\sin{u}$. Thus, the final result is:
$$I = \frac{\pi kq}{6}.~~~\blacksquare$$