Scattering length and effective range of microscopic two-body potentials

Abstract

Scattering processes are a fundamental way of experimentally probing distributions and properties of systems in several areas of physics. Considering two-body scattering at low energies, when the de Broglie wavelength is larger than the range of the potential, partial waves with high angular momentum are typically unimportant. The dominant contribution comes from $l=0$ partial waves, commonly known as s-wave scattering. This situation is very relevant in atomic physics, e.g. cold atomic gases, and nuclear physics, e.g. nuclear structure and matter. This manuscript is intended as a pedagogical introduction to the topic while covering a numerical approach to compute the desired quantities. We introduce low-energy scattering with particular attention to the concepts of scattering length and effective range. These two quantities appear in the effective-range approximation, which universally describes low-energy processes. We outline a numerical procedure for calculating the scattering length and effective range of spherically symmetric two-body potentials. As examples, we apply the method to the spherical well, modified Pöschl-Teller, Gaussian, and Lennard-Jones potentials. We hope to provide the tools so students can implement similar calculations and extend them to other potentials.

Keywords:
Low-energy scattering; scattering length; effective range

1. Introduction

Scattering is commonly described as a process in which the observed object interacts with a scattering center. The objects are typically classical particles, waves or quantum particles. As to the scatterer, we generally deal with another entity, represented as a scattering potential. From a classical point of view, a collision is the most common example, where physical parameters such as the cross-section and the scattering angle arise [¹, ²]. From a quantum point of view, scattering is related to an incident wave that “collides” with a scattering potential, resulting in a scattered wave function.

In quantum scattering theory [¹, ³], if the de Broglie wavelength is comparable to (or larger than) the range of the scattering potential, we are at the low-energy limit, in which interesting behavior appears. The properties of low-energy scattering can be described universally by two parameters: the scattering length and the effective range [⁴]. Physically, the scattering length can be understood as an “effective size” of the target potential [³, ⁵]. As the name suggests, the effective range may be defined as the “real” range of the scattering potential [⁶].

Low-energy quantum scattering theory is important in many areas, such as nuclear, atomic, and condensed matter physics. In most physical systems, the interaction, and consequently the scattering length and effective range, is fixed. Hence, the task becomes constructing short-range potentials that describe the properties of the system. However, cold atomic gases provide a richer scenario from the point of view of varying the scattering length. The interatomic interactions between two atoms can be tuned in the laboratory via Feshbach resonances [⁷, ⁸], which allows us to explore how physical properties change as the interaction is varied. The universal behavior in the low-energy limit enables us to draw a parallel between cold atomic gases and nuclear systems, for example [⁹,¹⁰,¹¹].

Besides textbooks, several references aim to introduce aspects of scattering theory pedagogically. Detailed studies regarding analytical solutions in simple one-dimensional (1D) potentials are discussed in Refs. [¹²,¹³,¹⁴]. For 1D potentials that do not support bound states, quantities such as the reflection and transmission coefficients are typically calculated. In the tridimensional (3D) case, the quantum theory of scattering is often needed, where quantities such as the $S$-matrix, $T$-matrix, and Green’s functions receive attention. Concerning analytical solutions, Ref. [¹⁵] introduces electron-positron scattering by using a covariant form of the $S$-matrix. Reference [¹⁶] gives an analytical treatment to Green’s functions of delta potentials with a parallel with the scattering amplitude. Scattering length and effective range are discussed for an arbitrary angular momentum in Ref. [¹⁷], which also presents analytical calculations for the hard-sphere, soft-sphere, spherical well and a combination of the last two potentials. As to numerical solutions, Ref. [¹⁸] provides a numerical treatment to the 1D well potential by using Numerov’s method, while Ref. [¹⁹] uses the Runge-Kutta integration method to study the variable phase equation in 3D scattering.

This work aims to construct microscopic two-body potentials to describe interactions with a desired scattering length and effective range. Since we are tuning the potential to reproduce two values, the functional form of the potential must have at least two parameters. Typically, one describes the strength of the attraction, and another is responsible for the range.

Although there are known analytical results for low-energy scattering by some potentials, it is impossible to find closed expressions for arbitrary potentials. Hence, we employed a numerical approach. This work shows how to compute the scattering length and the effective range of spherically symmetric two-body potentials by numerically solving the Schrödinger equation. Using Numerov’s method [²⁰], we integrate the equation, and the scattering theory gives the boundary conditions. Once we have the wave function, calculating the scattering length and effective range is straightforward.

We devised the structure of this work to encompass readers with several distinct goals. The main text is aimed at readers who wish to cover how the scattering length and the effective range arise in low-energy scattering theory and how to compute them numerically, but without going into the formal scattering theory. For those who want to cover scattering in more detail, in Appendix ^A, we show how to obtain the Lippmann-Schwinger equation, employ Green’s functions, and contour integrations in the complex plane to get the asymptotic behavior of the wave function. Covering the scattering length in graduate and undergraduate quantum mechanics courses is common, but this is not true for the effective range. Readers familiar with the scattering length may start with Sec. ^2.3, where we introduce the effective range and the shape-independent approximation. However, if the reader is proficient in low-energy scattering and seeks a numerical approach to calculate the scattering length and effective range of a given potential, Sec. ³ onwards provides the necessary steps with relevant examples.

This work is structured as follows. In Sec. ², we briefly review quantum scattering theory fundamentals focusing on the low-energy limit. To keep this section concise, we start with a qualitative discussion of the asymptotic behavior of a plane wave scattered by a potential. For the readers interested in the details of scattering theory and the derivation of the equations of Sec. ², we provide Appendix ^A. Section ³ contains the numerical procedure to calculate the scattering length and the effective range from the zero-energy solutions of the Schrödinger equation. In Sec. ⁴, we apply the method to the spherical well, modified Pöschl-Teller (mPT), Gaussian, and Lennard-Jones (LJ) potentials as examples to illustrate the method. Finally, in Sec. ⁵, we present our final considerations.

2. Scattering theory

We briefly present some quantum scattering theory results to define the scattering length, following textbook references [¹, ³]. We are also interested in the effective-range expansion, which is not typically covered in undergraduate quantum mechanics courses. For this reason, we employ arguments from references [⁵, ⁶]. The reader is referred to the references mentioned above for a detailed derivation of the results.

Formally, a scattering process is described as a transition from one quantum state to another. We are interested in a particle, initially far away from the region where it will be scattered, moving toward the scattering center. The initial state $|i\mathrm{\rangle }$ is assumed to be a plane wave $|\mathbf{k}\mathrm{\rangle }$ since it is a free particle. The final state $|f\mathrm{\rangle }$ is the result of the action of a scattering potential on $|i\mathrm{\rangle }$ which, at large distances, is an outgoing spherical wave as illustrated in Fig. 1. The initial state $|i\mathrm{\rangle }$ satisfies the eigenvalue equation for the free-particle Hamiltonian $H_0={\mathbf{p}}^2/2m$,

Figure 1
Scheme of a scattering process. We chose our coordinate system such that the incoming particle has its momentum in the $z$-direction, $\mathbf{k}=k\widehat{z}$. A plane wave $e^{i\mathbf{k}\mathbf{\cdot }\mathbf{r}}=e^{ikz}$ travels in space until it encounters a scattering potential. It produces an outgoing spherical wave $e^{ikr}/r$ at large distances. The scattering angle $\theta$ defines the direction in which the momentum was transferred after the scattering process.

(1) $$H_0\text{|}i\mathrm{\rangle }=E_i\text{|}i\mathrm{\rangle }=\frac{{\hslash }^2{k}^2}{2m}\text{|}k\mathrm{\rangle }.$$

We choose our coordinate system such that the scattering center is at the origin, and the position vector is given by $\mathbf{r}=x\widehat{\mathbf{x}}+y\widehat{\mathbf{y}}+z\widehat{\mathbf{z}}$. Scattering is taken into account by introducing a potential $V(\mathbf{r})$ that acts on a finite range $r<R$, where $r\equiv |\mathbf{r}|$. The Hamiltonian in this region is

(2) $$H=H_0+V(\mathbf{r}).$$

Our goal is to understand the action of $H$ on free-particle states and how the scattering occurs. In other words, we want to know how the initial state $|i\mathrm{\rangle }$ transitions to a final state $|f\mathrm{\rangle }$.

The initial state is given by a plane wave, $\mathrm{|}i\mathrm{\rangle }=|\mathbf{k}\mathrm{\rangle }$. We consider the particle to be inside a cubic box of side $L$, such that the free particle is represented in the position basis as

(3) $$\mathrm{\langle }\mathrm{}\mathbf{r}|\mathbf{k}\mathrm{\rangle }=𝒩\mathrm{}{e}^{i\mathrm{}\mathrm{}\mathbf{k}\mathrm{\cdot }\mathrm{}\mathbf{r}}=\frac{e^{i\mathrm{}\mathrm{}\mathbf{k}\mathrm{\cdot }\mathrm{}\mathbf{r}}}{L^{3/2}}.$$

The factor $𝒩=L^{-3/2}$ guarantees the normalization of the plane wave, and the dot product $\mathbf{k}\cdot \mathbf{r}$ defines the scalar product between the momentum $\mathbf{k}=\mathbf{p}/\mathrm{\hslash }=k_x\widehat{\mathbf{x}}+k_y\widehat{\mathbf{y}}+k_z\widehat{\mathbf{z}}$ and the position. At the end of our calculations, we must take $L\to \mathrm{\infty }$ to guarantee the continuum character of the state.

Although $V$ is time-independent, we may tackle our problem with a time-dependent formalism by considering the potential as an “adiabatic switch”. The potential is “turned on” while $r<R$, inside the scattering region, and rapidly “turned off” when $r>R$. Furthermore, the scattering may be assumed to be elastic. The transition is from $|i\mathrm{\rangle }$ to a group of continuum states $|f\mathrm{\rangle }$ with energies $E_f$, where both $|i\mathrm{\rangle }$ and $|f\mathrm{\rangle }$ are eigenstates of $H_0$.

The procedure outlined in the last paragraph is carried out in detail in ^A. Here, we reproduce the result for the asymptotic behavior of the wave function,

(4) $$\psi (\mathbf{r})\stackrel{\text{large}r}{\to }\frac{1}{L^{3/2}}\left[e^{i\mathbf{k}\cdot \mathbf{r}}+\frac{e^{ikr}}{r}f({\mathbf{k}}^{\prime },\mathbf{k})\right].$$

This is a quantitative description of what we saw in Fig.1: at large distances, the wave function combines the incident plane wave and a scattered spherical wave. Notice that the factor $f({\mathbf{k}}^{\prime },\mathbf{k})$ is multiplying the spherical wave in Eq. (4). This term is called scattering amplitude, and it indicates how much of the incident wave was scattered. It depends on the initial $\mathbf{k}$ and final${\mathbf{k}}^{\prime }$ momenta,

(5) $$f(k^{\prime },k)=-\frac{m{L}^3}{2\pi {\hslash }^2}{{\displaystyle \int }}^{\text{}}{d}^3{r}^{\mathrm{\prime }}\mathrm{\langle }{k}^{\mathrm{\prime }}\mathrm{|}{r}^{\mathrm{\prime }}\mathrm{\rangle }\mathrm{}V(r^{\mathrm{\prime }})\mathrm{}\mathrm{\langle }{r}^{\prime }|\psi \mathrm{\rangle }.$$

Equation (4) contains a large amount of information. Without solving the Schrödinger equation explicitly, we arrived at a result for the wave function at large distances with only a few assumptions about the scattering potential.

2.1. Partial waves expansion

To obtain Eq. (4), we assumed the scattering potential $V(r)$ to be real, finite-ranged, and local. Now, we consider one last restriction: that it is spherically symmetric, i.e., $V$ is invariant under rotations. Hence, we can write $V(r)=V(r)$. We must then solve Schrödinger’s equation for a potential $V(r)$ such that $V(0<r<R)\ne 0$ and $V(r>R)=0$. In the scattering region ($0<r<R$), we write

(6) $$-\frac{{\mathrm{\hslash }}^2}{2m}{\nabla }^2\psi +V(r)\psi =E\psi .$$

The total energy is $E={\mathrm{\hslash }}^2\mathrm{}{k}^2\mathrm{/}2m$. A straightforward way of seeing this is to notice that this is the energy of the particle in the region $r>R$, where $V(r)=0$. Since the scattering process is elastic, the total energy is conserved.

The wave function $\psi ={\psi }_k\left(r,\theta \right)$ depends only on the momentum $\mathbf{k}$, the position $r$ and the scattering angle $\theta$, due to the spherical symmetry. We want the solution to obey the asymptotic behavior of Eq. (4), which we now write as

(7) $${\psi }_{\mathbf{k}}(r,\theta )\stackrel{\text{large}r}{\to }𝒩\left[e^{ikz}+\frac{e^{ikr}}{r}f(\theta )\right],$$

where $𝒩$ is a normalization constant. We also considered $\mathbf{k}\cdot \mathbf{r}=kz$, corresponding to an incident plane wave in the $z$ direction.

Due to the spherical symmetry of $V(r)$, it is convenient to employ spherical coordinates. We write the Laplace operator in spherical coordinates and $\psi =\psi (r,\theta ,\varphi )$, for which we derive the eigenvalue equation

(8) $$\left(-\frac{{\mathrm{\hslash }}^2}{2m}\frac{1}{r}\frac{{\partial }^2}{\partial r^2}r+\frac{L^2}{2m{r}^2}+V(r)\right)\mathrm{}\psi (r,\theta ,\varphi )=E\psi (r,\theta ,\varphi ).$$

The angular momentum operator $L$ is such that

(9) $$L^2=-{\mathrm{\hslash }}^2\left(\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\sin \theta \frac{\partial }{\partial \theta }+\frac{1}{{\sin }^2\theta }\frac{{\partial }^2}{\partial {\phi }^2}\right).$$

Its $z$-component is given by

(10) $$L_z=-i\mathrm{\hslash }\frac{\partial }{\partial \varphi }.$$

We then construct a complete set of eigenfunctions related to the commuting observables $H,L^2$, and $L_z$ with eigenvalues $E$, $l(l+1){\mathrm{\hslash }}^2$, and $m\mathrm{\hslash }$, respectively,

(11) $$H\psi (r,\theta ,\varphi )=E\psi (r,\theta ,\varphi ),$$ $$\mathrm{}{L}^2\psi (r,\theta ,\varphi )=l(l+1){\mathrm{\hslash }}^2\psi (r,\theta ,\varphi ),$$ $$\mathrm{}{L}_z\psi (r,\theta ,\varphi )=m\mathrm{\hslash }\psi (r,\theta ,\varphi ).$$

We propose a separable solution of the form

(12) $$\psi (r,\theta ,\varphi )=A_l(r)Y_l^m(\theta ,\varphi ),$$

where the $A_l(r)$ are radial functions and the $Y_l^m(\theta ,\varphi )$ are the spherical harmonics. The angular part of the wave function, which depends on $\theta$ and $\varphi$, for any spherically symmetric potential corresponds to the spherical harmonics. Typically, students first encounter this property while studying the particular case of the hydrogen atom. This is a consequence of the potential commuting with $L^2$ and $L_z$, $[V,L^2]=0$ and $[V,L_z]=0$, which is true for any potential that depends only on the radial coordinate$r$.

We perform a change of variables $A_l(r)=u_l(r)/r$, which is convenient given that there is a first derivative in the radial equation for $A_l(r)$, and this change removes it from the corresponding equation for $u_l(r)$. Substituting Eq. (12) into (8) yields an equation that depends only on the radial coordinate,

(13) $$\left(\frac{d^2}{d{r}^2}+k^2-U(r)-\frac{l(l+1)}{r^2}\right)u_l(r)=0,$$

