In classical physics, the sign of a potential (say, gravitational potential or electric potential) is merely a convention. For example, since the electric field $\mathbb{E}$ is conservative, it is the gradient of some potential, which we call the electric potential or Coulomb potential $V$, so it can be written as $\mathbb{E}=\nabla V$ or $\mathbb{E}=-\nabla V$. Mathematically, it doesn’t matter whichever you use. If you choose $\mathbb{E}=\nabla V$, Coulomb potential should be $V=-\frac{1}{4\pi\epsilon_0}\frac{Q}{r}$ and if you choose $\mathbb{E}=-\nabla V$, $V=\frac{1}{4\pi\epsilon_0}\frac{Q}{r}$. In physics, the usual convention is $\mathbb{E}=-\nabla V$ and it is interpreted as electric field pointing downhill towards lower voltages.

On the other hand, in quantum mechanics the sign of Coulomb potential matters and it does depend on the problem that you are studying. In a hydrogen atom or a hydrogen-like atom (for example, an ionized helium atom $\mathrm{He}^+$) an electron is trapped in the atom by Coulomb force between proton and electron. In this case, we say the electron is in a bound state. The bound state with the minimum energy is called the ground state and the minimum energy is called the ground state energy. A bound state with energy higher than the ground state energy is said to be unstable and the electron in an unstable bound state always tend to move to the ground state by emitting photons. Modeling bound state of a hydrogen atom or a hydrogen-like atom requires negative Coulomb potential $V(r)=-\frac{Z\alpha}{r}$ as seen in the following figure.

In this simple figure, the electron is trapped in the region $0<r<1$. The figure also clearly shows you that bound state energy must be negative. By drawing a picture, one can easily see that bound state cannot be modeled with positive Coulomb potential, but positive Coulomb potential is used to model scattering of a particle.

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|$$