This note is based on my friend Khin Maung’s short lecture.

Let us begin with Schrödinger equation
$$\left(\frac{\hat p^2}{2m}+\hat V(r)\right)|\psi\rangle=E|\psi\rangle$$
Use the completeness relation
$$1=\int|\vec{p}\rangle\langle\vec{p}|d\vec{p}$$
to get the momentum space representation of the Schrödinger equation
\begin{equation}
\label{eq:schrodingerms}
\frac{p^2}{2m}\psi(\vec{p})+\int\langle\vec{p}|\hat V(r)|\vec{p’}\rangle\psi(\vec{p’})d\vec{p’}=E\psi(\vec{p})
\end{equation}
where $\psi(\vec{p}):=\langle\vec{p}|\psi\rangle$. \eqref{eq:schrodingerms} is called the Schrödinger equation in momentum space. Using the completeness relation
$$1=\int|\vec{r}\rangle\langle\vec{r}|d\vec{r}$$
we obtain
$$\langle\vec{p}|\hat V(r)|\vec{p’}\rangle=\frac{1}{(2\pi\hbar)^3}\int e^{\frac{i}{\hbar}(\vec{p’}-\vec{p})\cdot\vec{r}}V(r)d\vec{r}$$
Here, recall that $\langle\vec{p}|\vec{r}\rangle=\frac{1}{(2\pi\hbar)^{\frac{3}{2}}}e^{-\frac{i}{\hbar}\vec{p}\cdot\vec{r}}$. Let $\vec{q}=\vec{p’}-\vec{p}$ and $V(\vec{q})=\langle\vec{p}|\hat V(r)|\vec{p’}\rangle$. Then we have
$$V(\vec{q})=\frac{1}{(2\pi\hbar)^3}\int e^{\frac{i}{\hbar}\vec{q}\cdot\vec{r}}V(r)d\vec{r}$$
This is just the Fourier transform of $V(r)$. For the Yukawa potential
$$V(r)=V_0\frac{e^{-\mu r}}{r}$$
\begin{align*} V(\vec{q})&=\frac{1}{(2\pi\hbar)^3}\int e^{\frac{i}{\hbar}qr\cos\theta}\int_0^{2\pi}\int_0^\pi\int_0^\infty V(r)r^2dr\sin\theta d\theta d\phi\\ &=\frac{1}{(2\pi)^2\hbar^3}\int_0^\infty V(r)r^2\int_{-1}^1e^{\frac{i}{\hbar}qru}dudr\\ &=\frac{1}{(2\pi\hbar)^2iq}\int_0^\infty V(r)r(e^{\frac{i}{\hbar}qr}-e^{-\frac{i}{\hbar}qr})dr\\ &=\frac{V_0}{(2\pi\hbar)^2iq}\int_0^\infty e^{-\mu r}(e^{\frac{i}{\hbar}qr}-e^{-\frac{i}{\hbar}qr})dr\\ &=\frac{V_0}{2\pi^2\hbar^3}\frac{1}{\mu^2+\frac{q^2}{\hbar^2}} \end{align*}
From here on, we assume that $\hbar=1$ for simplicity.

Let
$$\psi(\vec{p})=\psi_l(p)Y_l^m(\hat p)$$
Here, $\hat p$ stands for the unit vector in the momentum space in spherical coordinates $\hat p=\frac{\vec{p}}{p}=(\theta_l,\phi_p)$ where $\theta_l$ is the polar angle and $\phi_p$ is the azimuth angle corresponding to the momentum vector $\vec{p}$. $V(\vec{q})$ can be written as the series
\begin{align*} V(\vec{q})&=\sum_{l=0}^\infty\sum_{m=-l}^lV_l(p,p’)Y_l^m(\hat p)Y_l^{m}(\hat p’)\\
&=\sum_{l=0}^\infty\frac{2l+1}{4\pi}V_l(p,p’)P_l(x)
\end{align*} The last line is obtained by the addition theorem for spherical harmonics $$\frac{4\pi}{2l+1}\sum_{m=-l}^lY_l^m(\hat p)Y_l^{m}(\hat p’)=P_l(x)$$
where $x=\cos\theta_{pp’}$.
By the orthonormality of spherical harmonics
$$\int_0^{2\pi}\int_0^\pi Y_{l_1}^{m_1}(\hat p)Y_{l_2}^{m_2}(\hat p)\sin\theta d\theta d\phi=\delta_{l_1l_2}\delta_{m_1m_2}$$
the Schrödinger equation in momentum space \eqref{eq:schrodingerms} yields
$$\frac{p^2}{2m}\psi_l(p)+\int_0^\infty V_l(p,p’)\psi_l(p’)p’^2dp’=E\psi_l(p)$$
Using the orthogonality of Legendre polynomials
$$\int_{-1}^1P_{l’}(x)P_l(x)dx=\frac{2}{2l+1}\delta_{ll’}$$
we obtain
$$V_l(p,p’)=2\pi\int_{-1}^1V(\vec{q})P_l(x)dx$$
For the Yukawa potential, we have
$$V(\vec{q})=\frac{V_0}{2\pi^2}\frac{1}{\mu^2+q^2}=\frac{V_0}{2\pi^2}\frac{1}{p^2+p’^2+\mu^2-2pp’x}$$
The Neumann’s formula yields then
$$V_l(p,p’)=\frac{V_0}{pp’\pi}Q_l\left(\frac{p^2+p’^2+\mu^2}{2pp’}\right)$$
Note here that we require $\left|\frac{p^2+p’^2+\mu^2}{2pp’}\right|>1$. Recall that
$$Q_0(z)=\frac{1}{2}\ln\frac{z+1}{z-1}$$
for $|z|>1$. Hence for the Yukawa potential, $V_0(p,p’)$ can be written as
$$V_0(p,p’)=\frac{V_0}{2pp’\pi}\ln\frac{(p+p’)^2+\mu^2}{(p-p’)^2+\mu^2}$$
With $V_0=Z\alpha$ and $\mu=0$, the Yukawa potential reduces to the Coulomb potential $V(r)=\frac{Z\alpha}{r}$. $V_l(p,p’)$ is then given by$$V_l(p,p’)=\frac{Z\alpha}{pp’\pi}Q_l\left(\frac{p^2+p’^2}{2pp’}\right)$$
and accordingly,$$V_0(p,p’)=\frac{Z\alpha}{pp’\pi}\ln\left|\frac{p+p’}{p-p’}\right|$$