where $U(r)=2mV(r)/{\mathrm{\hslash }}^2$. At the origin, we require that $A_l(r)$ is finite. Since $A_l(r)=u_l(r)/r$, we need $u_l(0)=0$.

We have a free particle in the region $r>R$ since $U(r)=0$. In this case, the solutions of Eq. (13) can be written in terms of the spherical Bessel functions of the first and second kind, usually denoted by $j_l(x)$ and $n_l(x)$, respectively. The solution is of the form

(14) $$u_l(r)=c_l^{\prime }r{j}_l(kr)+c_l^{\prime \prime }r{n}_l(kr),$$

where $c_l^{\prime }$ and $c_l^{\prime \prime }$ are constants. It is convenient to write the solution in terms of a linear combination of spherical Bessel,

(15) $$h_l^{(1)}(x)=j_l(x)+i{n}_l(x),$$ $$h_l^{(2)}(x)=j_l(x)-i{n}_l(x),$$

where $h_l^{(1)}$ and $h_l^{(2)}$ are called spherical Hankel functions of the first and second kind, respectively. Thus, Eq. (14) can be written as

(16) $$u_l(r)=c_l^{(1)}r{h}_l^{(1)}(kr)+c_l^{(2)}r{h}_l^{(2)}(kr),$$

where $c_l^{(1)}$ and $c_l^{(2)}$ are constants.

We know that the solution for a free particle, a plane wave $e^{ikz}=e^{ikr\cos \theta }$, contains components with all possible values of the angular momentum, $l=0,1,2,\mathrm{\dots }$, and we wish to decompose it into a sum of different angular momentum components. This can be done with Rayleigh’s formula [²¹],

(17) $$e^{ikr\mathrm{}\cos \theta }=\sum _{l=0}^{\mathrm{\infty }}{i}^l(2l+1)j_l(kr)P_l(\cos \theta ),$$

where the $P_l$ are Legendre polynomials. Notice that the plane wave, and consequently its expansion, does not depend on $\varphi$. This is due to the spherical symmetry of the problem. If we think in terms of spherical harmonics, their $\varphi$ dependence is $Y_l^m\propto e^{im\varphi }$. The absence of $\varphi$ indicates that we must set $m=0$, which leads directly to the Legendre polynomials since $Y_l^0(\theta ,\varphi )\propto P_l(\cos \theta )$.

Let us analyze the asymptotic behavior of Eq.(17), when $r\to \mathrm{\infty }$. Equation15 allows us to write

(18) $$j_l(x)=\frac{h_l^{(1)}(x)+h_l^{(2)}(x)}{2}.$$

The spherical Hankel functions obey the following limits [²¹]

(19) $$h_l^{(1)}(x)\stackrel{\text{large}x}{\to }{(-i)}^{l+1}{\displaystyle \frac{e^{ix}}{x}},$$ $$h_l^{(2)}(x)\stackrel{\text{large}x}{\to }{i}^{l+1}{\displaystyle \frac{e^{-ix}}{x}}.$$

Gathering all this information, we have that

(20) $$e^{ikr\mathrm{}\cos \theta }\stackrel{\text{large}r}{\to }\sum _{l=0}^{\mathrm{\infty }}\frac{(2l+1)}{2ikr}\left[e^{ikr}-{(-1)}^l{e}^{-ikr}\right]P_l(\cos \theta ).$$

The first term inside the square brackets represents an outgoing spherical wave, while the second is related to an incoming spherical wave.

In the region $r>R$ we have not only the incident plane wave, but also the spherical wave produced by the scattering center. Motivated by Eq. (17) and its asymtoptic behavior, Eq. (20), we write the solution for every $r>R$ as

(21) $$\psi (r,\theta )=𝒩\sum _{l=0}^{\mathrm{\infty }}{i}^l(2l+1)\frac{u_l(r)}{r}{P}_l(\cos \theta ),$$

where $𝒩$ is a normalization constant. Let us analyze the behavior of $\psi (r,\theta )$ when $r\to \mathrm{\infty }$. Equations (16) and (19) yield

(22) $$\psi (r,\theta )\stackrel{\text{large}r}{\to }$$ $$𝒩{\displaystyle \sum _{l=0}^{\mathrm{\infty }}}{\displaystyle \frac{(2l+1)}{ikr}}\left[c_l^{(1)}{e}^{ikr}-{(-1)}^l{c}_l^{(2)}{e}^{-ikr}\right]P_l(\cos \theta ).$$

Let us pause to consider the implications of this result. Equation (20) describes the asymptotic behavior of the wave function for a plane wave without being scattered, while Eq. (22) does the same, but in a situation where scattering could have taken place. If we take $𝒩=1$ and $c_l^{(1)}=c_l^{(2)}=1/2$, then both equations are the same. This is not surprising since this particular choice makes the radial function the same as the one for a free particle, $u_l(r)/r=j_l(kr)$ [see Eqs.(16) and (18)]. However, this situation arises only because the coefficients in front of the outgoing spherical wave and the incoming one, $c_l^{(1)}$ and $c_l^{(2)}$ respectively, are the same. If we take $c_l^{(1)}\ne c_l^{(2)}$, then scattering certainly took place. Moreover, the proportion of outgoing to incoming spherical waves, given by the ratio of these two coefficients, can be used to quantify the impact of the scattering potential on the free particle.

This motivates us to introduce a new quantity related to the ratio between the constants in front of $e^{ikr}/r$ and $e^{-ikr}/r$,

(23) $$\frac{c_l^{(1)}}{c_l^{(2)}}=S_l(k)=e^{2i{\delta }_l(k)},$$

where the $S_l(k)$ are called scattering matrix components and the ${\delta }_l(k)$ are called phase shifts. Notice that ${\delta }_l(k)$ depends explicitly on the momentum or, equivalently, on the energy. The free particle case, $c_l^{(1)}=c_l^{(2)}=1/2$, corresponds to a phase shift ${\delta }_l(k)=0$. Writing the ratio of two constants as a complex exponential may seem strange, but the reason will be apparent shortly.

Let us express the asymptotic wave function, Eq. (22), in terms of the phase shifts, Eq. (23),

(24) $$\psi (r,\theta )\stackrel{\text{large}r}{\to }𝒩{\displaystyle \sum _{l=0}^{\mathrm{\infty }}}{\displaystyle \frac{(2l+1)}{ikr}}{c}_l^{(2)}\left[e^{2i{\delta }_l}{e}^{ikr}-{(-1)}^l{e}^{-ikr}\right]P_l(\cos \theta ).$$

Recall that our analysis started with the idea that a plane wave is scattered and produces a spherical wave, the latter multiplied by the scattering amplitude $f(\theta )$, as given by Eq. (7). Now, we can connect the scattering amplitude and the phase shifts. For that, we recast Eq. (7) using the asymptotic limit of Rayleigh’s formula, Eq.20,

(25) $$\psi (r,\theta )\stackrel{\text{large}r}{\to }𝒩\{[{\displaystyle \sum _{l=0}^{\mathrm{\infty }}}{\displaystyle \frac{(2l+1)}{2ikr}}(e^{ikr}-{(-1)}^l{e}^{-ikr})\times P_l(\cos \theta )]+f(\theta ){\displaystyle \frac{e^{ikr}}{r}}\}.$$

Comparison between Eqs. (24) and (25) allows us to write the scattering amplitude as a function of the phase shifts,

(26) $$f(\theta )=\sum _{l=0}^{\mathrm{\infty }}(2l+1)\frac{(e^{2i{\delta }_l}-1)}{2ik}{P}_l(\cos \theta ).$$

The factor $(e^{2i{\delta }_l}-1)/2ik$ is referred to as the partial wave amplitude $f_l(k)$, which may be rewritten as

(27) $$f_l(k)=\frac{e^{2i{\delta }_l}-1}{2ik}=\frac{e^{i{\delta }_l}\sin {\delta }_l}{k},$$

or

(28) $$f_l(k)=\frac{e^{i{\delta }_l}\sin {\delta }_l}{k}=\frac{1}{k\cot {\delta }_l-ik},$$

whichever is more convenient. Using Eq. (23), we may also express the partial wave amplitude in terms of $S_l(k)$as

(29) $$S_l(k)=1+2ik{f}_l(k).$$

Now we can attribute physical meaning to the partial wave scattering amplitude and the phase shifts. Conservation of the probability during scattering tells us that, at large distances, the only thing that can change is the phase of the wave function (with respect to the incident wave). The difference between the phases is the phase shift ${\delta }_l(k)$. When there is no scattering, $V=0$, ${\delta }_l(k)=0$ and $f_l(k)=0$, meaning that we have the free particle solution in all space, indicated as a dot-dashed line in Figure 2. For a potential $V\ne 0$, the radial solution for $r<R$ will depend on the details of the potential. However, we have a free particle solution outside the range $R$ of the potential, $V(r>R)=0$. Hence, what happens inside the range of the potential determines the phase shift observed outside of it. An attractive potential “pulls” the particle toward the origin, Fig.2a, while a repulsive potential “pushes” it away, Fig.2b. The mathematical advantage of this formulation is that we describe the whole process in terms of a real quantity ${\delta }_l(k)$, and we reduce our problem to calculating it.

Figure 2
Scattered radial solution $u_0(r)$ and the free ($V=0)$ radial solution $g_0(r)$ to highlight the phase shift. At low energies, an attractive potential ($V(r)<0$) pulls the wave function, resulting in ${\delta }_0(k)>0$. The repulsive potential ($V(r)>0$) pushes $u_0(r)$ away, resulting in ${\delta }_0(k)<0$.

To compute the phase shifts, we can use two properties of wave functions: they are continuous, and so are their first derivatives. Hence, we need to match the inside ($r<R$) and outside ($r>R$) solutions and their derivatives at $r=R$. We write the radial solution in the region $r>R$ as

(30) $$A_l(r)=\frac{u_l(r)}{r}=\frac{1}{2}{e}^{2i{\delta }_l}{h}_l^{(1)}(kr)+\frac{1}{2}{h}_l^{(2)}(kr)$$

in terms of spherical Hankel functions, and

(31) $$A_l(r)=e^{i{\delta }_l}(\cos {\delta }_l{j}_l(kr)-\sin {\delta }_l{n}_l(kr))$$

in terms of spherical Bessel functions. To avoid calculating the normalization of the wave function, we can work with the ratio $u_l^{\prime }(r)/u_l(r)$ so that it cancels out. This quantity is known as a logarithmic derivative because

(32) $$\frac{d}{dx}\ln f(x)=\frac{f^{\prime }(x)}{f(x)}.$$

We define a dimensionless logarithmic derivative by including a factor of $r$. At $r=R^{-}$ (at the range $R$ of the potential coming from inside), it is a constant which depends on the expression of $V(r)$,

(33) $${\beta }_l={\left[r\frac{u_l^{\prime }(r)}{u_l(r)}\right]}_{r=R^{-}}.$$

The phase shift can be determined by equating the logarithmic derivative at $r=R^{-}$ with the outside solution (at $r=R^{+}$) given by Eq. (31),

(34) $${\beta }_l={\left[r{\displaystyle \frac{u_l^{\prime }(r)}{u_l(r)}}\right]}_{r=R^{+}}=1+kR\left[{\displaystyle \frac{\cos {\delta }_l{j}_l^{\prime }(kR)-\sin {\delta }_l{n}_l^{\prime }(kR)}{\cos {\delta }_l{j}_l(kR)-\sin {\delta }_l{n}_l(kR)}}\right],$$

where the notation $j_l^{\prime }(kR)$ means the derivative of $j_l$ with respect to $kr$ evaluated at $kR$, and similarly for $n_l^{\prime }(kR)$. After some algebra, we arrive at

(35) $$\cot {\delta }_l(k)=\frac{kR{n}_l^{\prime }(kR)-({\beta }_l-1)n_l(kR)}{kR{j}_l^{\prime }(kR)-({\beta }_l-1)j_l(kR)}.$$

This is an analytic expression to calculate the $l$-th partial wave phase-shift ${\delta }_l(k)$ provided we know the inside solution to compute the constant ${\beta }_l$. This result is a cornerstone for the numerical procedure we will treat later in this article.

2.2. The low-energy limit and the scatteringlength

From Eq. (8), it is possible to see that the particle is subjected to an effective potential for the $l$-th partial wave of

(36) $$V_{\mathrm{eff}}(r)=V(r)+{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\displaystyle \frac{l(l+1)}{r^2}},$$

where the second term on the right-hand side is a repulsive centrifugal barrier. It is absent for $l=0$, becoming more repulsive as we increase the value $l$ of the considered partial wave. This work focuses on low-energy scattering, where the reduced wavelength $ƛ=\lambda /2\pi =1/k$ is much larger than the potential range, that is $ƛ\gg R$ or $kR\ll 1$. If the energy is close to zero, $E\approx 0$, then the particle cannot overcome the centrifugal barrier. In this case, the partial waves with $l>0$ are unimportant, and the $l=0$ component is dominant in understanding low-energy scattering.

In this low-energy scenario, we consider partial waves with $l\ne 0$ to vanish, and the resulting $l=0$ term is referred to as “$s$-wave”. Thus, the $s$-wave scattering radial component is given by Eq. (31),

(37) $$A_0(r)={\displaystyle \frac{u_0(r)}{r}}=e^{i{\delta }_0}(\cos {\delta }_0{j}_0(kr)-\sin {\delta }_0{n}_0(kr))=e^{i{\delta }_0}\left[{\displaystyle \frac{1}{kr}}\sin (kr+{\delta }_0)\right],$$

where we used

(38) $$j_0(x)={\displaystyle \frac{\sin x}{x}},$$ $$n_0(x)=-{\displaystyle \frac{\cos x}{x}}.$$

Schrödinger’s equation for the radial solution becomes very simple in this situation. Outside the range of the potential, $V(r>R)=0$. Moreover, there is no centrifugal barrier since $l=0$, and for low-energy scattering, $k\approx 0$. Hence, Eq. (13) reduces to

(39) $$u_0^{\prime \prime }(r)=0.$$

The solution can be written as

(40) $$u_0(r)=c(r-a),$$

where $c$ and $a$ are constants. Its logarithmic derivative, given by Eq. (33), is

(41) $$r\frac{u_0^{\prime }(r)}{u_0(r)}=\frac{r}{r-a}.$$

This needs to be equal to the logarithmic derivative of Eq. (37), $kr\cot (kr+{\delta }_0)$. Hence,

(42) $$kr\mathrm{}\cot (kr+{\delta }_0)=\frac{r}{r-a}.$$

In the limit $k\to 0$, and also setting $r=0$, we have

(43) $$\underset{k\to 0}{lim}k\mathrm{}\cot \mathrm{}{\delta }_0(k)=-\frac{1}{a},$$

where $a$ is called scattering length [³]. Let us consider the implications of this result. In the previous section, we reduced the scattering problem to calculating the phase shifts ${\delta }_l(k)$. Here we are considering only low-energy phenomena, to which we argued that the $l=0$ component is the dominant one. However, Eq. (43) reduces the problem even further: in the zero-energy limit, a single number, the scattering length, encodes all the information we need about scattering.

The scattering length has, as its name suggests, dimension of length. However, it may differ by orders of magnitude from the range $R$ of the potential, as we will see in several examples. A geometrical interpretation is possible if we choose $c=-1/a$ in Eq. (40) (remember that the normalization always cancels out when we work with ratios such as $u_0^{\prime }/u_0$),

(44) $$u_0(r>R)=1-\frac{r}{a}.$$

We see that the scattering length is simply the intercept of the outside wave function, or its extrapolation, as illustrated in Fig.3. Finally, the scattering amplitude for $l=0$, Eq. (28), in the low-energy limit is

Figure 3
Geometrical interpretation of the scattering length $a$ based on three possibilities.

