Numerical Solutions of the Non-Relativistic Yukawa Theory

The Yukawa potential (aka. screened Coulomb potential) $V(r) = -(\alpha/r) \exp(-\mu r)$ is a universal potential for interactions mediated by massive particles in the non-relativistic limit. It is also useful for describing interactions in the medium -- the "screened Coulomb interaction". Thus it is very important to find the solutions for the Yukawa potential in non-relativistic Schroedinger equation: $$\Big[ -\frac{\nabla^2}{2m} - \frac{\alpha}{r} e^{-\mu r} \Big] \psi(\vec r) = E \, \psi(\vec r),$$ where $m$ is the reduced mass, $\alpha \equiv g^2/(4\pi)$ is the strength of the interaction with $g$ the "charge", $\mu$ is the mass of the exchanged particle, and $\mu=0$ gives the Coulomb potential. Alas, the above eigenvalue problem does not admit analytic solution, except for the Coulomb case. This post describes a numerical method for solving Schroedinger equation with the Yukawa potential.

Fig. 1: Feynman diagram of a boson (wavy line) mediated interaction.

In the momentum space, the eigenvalue equation becomes an integral equation: $$\frac{\vec p^2}{2m} \psi(\vec p) - \int \frac{\mathrm{d}^3 p'}{(2\pi)^3} \frac{4\pi \alpha}{(\vec p-\vec p')^2+\mu^2} \psi(\vec p') = E\,\psi(\vec p).$$ Here $\psi(\vec p)$ represents the momentum-space wave functions, related to the coordinate-space wave function $\psi(\vec r)$ by a Fourier transformation, $$\psi(\vec p) = \int \mathrm d^3r \, e^{i \vec p\cdot \vec r} \psi(\vec r). \\ \psi(\vec r) = \int \frac{\mathrm d^3 p}{(2\pi)^3} \, e^{-i \vec p\cdot r} \psi(\vec p).$$ We have abused the notation $\psi$ for wave functions. The wave functions should be normalized according to $$\int \mathrm d^3 r\, \psi^*(\vec r) \psi(\vec r) = 1, \\ \int \frac{\mathrm d^3 p}{(2\pi)^3} \psi^*(\vec p) \psi(\vec p) = 1.$$

Dimension Reduction

Before employing numerical methods, let us first put the equation to a dimensionless form. The natural scale in the problem is the Bohr radius ($\hbar = c = 1$): $a \equiv 1/(\alpha m)$. We also define: $\lambda \equiv a\mu$, $\kappa \equiv 2 m E a^2 = 2 E/ (m \alpha^2) \equiv E / E_B$, with $E_B \equiv \frac{1}{2} \alpha^2 m$ being the ground state energy (binding energy) of the Coulomb theory (such as the Hydrogen atom, positronium and so on), and $\vec p \to a \vec p$, $\vec r \to \vec r/a$. Then, the equations becomes, $$\Big[ - \nabla^2 - \frac{2}{r} e^{-\lambda r} \Big] \psi(\vec r) = \kappa \, \psi(\vec r), \\ \vec p^2 \psi(\vec p) - \int \frac{\mathrm d^3 p'}{(2\pi)^3} \frac{8\pi}{(\vec p - \vec p')^2 + \lambda^2} \psi(\vec p') = \kappa \, \psi(\vec p).$$ The theory has a rotational symmetry and the wave functions can be written as, $$\psi(\vec r) = R(r) Y_{lm}(\hat r), \quad \psi(\vec p) = \phi(p) Y_{lm}(\hat p).$$ $Y_{lm}$ are the spherical harmonics. They are the eigen function of the angular momentum squared operator $\vec L^2$. Therefore, we only need to solve for the radial part $R(r)$ or $\phi(p)$. They are related by a Hankel transformation: $$R(r) = \frac{(-i)^\ell}{8\pi^3} \int_0^\infty \mathrm dp \, p^2 j_\ell(pr) \phi(p), \\ \phi(p) = i^\ell \int_0^\infty \mathrm dr \, r^2 j_\ell(pr) R(r).$$ $j_\ell(z)$ is the spherical Bessel function of the first kind. $Y_{lm}$ are normalized: $$\int \mathrm d^2\Omega(\hat r) Y^*_{l'm'}(\hat r) Y_{lm}(\hat r) = \delta_{ll'}\delta_{mm'}.$$ Similarly, the radial wave function are normalized according to: $$\int_0^\infty \mathrm d r \, r^2 R^2(r) = 1, \\ \frac{1}{8\pi^3} \int_0^\infty \mathrm d p\, p^2 \phi^2(p) = 1.$$

Coordinate Space

The coordinate-space Schroedinger equation can be written in the spherical coordinates: $$\Big[ - \frac{1}{r^2}\frac{\partial }{\partial r}\big( r^2 \frac{\partial }{\partial r} ) + \frac{l(l+1)}{r^2} - \frac{2}{r} e^{-\lambda r} \Big] R(r) = \kappa \, R(r). \\ \Leftrightarrow \Big[- \frac{\partial^2 }{\partial r^2} + \frac{l(l+1)}{r^2} - \frac{2}{r} e^{-\lambda r} \Big] u(r) = \kappa \, u(r),$$ where $u(r) = r R(r)$.

The above differential equations can be discretized and solved numerically. In this post, however, we will focus on the momentum space representation and its numerical solutions.

Momentum Space

In the momentum space, Schroedinger equation can be written as $$p^2 \phi(p) - \frac{4}{\pi} \int_0^\infty \mathrm dp' \, p'^2 K_\ell (p, p') \phi(p') = \kappa \, \phi(p),$$ where $$K_\ell(p, p') \equiv \int_0^\infty \mathrm dr \, r \exp(-\lambda r) j_\ell(pr) j_\ell(p' r).$$
Here are the analytical expressions for the first few $K_\ell$: $$\begin{split} (S\text{-wave})\quad K_0(p, p') =& \frac{1}{4 pp'}\log \Big[\frac{(p+p')^2+\lambda^2}{(p-p')^2+\lambda^2}\Big]; \\ (P\text{-wave})\quad K_1(p, p') =& \frac{p^2+p'^2+\lambda^2}{4 p^2p'^2}\log\Big[\frac{(p+p')^2+\lambda^2}{(p-p')^2+\lambda^2}\Big]-\frac{1}{2pp'}\\ (D\text{-wave})\quad K_2(p, p') =& \frac{3p^4+2p^2p'^2+3p'^4+6(p^2+p'^2)\lambda^2+3\lambda^4}{16p^3p'^3}\log\Big[\frac{(p+p')^2+\lambda^2}{(p-p')^2+\lambda^2}\Big] - \frac{3(p^2+p'^2+\lambda^2)}{8p^2p'^2}\\ \vdots& \end{split}$$ In the general case, $K_\ell$ can be evaluated numerically using quadrature method.

We can use the quadrature method to approximate the integral. For a general integral, $$\int \mathrm dx \, f(x) = \sum_{i=1}^N w_i f(x_i) + R[f^{(2N+1)}(\xi)]$$ where $x_i$ and $w_i$ are pre-chosen abscissas and weights. For the radial integral, we can use the Gauss-Legendre quadrature with a mapping function: $$\theta: (-1, +1) \to (0, \infty); \\ r_i = \theta(x^\text{gl}_i), \quad w_i = w^\text{gl}_i \theta'(x_i)$$ where $x^\text{gl}_i$, and $w^\text{gl}_i$ are the standard Gauss-Legendre quadrature abscissas and weights.

While in principle the converged result is independent of the mapping function, the choice of the mapping function is very important in practice, as it controls the rate of convergence. It should be chosen to cover the extent of the wave function. Because different states have different radial extent, it is often difficult to come up with a universal mapping function optimal for all excited states. Fortunately, in most practical cases, we are interested in the first lowest states.  Common choices include (first, $x$ is mapped to $(0,1)$ ): $r = x/(1-x)$, $r=a \tan (x \pi/2)$, $r = a (1-e^{-bx})/(e-e^x)$ ...

The integral equation becomes, $$p^2_i \phi(p_i) - \frac{4}{\pi} \sum_j w_j p^2_j K_\ell (p_i, p_j) \phi(p_j) = \kappa \, \phi(p_i).$$ It is convenient to define $v_i \equiv \sqrt{w_i} p_i \phi(p_i)$, with normalization $$\sum_i v_i^2 = 1.$$ Then, the discretized equation becomes, $$\sum_j H_{ij} v_j = \kappa \, v_i,$$ where the Hamiltonian matrix $H_{ij}$ reads, $$H_{ij} = \delta_{ij} p^2_i - (4/\pi) \sqrt{w_i w_j} \, p_i p_j \, K_\ell (p_i, p_j).$$ $H$ is obviously Hermitian. The energy eigenvalue and the wave functions can be obtained by numerically diagonalizing the Hamiltonian matrix.

Fig 2: Visualization of the Hamiltonian matrix, $\ell = 0$ (S-wave), quadrature order $N=512$.

The system can be solved in different quadrature order. The continuum limit can be reach by extrapolating the quadrature order $N\to\infty$.

Energy Eigenvalues 

The obtained energy eigenvalues as a function of the exchanged particle mass $\mu$ are shown in Fig. 3. At small $\mu$, the binding energies approach to the standard Coulomb values $E_n^{(\text{coul})}=-\alpha^2 m / (2n^2)$ where $n=1,2,3,4,\cdots$ is the principle quantum number. The ground state energy is also in good agreement with the upper bound obtained from the variational method, using the Coulomb wave function as trial functions. As $\mu$ increases, bound states gradually disappear. Even the ground state dissociates as $\mu \gtrsim \alpha m$.

Fig. 3: The S-wave energy eigenvalues as a function of $\mu$. The values are obtained by extrapolating results over quadrature grids of the order $N=16, 32, 64, 128, 512, 1024$. The red crosses represent Coulomb energy levels $E_n^{(\text{coul})} = - \alpha^2 m/(2n^2)$, $n=1,2,3,\cdots$. The red dashed line represents the upper bound obtained from variational method, using the Coulomb wave function as trial functions.

Wave Functions

Fig. 4: The obtained ground state S-wave momentum-space wave functions for selected $\mu$.

The momentum-space wave functions can be extracted from the obtained eigenvectors (see Fig. 4 and 5). At small $\mu$, the wave functions approach to the Coulomb result. As $\mu$ increases, the wave function within the momentum space becomes narrower. Furthermore, the change in excited wave functions are more dramatic comparing to the ground state wave function as $\mu$ increases.

Fig. 5: The obtained radially excited S-wave momentum-space wave functions for selected $\mu$.

.