My question is technically just a math question that came up in the context of a physics problem. In particular, it came up in a perturbation calculation for a two dimensional electron gas. Long story short, the integral is $$\int_{k'<k_F}\int_{k<k_F}\frac{d^2\mathbf{k}\,d^2\mathbf{k}'}{|\mathbf{k}-\mathbf{k}'|}$$ and I am not sure how to evaluate this. For clarity, these integrals are over the 2D Fermi Sphere (i.e., the Fermi disk) with radius $k_F$.
For those curious of the physical context, the situation that this turned up in is a two-dimensional electron gas (not ignoring spin) with a $1/r$ two-body interaction. Writing the potential $V(\mathbf{r}_i-\mathbf{r}_j) = A/|\mathbf{r}_i-\mathbf{r}_j|$, it turns out that the interaction operator can be written in second quantized form as $$ \frac{1}{2}\sum_{i\neq j}V(\mathbf{r}_i-\mathbf{r}_j) = \frac{1}{2}\sum_{\substack{\mathbf{k},\mathbf{k}',\mathbf{q} \\ s,s'}}V_\mathbf{q} c_{\mathbf{k}-\mathbf{q},s}^\dagger c_{\mathbf{k}'+\mathbf{q},s'}^\dagger c_{\mathbf{k'},s'}c_{\mathbf{k},s} $$ where the $c$ and $c^\dagger$'s are the annihilation and creation operators for the state indicated in their subscript. The sum is over allowed momenta and spins. The matrix elements are $$ V_\mathbf{q} = \frac{1}{\Omega}\int_{\mathbb{R}^2} V(\mathbf{r}) e^{i\mathbf{q}\cdot\mathbf{r}}\,d^2\mathbf{r} = \frac{2\pi A}{\Omega q}. $$ where $\Omega$ is the area of the box containing the gas.
Ignoring the direct scattering term $\mathbf{q}=0$ (it is being compensated in a Jellium-like manner), we wish to evaluate the exchange ($\mathbf{q}=\mathbf{k}-\mathbf{k}'$) energy shift to the ground state $$ E_\text{ex} = \frac{\pi A}{\Omega}\sum \frac{1}{|\mathbf{k}-\mathbf{k}'|} $$ where this sum proceeds only over the occupied states with same spin. Taking the continuum limit leads to the integral in question.
There are more details regarding the spin and an actual net magnetization, but they are irrelevant to the integral I need to solve.
I was originally going to post this on Physics Stackexchange, but I figured since it boils down to a math question, someone here might be able to help.
Thanks in advance for any advice you may be able to offer!
Let $k = (r_1\cos\theta_1,r_1\sin\theta_1), k' = (r_2\cos\theta_2,r_2\sin\theta_2)$.
In physicist's convention, the integral we want to evaluate has the form:
$$\mathcal{I} \stackrel{def}{=} \int_0^{k_F} r_1 d r_1 \int_0^{k_F}r_2 dr_2 \int_0^{2\pi} d\theta_1 \int_0^{2\pi} d\theta_2 \frac{1}{ \sqrt{r_1^2+r_2^2 - 2r_1r_2\cos(\theta_2-\theta_1)}} $$ Let $r = \max(r_1,r_2)$, $\displaystyle\;\lambda = \frac{\min(r_1,r_2)}{\max(r_1,r_2)}$, $\phi = \theta_2 - \theta_1 + \pi$. Let $x = \lambda\cos\phi$, $y = \lambda\sin\phi$.
Change variable from $r_1, r_2, \theta_1, \theta_2$ first to $r, \lambda, \theta_1, \phi$ and then to $r, \theta_1, x,y$, we have
$$\begin{align} \mathcal{I} &= 2\int_0^{k_F} r dr \int_0^1 r^2\lambda d\lambda \int_0^{2\pi}d\theta_1\int_0^{2\pi}d\phi\frac{1}{r\sqrt{1+\lambda^2+2\lambda\cos\phi}}\\ &= 2\left(\int_0^{k_F}r^2dr\right) \left(\int_0^{2\pi} d\theta_1\right)\int_0^1 \lambda d\lambda \int_0^{2\pi}d\phi\frac{1}{\sqrt{(\lambda\cos\phi + 1)^2 + (\lambda\sin\phi)^2}}\\ &= \frac{4\pi}{3}k_F^3\int_{x^2+y^2\le 1}\frac{dxdy}{\sqrt{(x+1)^2+y^2}} \end{align} $$ Introduce "polar" coordinate $(\rho,\psi)$ centered at $(-1,0)$, i.e, $(x,y) = (-1+\rho\cos\psi,\rho\sin\psi)$.
We can evaluate the integral in last line as
$$\int_{x^2+y^2\le 1}\frac{dxdy}{\sqrt{(x+1)^2+y^2}} = \int_{-\pi/2}^{\pi/2}\int_0^{2\cos\psi} \frac{\rho d\rho d\psi}{\rho} = \int_{-\pi/2}^{\pi/2} 2\cos\psi d\psi = 4$$ As a result, $\displaystyle\;\mathcal{I} = \frac{16\pi}{3} k_F^3$.