(45) $$f_0(k)=\underset{k\to 0}{lim}\frac{1}{k\mathrm{}\cot {\delta }_0-ik}=-a.$$

2.3. The effective range

Equation (43) and the concept of scattering length are remarkable. However, the underlying $kR\ll 1$ assumption makes them good approximations to physical situations only when the energy or the range of the potential goes to zero. We might wonder if we can modify Eq. (43) to consider a finite but small value of $kR$. In other words, we would like to express $k\cot {\delta }_0(k)$ as a series in powers of $k$. We already know that the first term is $-1/a$, so the task becomes computing the next.

We follow the procedure outlined in Ref. [⁶]. First, let us choose a different normalization for $u_0(r)$ in Eq. (37),

(46) $$u_0(r)=\cot {\delta }_0(k)\mathrm{}\sin (kr)+\cos (kr),$$

and the reason for this choice will be apparent shortly. Let us take the $l=0$ radial equation, Eq. (13), for two different wave functions $u_{k_1}(r)$ and $u_{k_2}(r)$, labeled by their wave vectors $k_1=\sqrt{2m{E}_1}/\mathrm{\hslash }$ and $k_2=\sqrt{2m{E}_2}/\mathrm{\hslash }$,

(47) $$u_{k_2}^{\prime \prime }(r)-U(r)u_{k_2}(r)+k_2^2{u}_{k_2}(r)=0.$$ $$u_{k_2}^{\prime \prime }(r)-U(r)u_{k_2}(r)+k_2^2{u}_{k_2}(r)=0.$$

Next, we multiply the first equation by $u_{k_2}$ and the second by $u_{k_1}$ and take their difference,

(48) $$u_{k_1}^{\prime \prime }(r)u_{k_2}(r)-u_{k_1}(r)u_{k_2}^{\prime \prime }(r)=(k_2^2-k_1^2)u_{k_1}(r)u_{k_2}(r).$$

We may write the left-hand side as

(49) $$u_{k_1}^{\prime \prime }(r)u_{k_2}(r)-u_{k_1}(r)u_{k_2}^{\prime \prime }(r)$$ $$={\displaystyle \frac{d}{dr}}\left[u_{k_1}^{\prime }(r)u_{k_2}(r)-u_{k_2}^{\prime }(r)u_{k_1}(r)\right].$$

Now we integrate Eq. (48) from $0$ to $R$,

(50) $${\left[u_{k_2}^{\prime }(r)u_{k_1}(r)-u_{k_1}^{\prime }(r)u_{k_2}(r)\right]}_0^R$$ $$=(k_2^2-k_1^2){\displaystyle {\int }_0^R}\mathit{d}r{u}_{k_1}(r)u_{k_2}(r).$$

The integral converges since $A_0(r)=u_0(r)/r$ is finite at the origin ($u_0(0)=0$ independently of the energy).

Next, we repeat the same procedure for the free-particle ($U=0$) radial equation with solutions denoted by $g_{k_1}(r)$ and $g_{k_2}(r)$. The result is the same as Eq. (50) if we replace $u$ by $g$. Finally, we take the difference between this result and Eq. (50),

(51) $${\left[g_{k_2}^{\prime }(r)g_{k_1}(r)-g_{k_1}^{\prime }(r)g_{k_2}(r)\right]}_0^R$$ $$-{\left[u_{k_2}^{\prime }(r)u_{k_1}(r)-u_{k_1}^{\prime }(r)u_{k_2}(r)\right]}_0^R$$ $$=(k_2^2-k_1^2){\displaystyle {\int }_0^R}\mathit{d}r[g_{k_1}(r)g_{k_2}(r)-u_{k_1}(r)u_{k_2}(r)].$$

The functions $g_{k_1}$ and $g_{k_2}$ are given by Eq. (46), and so are $u_{k_1}$ and $u_{k_2}$ for $r⩾R$. Then, Eq. (51) becomes

(52) $$k_2\cot {\delta }_0(k_2)-k_1\cot {\delta }_0(k_1)$$ $$=(k_2^2-k_1^2){\displaystyle {\int }_0^R}\mathit{d}r[g_{k_1}(r)g_{k_2}(r)-u_{k_1}(r)u_{k_2}(r)].$$

If we take the limit $k_1\to 0$, we can write $k_1\cot {\delta }_0(k_1)$ in terms of the scattering length, Eq. (43),

(53) $$k\cot {\delta }_0(k)=-\frac{1}{a}+k^2{\int }_0^R\mathit{d}r[g_0(r)g_k(r)-u_0(r)u_k(r)],$$

where we dropped the subscript in $k_2$ for simplicity. We define the quantity

(54) $$\rho (k)\equiv 2{\displaystyle {\int }_0^R}\mathit{d}r[g_0(r)g_k(r)-u_0(r)u_k(r)].$$

Finally, we have an expression for $k\cot {\delta }_0(k)$ at low-energies,

(55) $$k\mathrm{}\cot {\delta }_0(k)=-\frac{1}{a}+\frac{1}{2}{r}_0{k}^2+𝒪(k^4),$$

where the term $r_0$ is referred to as “effective range” and is defined as

(56) $$r_0\equiv \underset{k\to 0}{lim}\rho (k)=2{\int }_0^R\mathit{d}r[g_0^2(r)-u_0^2(r)],$$

where $g_0(r)$ is calculated by taking $k\to 0$ in Eq. (46). The result is exactly Eq. (44), justifying the normalization choice.

Comparing Eqs. (43) and (55), we see that we were able to include the next term in the expansion, which is proportional to $k^2$. Only even powers of $k$ appear on the right-hand side of Eq. (55) because of the symmetry with respect to changing the sign of $k$: the phase shift would also change sign and, since $\cot (-x)=-\cot (x)$, the left-hand side is unchanged. Hence, the right-hand side can only have even functions of $k$. To compute the effective range, we only need two zero-energy solutions of the radial equation: the free particle, $g_0$, and the solution with the potential $V(r)$, $u_0$.

Equation (55) asserts that the $s$-wave scattering phase shift does not depend on the shape of the potential $V(r)$ at low energies. Instead, two potentials with different forms will produce the same phase shift if their scattering lengths and effective ranges are the same. For this reason, Eq. (55) is commonly referred to as shape-independent approximation . Nevertheless, it is important to note the error order of $k^4$. At higher energies, outside the scope of this work, higher order terms become relevant, and Eq. (55) may not be appropriated.

2.4. Bound states

When we defined the phase shifts in terms of the ratio of the coefficients appearing in front of the outgoing and income spherical waves, Eq. (23), we also introduced the amplitude $S_l(k)$. The analytic properties of $S_l(k)$ contain information about bound states, as we will see in the following.

Let us rewrite Eq. (24) as

(57) $$\psi (r,\theta )\stackrel{\text{large}r}{\to }{\displaystyle \frac{1}{{(2\pi )}^{3/2}}}{\displaystyle \sum _{l=0}^{\mathrm{\infty }}}{\displaystyle \frac{(2l+1)}{2ik}}{P}_l(\cos \theta )$$ $$\times \left[S_l(k){\displaystyle \frac{e^{ikr}}{r}}-{\displaystyle \frac{e^{-i(kr-l\pi )}}{r}}\right].$$

For $l=0$ and large distances, the radial wave function is proportional to

(58) $$S_0(k)\frac{e^{ikr}}{r}-\frac{e^{-ikr}}{r}.$$

For an arbitrary finite-ranged potential $V$, the radial solution at $r>R$ for a bound state ($E<0$) obeys

(59) $$u^{\prime \prime }(r)=-\frac{2mE}{{\mathrm{\hslash }}^2}u(r)={\kappa }^{2}u(r),\kappa \equiv \frac{\sqrt{-2mE}}{\mathrm{\hslash }}.$$

The solution can be written as

(60) $$u(r>R)=A{e}^{\kappa r}+B{e}^{-\kappa r},$$

where $A$ and $B$ are constants. However, $u(r\to \mathrm{\infty })$ must be finite, so we set $A=0$. We conclude that the radial function for a bound state at large distances is

(61) $$A(r)=\frac{u(r)}{r}\propto \frac{e^{-\kappa r}}{r}\text{(large}r\text{)}.$$

The existence of a bound state is related to the Schrödinger equation admitting a solution with a discrete $E<0$ energy. By looking at Eq. (58), we may be tempted to substitute $k\to i\kappa$, with $k$ purely imaginary, to connect it to the bound state behavior of Eq. (61),

(62) $$\frac{e^{ikr}}{r}=\frac{e^{i(i\kappa )r}}{r}=\frac{e^{-\kappa r}}{r}.$$

However, the important difference is that the bound state is present even without the incident wave we have used for our scattering problem formulation. In Eq. (58), $S_0(k)$ controls the ratio of the outgoing to the incoming wave. Hence, in the bound state case, we have only the outgoing component, meaning an infinite ratio: $S_0(k)\to \mathrm{\infty }$. Since we are treating $S_0(k)$ as a function of a complex variable, this means that it has a pole at $k=i\kappa$. Hence, the conclusion is that a bound state appears as a pole in $S_0(k)$. It is convenient to represent the variable $k$ in the complex $k$-plane, i.e., where the abscissa is the real part of $k$, and the ordinate is the imaginary part of $k$. The pole at $k=i\kappa$ is a point on the $\mathrm{Im}(k)$ axis, while the scattering continuum is on the positive $\mathrm{Re}(k)$ axis, since $k>0$, as illustrated in Figure 4.

Figure 4
Complex $k$-plane. A bound state is represented by a pole in the $S$-matrix at $k=i\kappa$, while the scattering continuum corresponds to a real value of $k>0$.

Let us derive an expression for $S_0(k)$, so we can identify the $k=i\kappa$ pole. Eqs. (28) and (43) yield

(63) $$f_0(k)=\frac{1}{-1/a-ik}.$$

The scattering amplitude is related to $S_0(k)$ through Eq. (29),

(64) $$S_0(k)=1+2ik{f}_0(k)=\frac{-k-i/a}{k-i/a}.$$

This expression has a pole at $k=i\kappa$ if we identify

(65) $$\kappa =\frac{1}{a}.$$

Hence, in the zero-energy limit, the energy of a bound state and the scattering length are connected simply by

(66) $$E=\frac{{\mathrm{\hslash }}^2\mathrm{}{\kappa }^2}{2m}=\frac{{\mathrm{\hslash }}^2}{2m{a}^2}.$$

A single parameter originated from the potential determines the bound-state energy.

2.5. Two-body scattering

So far, we considered only the problem of a single particle being scattered by a finite-ranged potential $V(r)$ located at $r=0$. With a few modifications, we can use the results we obtained to describe two particles interacting through a pairwise potential which depends only on their spatial separation $r$.

Besides being separable in radial and angular coordinates, the Hamiltonian of a two-body system with a spherically symmetric potential is also separable in the center of mass (CM) and relative coordinates. The two-body Hamiltonian is given by

(67) $$H=-\frac{{\mathrm{\hslash }}^2}{2m_1}{\nabla }_{r_1}^2-\frac{{\mathrm{\hslash }}^2}{2m_2}{\nabla }_{r_2}^2+V(r_1-r_2),$$

where the indices $1$ and $2$ refer to different particles. We define the coordinates

(68) $$\mathbf{R}={\displaystyle \frac{m_1{\mathbf{r}}_1+m_2{\mathbf{r}}_2}{M}},$$ $$\mathbf{r}={\mathbf{r}}_1-{\mathbf{r}}_2,$$

where $M=m_1+m_2$. In this coordinate system, and using the fact that the potential depends only on $r=|\mathbf{r}|$,

(69) $$H=H_{\text{CM}}+H_r,$$ $$H_{\mathrm{CM}}=-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2M}}{\nabla }_R^2,$$ $$H_r=-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m_r}}{\nabla }_r^2+V(r),$$

where $m_r=m_1{m}_2/(m_1+m_2)$ is the reduced mass. The CM motion satisfies the free-particle equation, so its behavior is trivial and only adds a constant to the total energy. However, the relative movement is non-trivial and contains all the physics due to the scattering of the two particles. Notice that the relative motion Hamiltonian is exactly the one we used to discuss the one particle problem if we substitute $m\to m_r$. Thus we can apply our previous results to two-body scattering.

2.6. Some applications

Before presenting the numerical procedure to calculate the scattering length and effective range of two-body potentials, it is instructive to show some examples of analytical calculations.

2.6.1. Spherically symmetric finite well

For example, let us consider a two-body system with an interaction described by an attractive spherically symmetric finite well. It has a constant value inside the range $R$ and zero outside. We can write $V(r)$ as

(70) $$V(r)=\{\begin{array}{cc}-v_0\frac{{\hslash }^2}{m_r{R}^2},\hfill & \text{for}r<R,\hfill \\ 0,\hfill & \text{for}r>R,\end{array}$$

where $v_0>0$ is a dimensionless parameter related to the depth of the well. We generally want to tune the parameters $v_0$ and $R$ to achieve the physical conditions we are interested in studying. For a relatively shallow or short-ranged potential, we may only observe continuum scattering states ($E>0$), but increasing its depth or range may make it strong enough to produce a bound state ($E<0$).

Let us start with the $E>0$ case. For simplicity, let us consider the equal mass case, $m_1=m_2=m$. The $s$-wave ($l=0$) equation we need to solve is

(71) $$\left(\frac{d^2}{d{r}^2}-\frac{2m_r}{{\mathrm{\hslash }}^2}V(r)+\frac{2m_r}{{\mathrm{\hslash }}^2}E\right)u(r)=0,$$

where $u(r)=rA(r)$. Writing the explicit forms inside and outside the range of the potential yields

(72) $$\begin{array}{cc}\mathrm{}{u}^{\prime \prime }(r)+(k_0^2+k^2)u(r)=0& \text{for}r<R,\\ u^{\prime \prime }(r)+k^{2}u(r)=0& \text{for}r>R,\end{array}$$

where $k^2\equiv 2m_{r}E/{\mathrm{\hslash }}^2$ and $k_0^2\equiv 2v_0/R^2$. In the region $r<R$, the solution may be written as

(73) $$u(r)=A\sin \left(\sqrt{k^2+k_0^2}r\right)+B\cos \left(\sqrt{k^2+k_0^2}r\right).$$

Since $u(0)=0$, because the radial solution $A(r)=u(r)/r$ is regular at the origin, we set $B=0$. In the region $r>R$, $u(r)$ is given by Eq. (46). Hence, the solution is of the form

(74) $$u(r)=\{\begin{array}{cc}A\sin \left(\sqrt{k^2+k_0^2}r\right)\hfill & \text{for}r<R,\hfill \\ \cot {\delta }_0(k)\sin (kr)+\cos (kr)\hfill & \text{for}r>R,\hfill \end{array}$$

where ${\delta }_0(k)$ is the phase shift. Now we equate the logarithmic derivative at $r=R^{-}$ and $r=R^{+}$,

(75) $${\left[r\frac{u^{\prime }(r)}{u(r)}\right]}_{r=R^{-}}={\left[r\frac{u^{\prime }(r)}{u(r)}\right]}_{r=R^{+}},$$

which yields,

(76) $${\displaystyle \frac{\sqrt{k^2+k_0^2}\cos \left(\sqrt{k^2+k_0^2}R\right)}{\sin \left(\sqrt{k^2+k_0^2}R\right)}}$$ $$={\displaystyle \frac{k\cot {\delta }_0(k)\cos (kR)-k\sin (kR)}{\cot {\delta }_0(k)\sin (kR)+\cos (kR)}}.$$

Without any approximation, we may solve for the phase shift ${\delta }_0(k)$,

(77) $${\delta }_0(k)=-kR+\arctan \left[\frac{k\tan \left(\sqrt{k^2+k_0^2}R\right)}{\sqrt{k^2+k_0^2}}\right].$$

To calculate the scattering length $a$, we need to take the $k\to 0$ limit as prescribed in Eq. (43). One straightforward way to do this is to rearrange Eq. (76) so that we have factors of $k\cot {\delta }_0(k)$. It is important to keep track of the orders employed in the approximation: Eq. (43) ignores $𝒪(k^2)$ terms. Hence, to be consistent, we need to use the following approximations,

(78) $$\cos (kR)=1+𝒪(k^2),$$ $$\sin (kR)=kR+𝒪(k^3).$$

Performing these approximations, Eq. (76) becomes

(79) $$\sqrt{k_0^2}\cot \left(\sqrt{k_0^2}R\right)=\frac{-1/a}{-R/a+1}.$$

Finally, solving for $a$,

(80) $$a=R-\frac{\tan \left(\sqrt{k_0^2}R\right)}{\sqrt{k_0^2}}=R\left(1-\frac{\tan \left(\sqrt{2v_0}\right)}{\sqrt{2v_0}}\right),$$

where we used $k_0^2=2v_0/R^2$ in the last step. Equation (80) makes it apparent that the scattering length depends only on the parameters of the potential, its depth $v_0$ and range $R$. Further inspection of the equation leads us to conclude that the scattering length diverges for specific values of $v_0$. As we will see in the following, these are related to the appearance of bound states.

We could repeat the above procedure if we want to consider the $E<0$ cases. However, suppose we recall the discussion of bound states in Sec.^2.4. In that case, we conclude that we only need to replace $k=i\kappa$ in the solution inside the range of the potential and that the solution outside is proportional to $e^{-\kappa r}$. Thus, the equivalent of Eq. (74) for $E<0$ is

(81) $$u(r)=\{\begin{array}{cc}{A}^{\prime }\sin \left(\sqrt{k_0^2-{\kappa }^2}r\right)\hfill & \text{for}r<R,\hfill \\ B^{\prime }{e}^{-\kappa r}\hfill & \text{for}r>R,\hfill \end{array}$$

where $A^{\prime }$ and $B^{\prime }$ are constants. They can be determined by matching $u(r)$ at $r=R$ and normalizing the radial solution. Next, we match the logarithmic derivatives at $r=R$,

(82) $$\frac{\sqrt{k_0^2-{\kappa }^2}\cos \left(\sqrt{k_0^2-{\kappa }^2}R\right)}{\sin \left(\sqrt{k_0^2-{\kappa }^2}R\right)}=\frac{-\kappa e^{-\kappa R}}{e^{-\kappa R}}.$$

After some manipulations,

(83) $$\tan \left(\sqrt{k_0^2-{\kappa }^2}R\right)+\frac{\sqrt{k_0^2-{\kappa }^2}}{\kappa }=0.$$

This is a transcendental equation that shows where the bound-state energies are located. It cannot be solved analytically, although obtaining numerical solutions is straightforward. Nevertheless, the term $\sqrt{k_0^2-{\kappa }^2}/\kappa$ is always positive. The only way for Eq. (83) to be valid is if $\tan (\sqrt{k_0^2-{\kappa }^2}R)$ is negative. That is to say,

(84) $$\frac{\pi }{2}+n\pi <\sqrt{k_0^2-{\kappa }^2}R<\pi +n\pi ,$$

where $n$ is an integer. Since $\sqrt{k_0^2-{\kappa }^2}R$ is always positive, $n=0,1,\mathrm{\dots }$. The first bound state is $n=0$. Thus,

(85) $$\frac{\pi }{2R}<\sqrt{k_0^2-{\kappa }^2}<\frac{\pi }{R}.$$

Since $k_0>\sqrt{k_0^2-{\kappa }^2}$ then, substituting $k_0=\sqrt{2v_0}/R$ and using Eq. (85), we have

(86) $$v_0>\frac{{\pi }^2}{8}.$$

This result shows that there are no bound states if $v_0$ is not above a certain threshold value.

Now we have everything we need to connect the appearance of bound states and the scattering length. Equation (80) expresses the scattering length $a$ in terms of $v_0$. We know that $\tan x$ diverges for $x=\pi /2$, which implies that $a$ diverges for $v_0={\pi }^2/8$. This value is exactly the threshold we found for the appearance of a bound state, Eq. (86). The conclusion is that the scattering length diverges when a bound state appears. However, it depends on which side we approach since

(87) $$\underset{x\to {\frac{\pi }{2}}^{-}}{lim}\tan x=+\mathrm{\infty },$$ $$\underset{x\to {\frac{\pi }{2}}^{+}}{lim}\tan x=-\mathrm{\infty }.$$

Let us consider the situation where we fixed the range $R$ of the potential and start with a very shallow well. Then $a<0$, and no bound states are supported. If we increase the depth of the well, $|a|$ also increases to the point we are at the threshold for a bound state to appear. At this point, $|a|\to \mathrm{\infty }$ and a bound state with $E=0$ appears. If we continue to increase $v_0$, then $a$ changes sign, and it will decrease until the formation of a second bound state. This situation is illustrated in Fig.5. We centered our discussion around the appearance of the first bound state, but Eq. (84) implies that this behavior will happen whenever the system is at the threshold of forming an additional bound state. Figure 5 shows this since the points where $|a|\to \mathrm{\infty }$ are exactly where $\sqrt{2v_0}=\pi /2+n\pi$.

Figure 5
Behavior of the scattering length $a$ as function of $\sqrt{2\mathrm{}{v}_0}$ for the spherical well, Eq. (80). For $\sqrt{2\mathrm{}{v}_0}=\pi /2+n\mathrm{}\pi$ ($n=0,1,2,\mathrm{\dots }$), $a$ diverges. This is related to the potential admitting an additional bound state at those points.

The advantage of the low-energy scattering formulation is that we may study the behavior of the bound states of the potentials without explicitly solving Schrödinger’s equation. With the benefit of hindsight, Eq. (66) had already told us that $|a|\to \mathrm{\infty }$ corresponds to the formation of a new bound state since that, at the threshold, the state must have $E=0$.

Finally, to complete our analysis of the spherical well, we must compute the effective range as defined in Eq. (56). Notice that the normalization of the solution outside the range of the potential in Eq. (74) agrees with the one we used while deriving the effective range expansion. To determine the constant $A$ we impose the continuity of $u\mathrm{}(r)$ at $r=R$,

(88) $$A=\frac{\cot {\delta }_0\mathrm{}(k)\mathrm{}\sin (k\mathrm{}R)+\cos (k\mathrm{}R)}{\sin (\sqrt{k^2+k_0^2}\mathrm{}R)}.$$

The normalized solution is written as

(89) $$u(r)=\{\begin{array}{ll}\frac{\cot {\delta }_0(k)\sin (kR)+\cos (kR)}{\sin (\sqrt{k^2+k_0^2}R)}\sin \left(\sqrt{k^2+k_0^2}r\right)\hfill & \text{for}r<R,\hfill \\ \cot {\delta }_0(k)\sin (kr)+\cos (kr)\hfill & \text{for}r>R.\hfill \end{array}$$

The effective range, Eq. (56), is defined in the $k\to 0$ limit of $u\mathrm{}(r)$,

(90) $$u\mathrm{}(r)=\{\begin{array}{cc}\frac{\left(1-R/a\right)}{\sin (k_0\mathrm{}R)}\mathrm{}\sin (k_0\mathrm{}r)\hfill & \text{for}\mathrm{}r<R,\hfill \\ 1-r/a\hfill & \text{for}\mathrm{}r>R.\hfill \end{array}$$

Then, the effective range is given by the integral

(91) $$r_0=2\mathrm{}{\int }_0^{R}dr\mathrm{}\left[{\left(1-\frac{r}{a}\right)}^2-{\left(1-\frac{R}{a}\right)}^2\mathrm{}\frac{{\sin }^2(k_0\mathrm{}r)}{{\sin }^2(k_0\mathrm{}R)}\right].$$

After computing the integral, we use Eq. (80) to replace $a$ in favor of $R$ and $k_0$. The result is

(92) $$r_0=R(1-\frac{1}{3}{\left(\frac{k_{0}R}{\tan (k_{0}R)-k_{0}R}\right)}^2+\frac{1}{k_{0}R\tan (k_{0}R)-{(k_{0}R)}^2}),$$

which shows that $r_0$ is also dependent on parameters of the potential only. We illustrate the behavior of $r_0$ in Fig.6. Another characteristic is that the range of the potential and the effective range are equal ($r_0=R$) only when the scattering length diverges. This stresses that even for the case of a spherical well when there is no doubt that the range of the potential is $R$, a rigorous analysis should be carried out in terms of the effective range if we want to capture the low-energy properties of the system.

Figure 6
Behavior of the effective range $r_0$ as function of $\sqrt{2\mathrm{}{v}_0}$ for the spherical well, Eq. (92). The vertical dashed lines denote where the scattering length diverges, $|a|\to \mathrm{\infty }$, corresponding to an effective range equal to the range of the potential, $r_0/R=1$.

2.6.2. Zero-range and finite-rangeapproximations

Equation (65) has an immediate application when dealing with two-body bound states in the zero-energy limit. The energy can be readily calculated as follows:

(93) $$E_{zr}=-\frac{{\mathrm{\hslash }}^2{\kappa }^2}{2m_r}=-\frac{{\mathrm{\hslash }}^2}{2m_r{a}^2},$$

where we chose the subscript to denote the zero range, i.e., in deriving this result, $r_0$ was taken to be zero. If we want to modify the equation above to include finite-range corrections, we have to combine these results with Eq. (55) and make the substitution $k\to i\kappa$, which yields [²⁴, ²⁵]

(94) $$\frac{1}{a}=\kappa -\frac{1}{2}{r}_0{\kappa }^2.$$

We can solve this quadratic equation for $\kappa$, choose the appropriate root, and use it to compute the bound state energy,

(95) $$E_{fr}=-\frac{{\mathrm{\hslash }}^2{\kappa }^2}{2m_r}=-\frac{{\mathrm{\hslash }}^2}{2m_r{r}_0^2}{\left(1-\sqrt{1-\frac{2r_0}{a}}\right)}^2,$$

where the subscript $fr$ stands for finite-range.

Equations (93) and (95) seem very powerful in the sense that we only need one or two zero-energy parameters to compute the energy of a bound state, besides the reduced mass. However, the question of whether this model produces sensible results in realistic settings remains. At this point, applying the zero-range and finite-range results to physical systems is illustrative. To clarify, we do not compute $a$ and $r_0$ for the examples below; we just take values from the literature. Table1 contains a summary of the parameters we employed and the results obtained.

Table 1.
Summary of the low-energy properties of two physical systems: the ${}^4$He dimer and the deuteron. The values of the scattering length $a$, effective range $r_0$, and bound-state energy $E$ are taken from the indicated references. The zero-range and finite-range approximations, $E_{z\mathrm{}r}$ and $E_{f\mathrm{}r}$, were calculated using Eqs. (93) and (95).

System	$a$	$r_0$	$E$	Ref.	$E_{z\mathrm{}r}$	$E_{f\mathrm{}r}$
${}^4$He dimer	90.4 Å	8.0 Å	$-$1.62 mK	[22]	$-$1.48 mK	$-$1.63 mK
Deuteron	5.4112 fm	1.7436 fm	$-$2.224 MeV	[23]	$-$1.416 MeV	$-$2.223 MeV

We present two examples; the first is the bound state of two ${}^4$He atoms, called helium dimer. There are many realistic pairwise potentials for ${}^4$He. Let us consider the one of Ref. [²²]. It contains adiabatic, relativistic, and quantum electrodynamics interactions. Explaining the details of the potential is beyond the scope of this work, but we only need three values from this reference to illustrate how useful the effective-range approximation is. The authors reported that their potential produces a dimer energy of $E_d=-1.62$ mK, a scattering length of $a=90.4$ Å, and an effective range of $r_0=8.0$ Å.

Substituting $a$, and the mass of an ${}^4$He atom, in Eq. (93) yields $E_{zr}\approx -1.48$ mK, such that $E_{zr}/E_d\approx 0.92$. Let us pause to appreciate this result. We estimated the energy of a bound-state using only one parameter calculated from a zero-energy theory, and it differs only 8% from the full solution of the Schrödinger equation. If we include the finite-range correction, Eq. (95), then the result is $E_{fr}=-1.63$ mK, which is correct within 1%. Both the zero- and finite-range results successfully describe the physical system because $kR\sim 0.1$ in this case.

The second example we chose is from nuclear physics. The deuteron (the nucleus of deuterium) is the only bound state of two nucleons, and it is formed by a proton and a neutron. The details of nuclear interaction are much more complicated than what we present here. For example, the deuteron has a quadrupole moment due to a $l=2$ partial-wave component, which goes beyond $s$-wave ($l=0$) scattering. However, this component is small, and we will give an $s$-wave treatment to the deuteron.

Reference [²³] reports the values $a=5.4112$ fm and $r_0=1.7436$ fm for the proton-neutron scattering length and effective range, respectively. Unlike the dimer case, where the particles are identical, the proton and neutron masses are different, although this difference is small. The reduced mass is then $m_r=m_p{m}_n/(m_p+m_n)$. We employed the values of the proton and neutron masses from the 2018 CODATA recommended values [²⁶].

The experimental energy of the deuteron is $E_d=-2.224$ MeV. Substituting the appropriate values into Eq. (93) yields $E_{zr}=-1.416$ MeV, meaning that the zero-energy theory only accounts for $\approx$ 64% of the deuteron energy. However, if we include the finite-range correction, Eq. (95), then the energy is $E_{fr}=-2.223$ MeV, practically coinciding with the experimental value. Unlike the helium dimer, the range of the potential needed to be taken into account to reproduce the physical system because $kR\sim 0.4$ for the deuteron is larger than in the ${}^4$He case.

We should emphasize that the scales are very different in both examples. The ${}^4$He dimer, in the realm of atomic physics, is in a spatial scale of Å(${10}^{-10}$ m), and the energy is of the order of ${10}^{-7}$ eV. For the deuteron, the lengths are in the femtometer (${10}^{-15}$ m) scale, while the energy is of a few MeV (${10}^6$ eV). This exemplifies how universal are these low-energy scattering results.

3. Numerical procedure

Analytical expressions for low-energy scattering parameters are only available for a few potentials. Even in those cases, the calculations may be cumbersome, as we saw in Sec.^2.6.1 for one of the simplest potentials, the spherical well. In general, we need to calculate $a$ and $r_0$ numerically.

3.1. Numerical solutions of Schrödinger’sequation

3.1.1. Second-order central difference

We need two main ingredients to compute the low-energy scattering parameters, the radial wave function inside and outside the range of the potential. The latter has an analytical expression, but the former can only be computed numerically for most potentials. The radial equation for $u(r)=rA(r)$ resembles the one-dimensional Schrödinger’s equation, so we can employ any of the various numerical methods that are available in the literature [²⁷].

A common strategy in computational physics is to consider the function $u(r)$ in a discrete set of points $r_i=i\mathrm{\Delta }r$, where $i$ is an integer ranging from $i=0,\mathrm{\dots },N$ and $\mathrm{\Delta }r$ is a small quantity. Then, our goal becomes finding the value of $u(r_i)$ for every $i$. First, writing the expressions for the first and second numerical derivatives in this scenario is useful. For that, let us consider two Taylor expansions of $u(r)$ around the points $r\pm \mathrm{\Delta }r$,

(96) $$u(r+\mathrm{\Delta }r)=u(r)+(\mathrm{\Delta }r)u^{\prime }(r)+{\displaystyle \frac{{(\mathrm{\Delta }r)}^2}{2}}{u}^{\prime \prime }(r)$$ $$+{\displaystyle \frac{{(\mathrm{\Delta }r)}^3}{6}}{u}^{\prime \prime \prime }(r)+\mathrm{\cdots },$$ $$u(r-\mathrm{\Delta }r)=u(r)-(\mathrm{\Delta }r)u^{\prime }(r)+{\displaystyle \frac{{(\mathrm{\Delta }r)}^2}{2}}{u}^{\prime \prime }(r)$$ $$-{\displaystyle \frac{{(\mathrm{\Delta }r)}^3}{6}}{u}^{\prime \prime \prime }(r)+\mathrm{\cdots }$$

The difference of the two Taylor expansions yields an expression for the first derivative, while their sum results in the second derivative,

(97) $${{\displaystyle \frac{du}{dr}}|}_{r=r_i}\approx {\displaystyle \frac{u_{i+1}-u_{i-1}}{2\mathrm{\Delta }r}},$$

(98) $${{\displaystyle \frac{d^{2}u}{d{r}^2}}|}_{r=r_i}\approx {\displaystyle \frac{u_{i+1}-2u_i+u_{i-1}}{{(\mathrm{\Delta }r)}^2}},$$

where we introduced a compact notation $u(r_i)\equiv u_i$ and, since $\mathrm{\Delta }r$ is small, we neglected terms of $𝒪[{(\mathrm{\Delta }r)}^3]$. These expressions tell us that if we want to calculate the numerical derivative at some point $r_i$, we need to know the value of the function at the neighboring sites; for the second derivative, we also need the value of the function at the site $r_i$. Since the expression for the first derivative is symmetric with respect to $i$, it is called the central difference. The same is true for Eq. (98), called the second-order central difference. There are other numerical approximations to numerical differentiation, each with its advantages and drawbacks [²⁸].

Substituting Eq. (98) into the $s$-wave zero-energy radial equation, we have

(99) $$u_{i+1}=2u_i-u_{i-1}+\frac{2m_r{(\mathrm{\Delta }r)}^2}{{\mathrm{\hslash }}^2}V(r_i)u_i.$$

This equation tells us that if we know the value of the radial solution for two consecutive points, at $r_{i-1}$ and$r_i$, we can calculate the value for the next point $u_{i+1}$.

We know that $u(0)=0$; hence we need to know $u(\mathrm{\Delta }r)$ to begin the calculation. For attractive potentials, $u(\mathrm{\Delta }r)$ has a non-zero value. If we want to find a solution without worrying about the normalization, we can set $u(\mathrm{\Delta }r)=1$. In other words, Schrödinger’s equation is linear: if some $u(r)$ is a solution, then $Cu(r)$ is also a solution, with $C$ constant.

Now we have everything to write a program that solves $u(r)$ inside the range $R$ of the potential $V(r)$. For reasons that will be apparent when we compute the scattering length, it is helpful to find the solution up to $R+\mathrm{\Delta }r$. The inputs are typically the number of points $N$ (or the discretization $\mathrm{\Delta }r$), and the parameters of the potential (in the case of the spherical well, Eq.70, $v_0$ and $R$). Once we know the number of points $N$, we may initialize an array of the desired dimension to hold the values of $u_i$. We also need to write a function that returns the value of $V(r_i)$ given a position $r_i$. Once these functionalities are implemented, the following algorithm will provide the values of $u_i$ for the whole interval:

1. Set $u_0=0$, $u_1=1$, and $i=1$.

2. Use Eq.(99) to compute $u_{i+1}$.

3. If $r_i⩾R+\mathrm{\Delta }r$, stop. Else, increment $i$ by one.

4. Go to step 2.

3.1.2. Numerov’s method

Equation (98) is one possible discretization for a numerical second derivative. There are other alternatives if we want to improve the precision of our algorithm. For example, Numerov’s method is a numerical technique capable of solving differential equations of second order when the first-order term is not present [²⁰]. It applies to equations of the form

(100) $$\frac{d^{2}y}{d{x}^2}=-\xi (x)y(x)+s(x).$$

The method provides a relation between $y_i\equiv y(x_i)$ at three equally spaced consecutive points ($x_{i-1}$, $x_i$, and $x_{i+1}$),

(101) $$y_{i+1}={\displaystyle \frac{1}{\left(1+\frac{{(\mathrm{\Delta }x)}^2}{12}{\xi }_{i+1}\right)}}\{2y_i(1-{\displaystyle \frac{5{(\mathrm{\Delta }x)}^2}{12}}{\xi }_i)$$ $$-y_{i-1}\left(1+{\displaystyle \frac{{(\mathrm{\Delta }x)}^2}{12}}{\xi }_{i-1}\right)$$ $$+{\displaystyle \frac{{(\mathrm{\Delta }x)}^2}{12}}(s_{i+1}+10s_i+s_{i-1})\}+𝒪[{(\mathrm{\Delta }x)}^6],$$

where $\mathrm{\Delta }x$ is the spacing between the points, ${\xi }_i\equiv \xi (x_i)$, and $s_i\equiv s(x_i)$.

The appeal of this method to our case is immediate. The $s$-wave zero-energy radial equation is of the form of Eq. (100) with $y\to u$, $x\to r$, $s=0$, and

(102) $$\xi (r)=-\frac{2m_r}{{\mathrm{\hslash }}^2}V(r).$$

The algorithm presented in Sec.^3.1.1 is mostly unchanged if we use Numerov’s method instead of the second-order central difference. The key change is that substituting Eq. (99) by (101) increases the precision without significant technical complications.

3.1.3. Dimensionless quantities

Schrödinger’s equation contains relatively small quantities. Planck’s reduced constant is $\mathrm{\hslash }\sim {10}^{-34}$ J s (or $\sim {10}^{-15}$ eV s), while the typical masses, length, and energy scales are also small. This does not affect analytical calculations since the numerical values are commonly substituted at the end of the procedure. However, computers deal with real numbers differently. In computing, floating-point representation is an approximation that allows real numbers to be stored using an integer with a fixed precision, called the significand, scaled by an integer exponent of a fixed base. The number of bits dedicated to the significand and exponent determines the accuracy and range of the floating-point numbers that can be represented. However, due to the intrinsic limitations of representing real numbers in this format, rounding errors and other inaccuracies can occur, affecting the accuracy of computations involving floating-point numbers [²⁸].

For these reasons, we would like to work within our program with quantities of the order of one (or close to it) so that we do not have to deal with very large or very small numbers. At the end of the calculation, once the solution is found, we may substitute the physical constants to compute the values for a particular system of interest. The procedure we describe here is equivalent to what is commonly called “set $\mathrm{\hslash }=m_r=1$”, but we provide a step-by-step derivation so that it is more clear how to recover the units at the end of the calculation.

Let us choose a length scale $\mathrm{\ell }$. The convenient value of $\mathrm{\ell }$ depends on the system under study; for atomic physics, it may be 1 Å; for nuclear physics, we may use 1 fm or any other length scale that makes sense for a particular problem. In this section, we denote dimensionless quantities with bars over them. Then, the scaled distance

(103) $$\overline{r}=\frac{r}{\mathrm{\ell }}$$

is dimensionless. The radial function $u(r)$ has units of [length]${}^{-1/2}$ (a straightforward way to see this is to write the integral corresponding to its normalization). Hence the transformation to make it dimensionless is

(104) $$\overline{u}(\overline{r})=\frac{u(r)}{{\mathrm{\ell }}^{-1/2}}.$$

Equation (103) implies that $d^2/d{r}^2=(1/{\mathrm{\ell }}^2)d^2/d{\overline{r}}^2$, so that the $s$-wave zero-energy radial equation becomes

(105) $$-\frac{{\mathrm{\hslash }}^2}{2m_r{\mathrm{\ell }}^2}\frac{d^2\overline{u}}{d{\overline{r}}^2}+V(\overline{r})\overline{u}=0.$$

The quantity

(106) $$ϵ=\frac{{\mathrm{\hslash }}^2}{m_r{\mathrm{\ell }}^2}$$

has units of energy, so we may make the potential dimensionless by considering $\overline{V}=V/ϵ$. Then, the equation we want to solve becomes

(107) $$-\frac{1}{2}\frac{d^2\overline{u}}{d{\overline{r}}^2}+\overline{V}(\overline{r})\overline{u}=0.$$

The net result is the same as “set $\mathrm{\hslash }=m_r=1$”, but now it is clear what to do. We implement and solve Eq. (107) in our program. After we obtained the results in this dimensionless fashion, we use Eqs. (103), (104), and (106) to recover the solutions with the desired units.

3.2. Scattering length and effective range

After following the numerical procedure outlined in Sec.^3.1, we should have the solution to the $s$-wave zero-energy radial equation for the desired potential. That is, we calculated $u_i$ in the interval $r=0$ to $R+\mathrm{\Delta }r$ at points $r_i$ equally spaced by $\mathrm{\Delta }r$.

Equation (44) is the low-energy ($k\to 0$) behavior of the radial solution outside the range of the potential, which we denote by $g_0(r)$ to avoid confusion. Then, we have the analytical expression,

(108) $$\begin{array}{cc}\frac{{g^{\prime }}_0(r)}{g_0(r)}=\frac{1}{r-a}& \text{for}r>R.\end{array}$$

On the other hand, we can compute the same quantity at $r=R$ using the numerical solution we obtained. The numerical derivative $u_{\mathrm{num}}^{\prime }$ can be calculated using Eq. (97),

(109) $$u_{\mathrm{num}}^{\prime }(R)={\frac{du(r)}{dr}|}_{r=R}=\frac{u(R+\mathrm{\Delta }r)-u(R-\mathrm{\Delta }r)}{2\mathrm{\Delta }r}.$$

Now it becomes apparent why we included the point $R+\mathrm{\Delta }r$ in the interval where we find $u(r)$ numerically. Finally, $u_{\mathrm{num}}^{\prime }(R)/u_{\mathrm{num}}(R)$ must match Eq. (108) at $r=R$ (their logarithmic derivatives must be equal at this point), which gives us the expression

(110) $$a=R-\frac{2\mathrm{\Delta }ru(R)}{u(R+\mathrm{\Delta }r)-u(R-\mathrm{\Delta }r)}.$$

This expression relates the scattering length and our numerical solution.

Equation (110) depends on the ratio of radial solutions, so we ignored its normalization. In other words, if we multiply all $u_i$ by a constant $C$, Eq. (110) is unaltered. However, as we saw in Sec.^2.3, the expression we derived for the effective range assumes a particular normalization choice. The constant $C$ is obtained by matching the value of the numerically-obtained radial function and the analytical solution at $r=R$, that is $Cu(R)=g(R)$,

(111) $$C=\frac{g(R)}{u(R)}=\frac{\left(1-R/a\right)}{u(R)}.$$

So, in the following, we assume that all the $u_i$ have been multiplied by this constant $C$.

Finally, what is left is to compute the integral that defines the effective range, Eq. (56), numerically. Although the $V=0$ radial solution, denoted by $g(r)$, has an analytical expression, it is useful to discretize it at the same points where we determined $u_i$ numerically, $g_i\equiv g(r_i)$. Then, the task becomes computing an integral of the form

(112) $$I={\int }_{x_1}^{x_N}f(x)dx$$

when the function $f(x)$ is known only at a discrete set of equally spaced points, $f(x_i)\equiv f_i$, where $i=1,2,3,\mathrm{\dots },N$. This is a well-known problem in numerical calculus, and many methods accomplish the task [²⁸]. One of the simplest methods is called the trapezoidal rule because it approximates the area under the curve by a trapezoid,

(113) $${\displaystyle {\int }_{x_1}^{x_N}}f(x)dx=\mathrm{\Delta }x\left[{\displaystyle \frac{1}{2}}{f}_1+f_2+\mathrm{\cdots }+f_{N-1}+{\displaystyle \frac{1}{2}}{f}_N\right]$$ $$+𝒪\left({\displaystyle \frac{{(x_N-x_1)}^3{f}^{\prime \prime }}{N^2}}\right).$$

A more sophisticated method is called Simpson’s rule, which considers a quadratic interpolation between the points,

(114) $${\displaystyle {\int }_{x_1}^{x_N}}f(x)dx=\mathrm{\Delta }x[{\displaystyle \frac{1}{3}}{f}_1+{\displaystyle \frac{4}{3}}{f}_2+{\displaystyle \frac{2}{3}}{f}_3+{\displaystyle \frac{4}{3}}{f}_4+\mathrm{\cdots }$$ $$+{\displaystyle \frac{2}{3}}{f}_{N-2}+{\displaystyle \frac{4}{3}}{f}_{N-1}+{\displaystyle \frac{1}{3}}{f}_N]$$ $$+𝒪\left({\displaystyle \frac{{(x_N-x_1)}^5{f}^{(4)}}{N^4}}\right).$$

Either method will successfully calculate the integral necessary for the effective range.

4. Examples

We chose four potentials to illustrate the numerical procedure described in Sec.³. The reasoning behind these choices is the following. The first one is the spherical well. Since we derived analytical expressions for its scattering length and effective range, Eqs. (80) and (92), this represents an excellent opportunity to test our program. We can readily compare our numerical solutions to the expected results, and it is much more straightforward to check the correctness of the program.

However, as we stated before, numerical computations would be unnecessary if all potentials had analytical expressions for their low-energy parameters. Hence, next, we considered the modified Pöschl-Teller potential. In this case, there is an analytical expression for the scattering length, but the effective range has a closed form only when $|a|\to \mathrm{\infty }$. Next, we considered a Gaussian potential since, in this situation, everything has to be computed numerically.

In the last three examples, we only considered purely attractive potentials. However, many interesting and relevant physical situations involve potentials where there is competition between a strong short-range repulsion and a weak long-range attraction. For this reason, we also considered the Lennard-Jones potential.

4.1. Spherical well

This potential has been covered extensively in the previous sections. To make comparisons easier with the other potentials, we rewrite Eq. (70) as

(115) $$V_{\text{sw}}(r)=\{\begin{array}{cc}-v_{\text{sw}}\frac{{\hslash }^2{\mu }_{\text{sw}}^2}{m_r},& \text{for}r<R,\\ 0,& \text{for}r>R,\end{array}$$

where we introduced the quantity ${\mu }_{\mathrm{sw}}=1/R$. In Fig.7, we compare its shape with other purely attractive potentials.

Figure 7
Spherical well, modified Pöschl-Teller, and Gaussian potentials in the unitary limit ($|a|\to \mathrm{\infty }$). We make the axis dimensionless by rescaling the distance in effective range units, and the potential in ${\mathrm{\hslash }}^2/(m_r\mathrm{}{r}_0^2)$ units, such that $\overline{V}\equiv m_r\mathrm{}{r}_0^2\mathrm{}V/{\mathrm{\hslash }}^2$ is dimensionless.

4.2. Modified Pöschl-Teller

The modified Pöschl-Teller potential has been used to describe systems in nuclear [⁹] and atomic physics [²⁹]. We present it in the form

(116) $$V_{\mathrm{PT}}(r)=-v_{\mathrm{PT}}\frac{{\mathrm{\hslash }}^2}{m_r}\frac{{\mu }_{\mathrm{PT}}^2}{{\cosh }^2({\mu }_{\mathrm{PT}}r)}.$$

Figure 7 illustrates its shape. This potential contains two parameters, $v_{\mathrm{PT}}$ and ${\mu }_{\mathrm{PT}}$, related to its depth and range. However, associating a range of $R$ to this potential is not as immediate as for the spherical well case. Fortunately, we only need a point $R$ where $V(R)\approx 0$. So, our numerical approach is to look for a value of $R$ such that the potential is negligible $|V(R)|⩽\epsilon$. Using dimensionless quantities, we found that for the mPT, Gaussian, and LJ potentials, $\epsilon ={10}^{-15}$ produces the desired results.

The eigenfunctions of the mPT potential may be obtained analytically [³⁰], but their derivation is beyond the scope of this work. However, if we perform the substitution $v_{\mathrm{PT}}=\lambda (\lambda -1)/2$, there is an analytical expression for the scattering length in terms of the parameters of the potential [³¹],

(117) $$a{\mu }_{\mathrm{PT}}=\frac{\pi }{2}\cot \left(\frac{\pi \lambda }{2}\right)+\gamma +\mathrm{\Psi }(\lambda ),$$

where $\gamma$ is the Euler-Mascheroni constant [³²] and $\mathrm{\Psi }$ is the digamma function [³²].

The $|a|\to \mathrm{\infty }$ case corresponds to $\lambda =2$ ($\cot (\pi )$ diverges) or $\lambda =-1$ ($\mathrm{\Psi }(-1)$ diverges), that is, $v_{\text{PT}}=1$. For this particular case, the $s$-wave zero-energy radial function takes a relatively simple form [³³],

(118) $$u_0(r)=\frac{\tanh ({\mu }_{\text{PT}}r)}{\tanh ({\mu }_{\text{PT}}R)}.$$

Taking its second derivative and substituting it in the radial equation confirms it is the appropriate solution. Then, we can readily calculate the effective range by integrating Eq. (56). In this case, $g_0(r)=1-r/a=1$, so that

(119) $$r_0=2{\displaystyle {\int }_0^R}dr\left[1-{\displaystyle \frac{{\tanh }^2({\mu }_{\text{PT}}r)}{{\tanh }^2({\mu }_{\text{PT}}R)}}\right]$$ $$=2\left[R-{\displaystyle \frac{R}{{\tanh }^2({\mu }_{\text{PT}}R)}}+{\displaystyle \frac{1}{{\mu }_{\text{PT}}\tanh ({\mu }_{\text{PT}}R)}}\right].$$

Since $1/{\mu }_{\mathrm{PT}}\sim R$ and the $\tanh (x)$ function converges rapidly to 1 as we increase $x$, we may set $\tanh (\mu R)=1$. Thus we have that $r_0=2/{\mu }_{\text{PT}}$ for $v_{\text{PT}}=1$ (unitarity).

4.3. Gaussian

Many interactions in physical systems are modeled through Gaussians, and their convenience and usefulness cannot be overstated. We write the Gaussian potential as

(120) $$V_{\mathrm{g}}(r)=-v_{\mathrm{g}}\frac{{\mathrm{\hslash }}^2}{m_r}{\mu }_{\mathrm{g}}^2{e}^{-r^2{\mu }_{\mathrm{g}}^2}.$$

The comparison of this potential with the previous two is illustrated in Fig.7. Although there are no analytical expressions for the scattering length and effective range in this case, numerical investigations in the literature provide benchmarks [³⁴].

4.4. Lennard-Jones

To include in our analysis a potential that also contains repulsion, we considered the Lennard-Jones potential,

(121) $$V_{\mathrm{LJ}}(r)=\frac{{\mathrm{\hslash }}^2}{m_r}\left[\frac{C_{12}}{r^{12}}-\frac{C_6}{r^6}\right],$$

where $C_{12}$ and $C_6$ are constants responsible for the strength of the repulsive and attractive interactions, with units of [length]${}^{10}$ and [length]${}^4$, respectively. In Fig.8 we plot this potential for $|a|\to \mathrm{\infty }$. The reader may be more familiarized with an alternative way of writing this potential,

Figure 8
Same as Fig.7, but for the Lennard-Jones potential. Differently from the other potentials, we observe a strong repulsive component near the origin.

(122) $$V_{\mathrm{LJ}}(r)=4{\epsilon }_{\mathrm{LJ}}\left[{\left(\frac{{\sigma }_{\mathrm{LJ}}}{r}\right)}^{12}-{\left(\frac{{\sigma }_{\mathrm{LJ}}}{r}\right)}^6\right],$$

but the conversion between the different constants is straightforward.

Besides the repulsive part, the LJ potential differs from the others considered so far because it diverges for $r\to 0$. Physically, the hard core near the origin prevents two particles from coming close together. This poses a numerical difficulty since $V_{\mathrm{LJ}}(0)$ diverges. The boundary condition $u(0)=0$ certainly helps, but the potential at the neighboring site, $V_{\mathrm{LJ}}(\mathrm{\Delta }r)$, is still large and may lead to numerical instabilities. Hence, the first step of the algorithm presented in Sec.^3.1.1 must be slightly changed to circumvent this problem. Since the potential is very large and positive near $r\approx 0$, it prevents the particles from coming close together, i.e., $u\approx 0$ in this region. For this reason, we define a range $0⩽r<r_{\text{min}}$ where $u(r)=0$ and start our integration at $r=r_{\text{min}}$. After this is done, we may proceed as outlined in Sec.^3.1.1. Using dimensionless units, we found that taking $r_{\text{min}}$ such that $V(r_{\text{min}})\approx {10}^{10}$ produces the desired results.

4.5. Tuning the parameters

The four potentials we presented have two parameters. In the case of the purely attractive ones, one parameter is associated with the depth of the potential ($v_{\mathrm{sw}}$, $v_{\mathrm{PT}}$, and $v_{\mathrm{g}}$) and another with its range (${\mu }_{\mathrm{sw}}$, ${\mu }_{\mathrm{PT}}$, and${\mu }_{\mathrm{g}}$). In contrast, for the LJ potential, they control the attractive ($C_6$) and repulsive ($C_{12}$) components. The numerical procedure we described so far works if the two parameters are known. For example, they could be listed in the literature. However, the situation is often different: the scattering length and effective range are known, and we want to tune the parameters of a particular potential to reproduce the desired $a$ and $r_0$ values. Since we want to match two values and have two free parameters, the correspondence is one-to-one (with the restriction of how many bound states we want).

The Sec.^3.1.1 algorithm can still be used in this case, but we need to apply it repeatedly as we vary the two parameters. Let us denote the values of the two parameters by $v_1$ and $v_2$. The first step is to start with a guess for $v_1$ and $v_2$ and then perform the procedure as described in Sec.^3.1.1 to obtain the scattering length and effective range. They will probably not be the desired values, but it is important to check the radial wave function to see if at least the number of bound states is what we expect. Typically, we want no bound states for $a<0$ and one for $a>0$, in other words, a nodeless $u(r)$ for the first and one node for the latter, as depicted in Fig.3. If the initial guess contains more bound states than desired, then we need to decrease the potential depth until we reach the target number. Then, a possible procedure is:

1. Start with a guess $(v_1,v_2)$.

2. Compute $a$ and $r_0$ as in Sec.^3.1.1.

3. Keep $v_2$ fixed. Vary $v_1$ until $a$ has the desired value. Increasing the depth of the potential will decrease the value of the scattering length (until it diverges and changes from $-\mathrm{\infty }$ to $+\mathrm{\infty }$).

4. Keep $v_1$ fixed at the value found in step 3. Vary $v_2$ until $r_0$ has the desired value. Increasing the range of the potential will increase $r_0$.

5. If $a$ and $r_0$ match the desired values, stop. Else, go to step 3 and use the value of $v_2$ found in step 4.

We may implement this algorithm using two nested loops: one varies $v_1$, and the other runs $v_2$. It is helpful to define a tolerance for the numerical values of $a^{\mathrm{num}}$ and $r_0^{\mathrm{num}}$, such that $|a^{\mathrm{num}}-a|$ and $|r_0^{\mathrm{num}}-r_0|$ are small quantities. After the procedure is done, checking the number of nodes of the radial function is necessary to guarantee that it is indeed the situation we wanted to reproduce.

4.6. Results

Finally, applying the numerical procedure described so far to some concrete examples is instructive. The objective is twofold: we want to illustrate the method and provide some results so the reader can use them to check their program. We present 3 cases: $a<0$, $|a|\to \mathrm{\infty }$, and $a>0$, which correspond to three very distinct physical situations.

As we pointed out in Sec.^2.6.2, the universality in the low-energy scattering results is not specific to a single physics area so we can choose from several instances. We picked nuclear physics since two nucleon systems are excellent for illustrating low-energy scattering. There is no bound state in the region where $a$ is negative. That is the case of two neutrons, with a scattering length of $a=-18.5$ fm and an effective range of $r_0=2.7$ fm [⁹]. The region where $a>0$ is interesting due to the appearance of the first bound state of the system, where lies the deuteron, a proton and a neutron, with $a=5.4$ fm and $r_0=1.7$ fm. The region where $|a|\to \mathrm{\infty }$, often called unitary limit, is of great interest [⁸]. It does not occur naturally in nuclear physics, but we can still study its properties. Table2 summarizes the physical situations we cover in this section.

Table 2.
Scattering length and effective range values that we want to reproduce with four forms of two-body potentials.

System	$a$ (fm)	$r_0$ (fm)
Neutron-neutron	$-18.5$	$2.7$
Unitarity	$\pm \mathrm{\infty }$	$1.0$
Deuteron	$5.4$	$1.7$

As a side note, it is worth mentioning that most systems have their interactions fixed by what is observed in nature, such as the case of two neutrons, the deuteron, and countless other examples. Although it is possible to theoretically investigate the continuous variation of the scattering length from negative to positive values, we would like to see this happening in a physical system. In experiments with cold atoms, this is possible. A clever method called Feshbach resonances [⁷] uses the manipulation of magnetic fields to vary the scattering length of the potential between two atoms. This allows the study of systems with scattering lengths of both signs and the unitarity limit between them.

We performed the procedure described in Sec.^4.5 to obtain the parameters, scattering lengths, and effective ranges presented in Tables3 and4. It is worth remembering that we can check the results for the spherical well against the analytical expressions, Eqs. (80) and (92). For the mPT, Eq. (117) provides a benchmark for the scattering length, and in Sec.^4.2 we show that $r_0=2/{\mu }_{\mathrm{PT}}$ when $|a|\to \mathrm{\infty }$.

Table 3.
Parameters for the spherical well, modified Pöschl-Teller, and Gaussian potentials that reproduce the desired scattering lengths and effective range.

Potential	$v$	$\mu$ (fm${}^{-1}$)	$a$ (fm)	$r_0$ (fm)
Neutron-neutron
Well	$1.1096$	$0.3918$	$-18.52$	$2.7$
mPT	$0.9071$	$0.7991$	$-18.51$	$2.7$
Gaussian	$1.2121$	$0.5672$	$-18.55$	$2.7$
Unitarity
Well	$1.2337$	$1.0000$	$\sim -{10}^5$	$1.0$
mPT	$1.0000$	$2.0000$	$\sim {10}^9$	$1.0$
Gaussian	$1.3420$	$1.4349$	$\sim -{10}^5$	$1.0$
Deuteron
Well	$1.7575$	$0.5000$	$5.4$	$1.70$
mPT	$1.4388$	$0.8631$	$5.4$	$1.73$
Gaussian	$1.9102$	$0.6754$	$5.4$	$1.70$

Table 4.
Numerical results for the Lennard-Jones potential.

System	$C_{12}$ (fm${}^{10}$)	$C_6$ (fm${}^4$)	$a$ (fm)	$r_0$ (fm)
Neutron-neutron	$3.08836698$	$9.86668911$	$-18.5$	$2.71$
Unitarity	$0.00034068$	$0.26462461$	$\sim -{10}^5$	$1.00$
Deuteron	$0.90485319$	$6.81472000$	$5.4$	$1.70$

The two-body potentials are very different in shape, as shown in Figs.7 and8. The fact that we can tune them to produce the same scattering length $a$ and effective range $r_0$, although remarkable, is directly related to the shape-independent approximation, Eq. (55). As the name suggests, this equation describes the $s$-wave phase shift ${\delta }_0(k)$ as independent of the shape of the scattering potential.

Besides the values of $a$ and $r_0$, the numerical procedure produces the reduced radial wave functions $u(r)$. In Fig.9, we compare $u(r)$ solved with the Numerov’s method for the spherical well, modified Pöschl-Teller, Gaussian, and Lennard-Jones potentials. They are in excellent agreement with analytical solutions, Eq. (90) for the spherical well, and Eq. (118) for the mPT at unitarity. For the same physical system, the functions differ inside the range of the potentials but are all the same outside, $1-r/a$.

Figure 9
Reduced radial solutions, $u(r)=rA\mathrm{}(r)$, of Schrödinger’s equation with the proposed potentials for three different situations. The dashed curves correspond to the analytical solutions. Outside the range of the potential, large$r$, all solutions must be the same.

It is also possible to see the geometric interpretation of the scattering length alluded to in Figs.3a and3b. For $a<0$, Fig.9a, $u(r)$ does not intercept the $r$-axis, and the scattering length is negative because the extrapolation of the reduced radial function intercepts the axis at a negative value. In Fig.9c, $u(r)$ intercepts the $r$-axis at the scattering length, thus $a>0$. Neither the reduced radial function nor its extrapolation intercept the axis at unitarity, since $|a|\to \mathrm{\infty }$, Fig.9b.

Finally, after obtaining the parameters for these three illustrative physical systems, it is interesting to investigate how the scattering length varies as we continuously vary the parameters of the potentials. In Fig.10a, we present the ratio $a/r_0$ as a function of the potential depth of the spherical well, modified Pöschl-Teller, and Gaussian potentials ($v_{\mathrm{sw}}$, $v_{\mathrm{PT}}$, and $v_{\mathrm{g}}$). Analytical expressions are only available for the first two, and it is possible to see that our numerical investigations are in excellent agreement with them.

Figure 10
Scattering length behavior as a function of the interaction strength for the spherical well, modified Pöschl-Teller, Gaussian, and Lennard-Jones potentials. The dashed curves in (a) are the analytical solutions for the spherical well, Eq. (80), and for the modified Pöschl-Teller potential, Eq. (117).

To produce a similar figure for the Lennard-Jones case, we fixed the value of $C_{12}/r_0^{10}=0.000341$ and varied$C_6$. We present the results in Fig.10b. Although the behavior is qualitatively similar to the purely attractive potentials, the main difference is the region close to $C_6\approx 0$ where the scattering length is positive. This is because, for $C_6=0$, the Lennard-Jones potential is purely repulsive. Then the scattering length is positive not due to a bound state but because the repulsive part “pushes” the wave function away from the origin, as depicted in Fig.3c.

5. Conclusions

This work presented quantum scattering theory fundamentals focusing on the low-energy limit. In this context, we aimed to introduce two significant quantities: the scattering length and the effective range. To illustrate how these two parameters behave in a concrete example, we derived analytical expressions for both in the case of the spherical well. We also showed how the energy of a bound state could be calculated using zero- and finite-range expressions applied to a ${}^4$He dimer and the deuteron.

After introducing the main results of low-energy scattering theory, we proceeded with our main goal: describing a numerical procedure that can be used to compute the scattering length and effective range of any spherically symmetric finite-ranged two-body potential. We chose the spherical well, modified Pöschl-Teller, Gaussian, and Lennard-Jones potentials as examples. We then performed calculations related to three physical systems: two neutrons, unitarity, and the deuteron.

This work shows that these calculations are relevant and can be readily applied in atomic and nuclear physics. Hopefully, students can follow the procedure outlined in this manuscript, reproduce our results, extend them to their choice of physical systems, and apply the method to other potentials.

Acknowledgments

This work was partially supported by the São Paulo Research Foundation (FAPESP) grant 2018/09191-7 (L.M.), and by the National Council for Scientific and Technological Development (CNPq) grants 383501/2022-9 (L.M.) and 131025/2022-8 (M.M.L.).

Appendix A: Scattering theory

Appendix A.1: Time-dependent formalism

We may treat interactions more simply in a formalism called interaction picture. We denote the time evolution of a state ket at an initial time $t_0$ in the Schrödinger picture as

(A1) $$\text{|}\varphi {(t\mathrm{)\rangle }}_S=U_S(t,t_0)\text{|}\varphi {(t_0\mathrm{)\rangle }}_S,$$

where $|⟩{}_S$ is a ket in the Schrödinger picture and $U_S(t_0,t)$ is the time-evolution operator. If $H$ is time independent, then $U_S(t_0,t)=e^{-iH(t-t_0)/\mathrm{\hslash }}$. The interaction-picture state ket, $|⟩{}_I$, is defined as

(A2) $$\text{|}\varphi {(t\mathrm{)\rangle }}_I=e^{i{H}_{0}t/\mathrm{\hslash }}\text{|}\varphi {(t\mathrm{)\rangle }}_S.$$

The operators are defined as

(A3) $$A_I=e^{i{H}_{0}t/\mathrm{\hslash }}{A}_S{e}^{-i{H}_{0}t/\mathrm{\hslash }}.$$

The Schrödinger equation $i\mathrm{\hslash }\frac{\partial }{\partial t}\text{|}\varphi {(t\mathrm{)\rangle }}_S=H\text{|}\varphi {(t\mathrm{)\rangle }}_S$ may be written as

(A4) $$i\mathrm{\hslash }\frac{\partial }{\partial t}\text{|}\varphi {(t\mathrm{)\rangle }}_I=V_I\text{|}\varphi {(t\mathrm{)\rangle }}_I,$$

where $V_I=e^{i{H}_0(t-t_0)/\mathrm{\hslash }}V{e}^{-i{H}_0(t-t_0)/\mathrm{\hslash }}$ is the potential in the interaction picture. Equation (A4) is a Schrödinger-like equation. The advantage is that we remove $H_0$ from our calculations to consider the interaction. If $V=0$, $\text{|}\varphi {(t\mathrm{)\rangle }}_I$ is constant in time (and equal to $\text{|}\varphi {(t_0\mathrm{)\rangle }}_S$). Similar to Eq.A1, $\text{|}\varphi {(t\mathrm{)\rangle }}_I$ evolves in time as

(A5) $$\text{|}\varphi {(t\mathrm{)\rangle }}_I=U_I(t,t_0)\text{|}\varphi {(t_0\mathrm{)\rangle }}_I,$$

where $U_I(t,t_0)$ is defined as

(A6) $$U_I(t,t_0)=e^{i{H}_{0}t/\mathrm{\hslash }}{U}_S(t,t_0)e^{-i{H}_0{t}_0/\mathrm{\hslash }}$$

and obeys the Schrödinger-like equation

(A7) $$i\mathrm{\hslash }\frac{\partial }{\partial t}{U}_I(t,t_0)=V_I(t)U_I(t,t_0).$$

The solution of Eq.A7 may be formally written as

(A8) $$U_I(t,t_0)=1-\frac{i}{\mathrm{\hslash }}{\int }_{t_0}^{t}d{t}^{\prime }{V}_I(t^{\prime })U_I(t^{\prime },t_0),$$

for the required initial condition $U_I(t_0,t_0)=1$. Our goal is to calculate the evolution of the state in a distant past $t_0\to -\mathrm{\infty }$, that is, $\text{|}\varphi (t\to -\mathrm{\infty }\mathrm{)\rangle }=\text{|}i\mathrm{\rangle }$. Nevertheless, equation (A8) is only valid for finite times $t$ and $t_0$ and so setting $U_I(t,-\mathrm{\infty })$ would lead to convergence problems. To avoid it, we will give meaning to $U_I(t,-\mathrm{\infty })$ by writing it as [³⁵] (proof in Appendix ^B),

(A9) $$U_I(t,-\mathrm{\infty })=\underset{t_0\to -\mathrm{\infty }}{lim}{U}_I(t,t_0)=\underset{ϵ\to 0}{lim}ϵ{\int }_{-\mathrm{\infty }}^{0}d{t}^{\prime }{e}^{ϵt^{\prime }}{U}_I(t,t^{\prime }),$$

where $U_I(t,t^{\prime })$ is defined as in eq. (A6). Despite our time-dependent treatment, we recall that $H$ is time-independent. Thus $U_S(t,t^{\prime })=e^{-iH(t-t^{\prime })/\mathrm{\hslash }}$. The state vector at a time $t=0$ is given by

(A10) $$\text{|}\varphi {(t=0\mathrm{)\rangle }}_I=U_I(0,-\mathrm{\infty })\text{|}i\mathrm{\rangle },$$

where

(A11) $$U_I(0,-\mathrm{\infty })=\underset{ϵ\to 0}{lim}ϵ{\int }_{-\mathrm{\infty }}^{0}d{t}^{\prime }{e}^{ϵt^{\prime }}{e}^{iH{t}^{\prime }/\mathrm{\hslash }}{e}^{-i{H}_0{t}^{\prime }/\mathrm{\hslash }}.$$

Now applying it to equation (A10):

(A12) $$\text{|}\varphi {(t=0\mathrm{)\rangle }}_I=\underset{ϵ\to 0}{lim}ϵ{\displaystyle {\int }_{-\mathrm{\infty }}^0}d{t}^{\prime }{e}^{ϵt^{\prime }}{e}^{i(H-E_i)t^{\prime }/\mathrm{\hslash }}\text{|}i\mathrm{\rangle }$$ $$=\underset{ϵ\to 0}{lim}{\displaystyle \frac{iϵ}{E_i-H+iϵ}}\text{|}i\mathrm{\rangle }.$$

Using the identity (proof in ^{Appendix C})

(A13) $${\displaystyle \frac{1}{E_i-H+iϵ}}-{\displaystyle \frac{1}{E_i-H_0+iϵ}}$$ $$={\displaystyle \frac{1}{E_i-H_0+iϵ}}V{\displaystyle \frac{1}{E_i-H+iϵ}},$$

we rewrite Eq. (A12) as

(A14) $$\text{|}\varphi {(t=0\mathrm{)\rangle }}_I=\underset{ϵ\to 0}{lim}{\displaystyle \frac{iϵ}{E_i-H_0+iϵ}}\text{|}i\mathrm{\rangle }$$ $$+{\displaystyle \frac{1}{E_i-H_0+iϵ}}V{\displaystyle \frac{iϵ}{E_i-H+iϵ}}\text{|}i\mathrm{\rangle }$$ $$=\underset{ϵ\to 0}{lim}{\displaystyle \frac{iϵ}{E_i-H_0+iϵ}}\text{|}i\mathrm{\rangle }$$ $$+{\displaystyle \frac{1}{E_i-H_0+iϵ}}V\text{|}\varphi {(t=0\mathrm{)\rangle }}_I$$

Since $\text{|}i$〉 is an eigenfunction of $H_0$ (it satisfies $H_0\text{|}i\mathrm{\rangle }=E_i\text{|}i\mathrm{\rangle }$), the first term in the RHS is the identity operator. Then we can write

(A15) $$\text{|}\psi \mathrm{\rangle }=\text{|}i\mathrm{\rangle }+\frac{1}{E_i-H_0+iϵ}V\text{|}\psi \mathrm{\rangle },$$

where we left off the notation $\text{|}\varphi (t=0)$〉 to emphasize that this is an actual time-independent problem. This is know as the Lippmann-Schwinger equation.

Back to the time-dependent formulation, the state vector at any time $t$ could be obtained by the sequential relation of the time translation operator: $U_I(t,t_0)=U_I(t,t^{\prime })U_I(t^{\prime },t_0)$, so

(A16) $$\text{|}\varphi (t\mathrm{)\rangle }=U_I(t,0\mathrm{)\rangle }\text{|}\varphi (0)=U_I(t,-\mathrm{\infty })\text{|}i\mathrm{\rangle }.$$

In a distant future, that is, setting $t\to +\mathrm{\infty }$, we write

(A17) $$\text{|}f\mathrm{\rangle }=S\text{|}i\mathrm{\rangle },$$

where $\text{|}\varphi (t\to +\mathrm{\infty })$〉 is our asymptotic final state $\text{|}f$ 〉 and the scattering matrix (or $S$ matrix) is defined as

(A18) $$S\equiv U_I(+\mathrm{\infty },-\mathrm{\infty }).$$

Equation (A17) asserts that the action of the $S$ matrix on an initial state (that exists asymptotically for $t_0\to -\mathrm{\infty }$) transforms the ket $\text{|}i$〉 into a final state that exists in a distant future $t\to +\mathrm{\infty }$. That is, some sort of scattering occurred between $t$ and $t_0$. Furthermore, Eq. (A15) may be written as the power series expansion

$$\text{|}\psi \mathrm{\rangle }=\text{|}i\mathrm{\rangle }+G_{+}V\text{|}i\mathrm{\rangle }+G_{+}V{G}_{+}V\text{|}i\mathrm{\rangle }+\mathrm{\dots }$$ $$=\text{|}i\mathrm{\rangle }+G_{+}\left(V+V{G}_{+}V+\mathrm{\dots }\right)\text{|}i\mathrm{\rangle },$$

where $G_{+}\equiv {(E_i-H_0+iϵ)}^{-1}$.

We define the transition matrix $T$ as the perturbative series

(A19) $$T\equiv V+V{G}_{+}V+V{G}_{+}V{G}_{+}V+\mathrm{\dots },$$

which is often called the Dyson series. We may write the scattering state as

(A20) $$\text{|}\psi \mathrm{\rangle }=\text{|}i\mathrm{\rangle }+G_{+}T\text{|}i\mathrm{\rangle }.$$

The consequence is that

(A21) $$V\text{|}\psi \mathrm{\rangle }=T\text{|}i\mathrm{\rangle }.$$

The perturbative character of Eq. (A19) relates the $T$-matrix to being a species of a generalized potential, where in a first-order perturbation $T$ and $V$ are equivalent and thus $\text{|}\psi \mathrm{\rangle }\approx \text{|}i$〉. This is known as the first-order Born approximation [³].

Appendix A.2: Scattering theory integral equations

Now we focus on the physical implications of Eq. (A15). Writing $\text{|}\psi$〉 in the position basis and inserting a complete set of position states between $G_{+}$ and $V$ leads to the integral equation

(A22) $$\mathrm{\langle }\mathbf{r}\mathrm{|}\psi \mathrm{\rangle }=\mathrm{\langle }\mathbf{r}\text{|}i\mathrm{\rangle }+\int d^3{\mathbf{r}}^{\mathbf{\prime }}\mathrm{\langle }\mathbf{r}\text{|}{G}_{+}\text{|}{\mathbf{r}}^{\prime }\mathrm{\rangle }\mathrm{\langle }{\mathbf{r}}^{\prime }\mathrm{|}V\mathrm{|}\psi \mathrm{\rangle },$$

where we must calculate the function

(A23) $$G_{+}(\mathbf{r},{\mathbf{r}}^{\prime })\equiv \frac{{\mathrm{\hslash }}^2}{2m}⟨\mathbf{r}\left|\frac{1}{E-H_0+iϵ}\right|{\mathbf{r}}^{\prime }⟩.$$

Recalling that the momentum basis $\{\text{|}k\mathrm{\rangle \}}$ set elements are eigenstates of $H_0$ with eigenvalues ${\mathrm{\hslash }}^2{\mathbf{k}}^2/2m$, Eq. (1), we may insert this discrete basis by writing

(A24) $$G_{+}(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{{\mathrm{\hslash }}^2}{2m}\sum _{{\mathbf{k}}^{\prime },{\mathbf{k}}^{\prime \prime }}\mathrm{\langle }\mathbf{r}\mathrm{|}{\mathbf{k}}^{\prime }\mathrm{\rangle }⟨{\mathbf{k}}^{\prime }\left|\frac{1}{E-H_0+iϵ}\right|{\mathbf{k}}^{\prime \prime }⟩\mathrm{\langle }{\mathbf{k}}^{\prime \prime }\mathrm{|}{\mathbf{r}}^{\prime }\mathrm{\rangle }.$$

The quantized plane waves in the position representation are written as

(A25) $$\mathbf{r}\mathrm{|}\mathbf{k}\mathrm{\rangle }=\frac{e^{i\mathbf{k}\cdot \mathbf{r}}}{L^{3/2}}\text{and}\mathrm{\langle }\mathbf{k}\mathrm{|}\mathbf{r}\mathrm{\rangle }=\frac{e^{-i\mathbf{k}\cdot \mathbf{r}}}{L^{3/2}}.$$

Now by letting $H_0$ act on $\mathrm{\langle }{\mathbf{k}}^{\prime }$| and noting that $\mathrm{\rangle }{\mathbf{k}}^{\prime }{\mathbf{k}}^{\prime \prime }={\delta }_{{\mathbf{k}}^{\prime }{\mathbf{k}}^{\prime \prime }}$ and $E={\mathrm{\hslash }}^2{k}^2/2m$,

(A26) $$G_{+}(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{L^3}\sum _{{\mathbf{k}}^{\prime }}\frac{e^{i{\mathbf{k}}^{\prime }\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })}}{k^2-k^{\prime 2}+iϵ},$$

where we absorbed the factor $2m/{\mathrm{\hslash }}^2$ into $ϵ$.

We are left with a sum in $k$-space, a discrete space. Due to periodic boundary conditions, each $k_i$ ($i=x,y,z$) is located at $k_i=2\pi n_i/L$, where $n_i=0,1,2,3\mathrm{\dots }$. That is, $\{k\}=\{0,2\pi /L,4\pi /L,\mathrm{\dots }\}$. The distance between each point in the $k$-space is $\mathrm{\Delta }k=2\pi /L$. We must take $L\to \mathrm{\infty }$ to guarantee the required continuous character. For this reason, the separation between points is very small compared to $L$ ($\mathrm{\Delta }k\approx 0$). The summation in $\mathbf{k}$ may be taken to be over a continuous space. We may substitute the sum with an integral,

(A27) $$\sum _{{\mathbf{k}}^{\prime }}\to \int \rho (k)d^3{\mathbf{k}}^{\prime }=\frac{L^3}{{(2\pi )}^3}\int d^3{\mathbf{k}}^{\prime },$$

where the factor $\rho (k)=L^3/{(2\pi )}^3$ is the $k$-density in 3 dimensions. The function becomes

(A28) $$G_{+}(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{{(2\pi )}^3}\int d^3{\mathbf{k}}^{\prime }\frac{e^{i{\mathbf{k}}^{\prime }\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })}}{k^2-k^{\prime 2}+iϵ}.$$

We may perform this integration by choosing spherical coordinates $(k^{\prime },\theta ,\phi )$ and, without loss of generalization, we may let the vector $(\mathbf{r}-{\mathbf{r}}^{\prime })$ lie along the $k_z^{\prime }$ axis so that ${\mathbf{k}}^{\prime }\cdot (\mathbf{r}-{\mathbf{r}}^{\prime })=k^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|\cos \theta$. We then write

(A29) $$G_{+}(\mathbf{r},{\mathbf{r}}^{\prime })={\displaystyle \frac{1}{{(2\pi )}^2}}{\displaystyle {\int }_0^{\mathrm{\infty }}}d{k}^{\prime }{k}^{\prime 2}{\displaystyle {\int }_0^{\pi }}d\theta {\displaystyle \frac{e^{i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|\cos \theta }}{k^2-k^{\prime 2}+iϵ}}$$ $$={\displaystyle \frac{1}{4{\pi }^2|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}d{k}^{\prime }{\displaystyle \frac{k^{\prime }\sin (k^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|)}{k^2-k^{\prime 2}+iϵ}},$$

where we note that the integrand is even. We also note a pole in $k^{\prime }=\pm \sqrt{k^2+iϵ}$. It would make it ambiguous if we chose to take $ϵ\to 0$ before evaluating the integral. Luckily, the pole is not on the real axis as $\pm \sqrt{k^2+iϵ}\approx k+iϵ/2k+{ϵ}^2/8k^3+\mathrm{\dots }$, which gives meaning to the integration process. Ignoring terms higher than ${ϵ}^2$ and redefining $ϵ/2k\to ϵ$, we may factorize $(k^2-k^{\prime 2}+iϵ)=-(k^{\prime }-k-iϵ)(k^{\prime }+k+iϵ)$. The integral is then

(A30) $$G_{+}(\mathbf{r},{\mathbf{r}}^{\prime })={\displaystyle \frac{1}{8i{\pi }^2|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}d{k}^{\prime }{\displaystyle \frac{k^{\prime }(e^{-i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}-e^{i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|})}{(k^{\prime }-k-iϵ)(k^{\prime }+k+iϵ)}}.$$

We should let $k^{\prime }$ momentarily be a complex variable to carry out this integration and use the residue theorem [²¹]. Consider the paths in Fig.11. The closed path integral may be written as the sum

Figure 11
Paths to calculate the integral (A11). We must choose the upper plane path ${\mathrm{\Gamma }}_R^{+}$ for $e^{+i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}$ because $e^{+i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}\to 0$ as $\mathrm{Im}{k}^{\prime }$ takes $+\mathrm{\infty }$ values. Similarly, we choose the lower plane path ${\mathrm{\Gamma }}_R^{-}$ for $e^{-i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}$ because $e^{-i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}\to 0$ as $\mathrm{Im}{k}^{\prime }$ takes $-\mathrm{\infty }$ values. This makes the integration in ${\gamma }_R^{\pm }$ go to $0$, and we are left with the path $C_R$, which lies in the axis $\mathrm{Re}{k}^{\prime }$. That is, the closed path ${\mathrm{\Gamma }}_R^{\pm }$ is equal to the real integral (A30).

(A31) $${\oint }_{{\mathrm{\Gamma }}_R}={\int }_{{\gamma }_R}+{\int }_{C_R}=2\pi i\times \sum _j\text{Res}\{k^{\prime };j\},$$

where the integral in the closed contour ${\mathrm{\Gamma }}_R$ is $2\pi i$ times the sum of the residues $\text{Res}\{k^{\prime };j\}$ due to poles $j$ enclosed by the curve, and the integral through the path ${\gamma }_R$ is zero due to Jordan’s lemma. The only pole inside ${\mathrm{\Gamma }}_R^{\pm }$ is $k\pm iϵ$. The residues may be calculated as follows

(A32) $$\text{Res}\{k^{\prime };k+iϵ\}$$ $$=\underset{\begin{array}{c}\hfill k^{\prime }\to k+iϵ\hfill \\ \hfill ϵ\to 0\hfill \end{array}}{lim}(k^{\prime }-k-iϵ){\displaystyle \frac{k^{\prime }{e}^{i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{(k^{\prime }-k-iϵ)(k^{\prime }+k+iϵ)}}$$ $$={\displaystyle \frac{e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{2}},$$

(A33) $$\text{Res}\{k^{\prime };-k-iϵ\}$$ $$=\underset{\begin{array}{c}\hfill k^{\prime }\to -k-iϵ\hfill \\ \hfill ϵ\to 0\hfill \end{array}}{lim}(k^{\prime }+k+iϵ){\displaystyle \frac{k^{\prime }{e}^{-i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{(k^{\prime }-k-iϵ)(k^{\prime }+k+iϵ)}}$$ $$={\displaystyle \frac{e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{2}}.$$

The integral through the real axis $C_R$ is

(A34) $$\underset{R\to \mathrm{\infty }}{lim}{\int }_{-R}^{R}d{k}^{\prime }\frac{k^{\prime }(e^{-i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|}-e^{i{k}^{\prime }|\mathbf{r}-{\mathbf{r}}^{\prime }|})}{(k^{\prime }-k-iϵ)(k^{\prime }+k+iϵ)}=2\pi i\times e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

Thus

(A35) $$G_{+}(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{4\pi }\frac{e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}.$$

Returning to Eq. (A22):

(A36) $$\mathrm{\langle }\mathbf{r}\mathrm{|}\psi \mathrm{\rangle }=\mathrm{\langle }\mathbf{r}\mathrm{|}i\mathrm{\rangle }-\frac{2m}{{\mathrm{\hslash }}^2}\int d^3{\mathbf{r}}^{\prime }\frac{1}{4\pi }\frac{e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}\mathrm{\langle }{\mathbf{r}}^{\prime }\mathrm{|}V\mathrm{|}\psi \mathrm{\rangle }.$$

Considering the potential $V({\mathbf{r}}^{\prime })$ to be local, that is,

(A37) $$\mathrm{\langle }{\mathbf{r}}^{\prime }\mathrm{|}V\mathrm{|}{\mathbf{r}}^{\prime \prime }\mathrm{\rangle }=V({\mathbf{r}}^{\prime })\delta ({\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }).$$

This allows us to write

(A38) $$\mathrm{\langle }\mathbf{r}\mathrm{|}\psi \mathrm{\rangle }=\mathrm{\langle }\mathbf{r}\mathrm{|}i\mathrm{\rangle }-\frac{2m}{{\mathrm{\hslash }}^2}\int d^3{\mathbf{r}}^{\prime }\frac{1}{4\pi }\frac{e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}V({\mathbf{r}}^{\prime })\mathrm{\langle }{\mathbf{r}}^{\prime }\mathrm{|}\psi \mathrm{\rangle }.$$

We now restrict ourselves to cases where the potential is of a finite range. This condition is important because the scattering is observed far away from the scattering center. We are interested in measurements at large distances $|\mathbf{r}|\gg |{\mathbf{r}}^{\prime }|$. In this regime,

(A39) $$e^{ik|\mathbf{r}-{\mathbf{r}}^{\prime }|}\approx e^{ikr}{e}^{-i{\mathbf{k}}^{\prime }\cdot {\mathbf{r}}^{\prime }},$$

where $r\equiv |\mathbf{r}|$. Furthermore, we now specify that our initial state is $\text{|}i\mathrm{\rangle }=\text{|}\mathbf{k}\mathrm{\rangle }$ and $\mathrm{\langle }\mathbf{r}\mathrm{|}\mathbf{k}\mathrm{\rangle }=e^{i\mathbf{k}\cdot \mathbf{r}}/L^{3/2}$. Finally, we arrive at the form

(A40) $$\psi (\mathbf{r},\theta )\stackrel{\text{large}r}{\to }\frac{1}{L^{3/2}}\left[e^{i\mathbf{k}\cdot \mathbf{r}}+\frac{e^{ikr}}{r}f({\mathbf{k}}^{\prime },\mathbf{k})\right],$$

where

(A41) $$f({\mathbf{k}}^{\prime },\mathbf{k})=-\frac{m{L}^3}{2\pi {\mathrm{\hslash }}^2}\int d^3{x}^{\prime }\mathrm{\langle }{\mathbf{k}}^{\prime }\mathrm{|}{\mathbf{r}}^{\prime }\mathrm{\rangle }V({\mathbf{r}}^{\prime })\mathrm{\langle }{\mathbf{r}}^{\prime }\mathrm{|}\psi \mathrm{\rangle }$$

is referred to as the scattering amplitude.

Appendix B:

We follow a procedure present in Ref. [³⁶]. Consider a well-behaved function $f(t)$ such that

(A42) $$\underset{t_0\to -\mathrm{\infty }}{lim}f(t_0)=f(-\mathrm{\infty }).$$

Now consider

(A43) $$\underset{ϵ\to 0}{lim}ϵ{\int }_{-\mathrm{\infty }}^{0}d{t}^{\prime }{e}^{ϵt^{\prime }}f(t^{\prime }),$$

where $ϵ\ll 1$. Integrating by parts leads to

(A44) $$\underset{ϵ\to 0}{lim}ϵ{\displaystyle {\int }_{-\mathrm{\infty }}^0}d{t}^{\prime }{e}^{ϵt^{\prime }}f(t^{\prime })=$$ $${\left[f(t^{\prime })e^{ϵt^{\prime }}\right]}_{-\mathrm{\infty }}^0-\underset{ϵ\to 0}{lim}{\displaystyle {\int }_{-\mathrm{\infty }}^0}d{t}^{\prime }{\displaystyle \frac{df}{d{t}^{\prime }}}{e}^{ϵt^{\prime }}$$ $$=f(0)-{\displaystyle {\int }_{-\mathrm{\infty }}^0}d{t}^{\prime }{\displaystyle \frac{df}{d{t}^{\prime }}}=f(-\mathrm{\infty }).$$

Thus, we have the equality

(A45) $$\underset{t_0\to -\mathrm{\infty }}{lim}f(t_0)=\underset{ϵ\to 0}{lim}ϵ{\int }_{-\mathrm{\infty }}^{0}d{t}^{\prime }{e}^{ϵt^{\prime }}f(t^{\prime }).$$

We may apply this prescription to the time evolution operator to give a mathematical meaning to $U_I(t,-\mathrm{\infty })$ or $U_I(+\mathrm{\infty },t)$. Thus

(A46) $$U_I(t,-\mathrm{\infty })=\underset{ϵ\to 0}{lim}ϵ{\displaystyle {\int }_{-\mathrm{\infty }}^0}d{t}^{\prime }{e}^{ϵt^{\prime }}{U}_I(t,t^{\prime }),$$

(A47) $$U_I(+\mathrm{\infty },t_0)=\underset{ϵ\to 0}{lim}ϵ{\displaystyle {\int }_0^{+\mathrm{\infty }}}d{t}^{\prime }{e}^{-ϵt^{\prime }}{U}_I(t^{\prime },t_0).$$

Appendix C

By defining $G={(E_i-H+iϵ)}^{-1}$ and $G_0={(E_i-H_0+iϵ)}^{-1}$ and noting that $G_0^{-1}-G^{-1}=H-H_0=V$, we take the difference $G-G_0$ as

(A48) $$G-G_0=(G_0{G}_0^{-1})G-G_0(G^{-1}G)$$ $$=G_0(G_0^{-1}-G^{-1})G=G_{0}VG,$$

which is the identity in Eq. (A13):

(A49) $${\displaystyle \frac{1}{E_i-H+iϵ}}-{\displaystyle \frac{1}{E_i-H_0+iϵ}}$$ $$={\displaystyle \frac{1}{E_i-H_0+iϵ}}V{\displaystyle \frac{1}{E_i-H+iϵ}}.$$

Alternatively, we can also write

(A50) $$G-G_0=G(G_0^{-1}{G}_0)-(G{G}^{-1})G_0$$ $$=G(G_0^{-1}-G^{-1})G_0=GV{G}_0,$$

which results in

(A51) $${\displaystyle \frac{1}{E_i-H+iϵ}}-{\displaystyle \frac{1}{E_i-H_0+iϵ}}$$ $$={\displaystyle \frac{1}{E_i-H+iϵ}}V{\displaystyle \frac{1}{E_i-H_0+iϵ}}.$$

References

• ^[1]
  D.J. Griffiths, Introduction to Quantum Mechanics (Cambridge University Press, Cambridge, 2018).
• ^[2]
  R.G. Newton, Scattering Theory of Waves and Particles (Dover Publications, Mineola, 2013).
• ^[3]
  J.J. Sakurai and J. Napolitano, Modern Quantum Mechanics (Oxford University Press, Oxford, 2014).
• ^[4]
  H.A. Bethe, Physical Review 76, 1 (1949).
• ^[5]
  L.S. Rodberg and R.M. Thaler, Introduction to the Quantum Theory of Scattering (Academic Press, Cambridge, 1967).
• ^[6]
  L.B. Madsen, American Journal of Physics 70, 8 (2002).
• ^[7]
  A.J. Moerdijk, B.J. Verhaar and A. Axelsson, Physical Review A 51, 4852 (1995).
• ^[8]
  M. Randeria and E. Taylor, Annual Review of Condensed Matter Physics 5, 209 (2014).
• ^[9]
  A. Gezerlis and J. Carlson, Physical Review C 77, 025803 (2008).
• ^[10]
  G.C. Strinati, P. Pieri, G.Röpke, P. Schuck and M. Urban, Physics Reports 738, 1 (2018).
• ^[11]
  L. Madeira and V.S. Bagnato, Brazilian Journal of Physics 51, 170 (2021).
• ^[12]
  M.A.C. Ribeiro, V.C Franzoni, W.R. Passos, E.C. Silva and A.N.F. Aleixo, Revista Brasileira de Ensino de Fisica 26, 1 (2004).
• ^[13]
  A.S. Castro, Revista Brasileira de Ensino de Fisica 33, 4312 (2011).
• ^[14]
  L. Rizzi, O.F. Piattella, S.L. Cacciatori and V. Gorini, Revista Brasileira de Ensino de Fisica 38, e2302 (2016).
• ^[15]
  E.O.A Landim and A.F.R. Rodrigues, Revista Brasileira de Ensino de Fisica 44, e20210338 (2022).
• ^[16]
  L.M. Cavalcanti, Revista Brasileira de Ensino de Fisica 21, 336 (1999).
• ^[17]
  J. Pera and J. Boronat, American Journal of Physics 91, 90 (2023).
• ^[18]
  L.B. Carvalho, W.C. Santos and E.A. Correa, Revista Brasileira de Ensino de Fisica 41, e20180359 (2019).
• ^[19]
  V.D. Viterbo, N.H.T. Lemes and J.P. Braga, Revista Brasileira de Ensino de Fisica 36, 1310 (2014).
• ^[20]
  F. Caruso and V. Oguri, Revista Brasileira de Ensino de Física 36, 2310 (2014).
• ^[21]
  G.B. Arfken and H.J. Weber, Mathematical Methods for Physicists (Elsevier Academic Press, Cambridge, 2005).
• ^[22]
  W. Cencek, M. Przybytek, J. Komasa, J.B. Mehl, B. Jeziorski and K. Szalewicz, The Journal of Chemical Physics 136, 22 (2012).
• ^[23]
  R.W. Hackenburg, Physical Review C 73, 044002 (2006).
• ^[24]
  M.J. Jamieson, A. Dalgarno and M. Kimura, Physical Review A 51, 2626 (1995).
• ^[25]
  A.R. Janzen and R.A. Aziz, The Journal of Chemical Physics 103, 22 (1995).
• ^[26]
  E. Tiesinga, P.J. Mohr, D.B. Newell and B.N. Taylor, CODATA Internationally recommended 2018 values of the Fundamental Physical Constants, available in: http://physics.nist.gov/constants
  » http://physics.nist.gov/constants
• ^[27]
  N. Giordano and H. Nakanishi, Computational Physics (Prentice Hall, London, 2005).
• ^[28]
  W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery, Numerical Recipes in Fortran 77: The Art of Scientific Computing (Cambridge University Press, Cambridge, 1996).
• ^[29]
  J. Carlson, S. Chang, V. Pandharipande and K. Schmidt, Physical Review Letters 91, 050401 (2003).
• ^[30]
  S. Függe, Practical Quantum Mechanics (Springer, Berlin, 1994).
• ^[31]
  L. Madeira, S. Gandolfi, K.E. Schmidt and V.S. Bagnato, Physical Review C 100, 014001 (2019).
• ^[32]
  M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions (Dover Publications, Mineola, 1965).
• ^[33]
  L.D. Landau and E.M. Lifschitz, Quantum Mechanics (Non-Relativistic Theory) (Pergamon Press, Oxford, 1977).
• ^[34]
  P. Jeszenszki, A.Y. Cherny and J. Brand, Physical Review A 97, 042708 (2018).
• ^[35]
  M. Gell-Mann and M.L. Goldberger, Physical Review 91, 398 (1953).
• ^[36]
  P. Roman, Advanced Quantum Theory (Addison-Wesley, Reading, 1965).

Publication Dates

• Publication in this collection
  21 July 2023
• Date of issue
  2023

History

• Received
  18 Mar 2023
• Accepted
  24 May 2023