\chapter{Resonances}
\label{chap:Resonances}

\section{Introduction}
\label{sec:res_intro}

\section{Complex-scaling method}
\label{sec:CSM}

For most physically relevant cases, we are work with finite-ranged potentials, which are rigorously constrained such that,
\begin{subalign}
    \lim_{r \to 0} r^2 V(r) &= 0 \, , \\
    \lim_{r \to \infty} r^3 V(r) &= 0 \, . \label{eq:finite_ranged}
\end{subalign}
It can be shown that the condition Eq.~\eqref{eq:finite_ranged} implies that the wave function asymptotically converges to Eq.~\eqref{eq:free_wf} as $r \to \infty$. Then we can write,
\begin{equation}
    u(r) \xrightarrow[r \rightarrow \infty]{} \frac{\ii}{2} \left[ \hat{h}^-_l(pr) - s_l(p) \hat{h}^+_l(pr) \right] \, .
    \label{eq:asymptotic_wf}
\end{equation}
Here, $s_l(p)$ can be interpreted as the partial-wave $S$ matrix, since it appears as the ratio between the outgoing and incoming waves. Equation~\eqref{eq:asymptotic_wf} implies that wherever $s_l(p)$ has a pole, we have
\begin{equation}
\label{eq:pobc}
 u(r) \xrightarrow[r \rightarrow \infty]{} \ii^l \gamma \, \hat{h}^+_l(pr) \,,
\end{equation}
where $\gamma$ is the asymptotic normalization constant and the $\ii^l$ factor is conventional. An illustration of the analytic structure of the $S$ matrix and its poles is shown in Fig.~\ref{fig:s_matrix}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
    \centering
    \includegraphics[width=0.65\textwidth]{Chapter-4/s_matrix.pdf}
    \caption{The analytic structure of the $S$ matrix indicating its poles in the $p$ plane corresponding to bound states, virtual states, resonances, and antiresonances (capturing resonances).
    Note especially how bound states lie on the positive imaginary axis while resonances are located in the \nth{4} quadrant.}
    \label{fig:s_matrix}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

As shown in in Fig.~\ref{fig:s_matrix}, bound states have associated imaginary momenta $p=\ii\kappa$ with real $\kappa>0$, whereas resonances are described by complex $p$ with $\Ip(p)<0$. This means that asymptotically, bound state wave functions decay exponentially with $r$ (hence the term ``bound'') whereas resonance wave functions grow exponentially and are therefore---like scattering states but in some sense even more so---not square-integrable; \ie, they do not correspond to normalizable states in the ordinary Hilbert space. While the rigged Hilbert space construction offers a rigorous mathematical formalism to deal with this difficulty (see, for example, \cite{delaMadrid:2012aa} for an introduction), for practical calculations there exists a much simpler alternative.

The so-called (uniform) complex-scaling method~\citep{Aguilar:1971ve,Balslev:1971vb, Moiseyev:1978aa,Moiseyev:1998aa,Afnan:1991kb,ho_method_1983,Reinhardt:1982aa} enables a description of resonances with, essentially, bound-state techniques. This is achieved by expressing the wave function not as usual along the real $r$ axis, but on a contour rotated into the complex-$r$ plane. This can be achieved by applying the transformation
%
\begin{equation}
 r \to r e^{\ii \phi}
\label{eq:r-scaled}
\end{equation}
%
to Eq.~\eqref{eq:SE_rad}, with some angle $\phi$. The proper choice of $\phi$ in general depends on the position of the resonance one wishes to study. If the state of interest has a complex energy $E$, then it is necessary to ensure that $\phi > {-}\arg(E)/2$. As $E$ is usually not known beforehand, one might repeat the calculation while increasing $\phi$ until a resonance is found.
\footnote{One might think of simply setting $\phi=\pi/4$ to accommodate all possible resonances. However, in most cases, large $\phi$ angles lead to potentials not vanishing fast enough along the contour, thereby demanding higher momentum cutoffs, or even potentials becoming completely divergent.}

With the convention in Eq.~\eqref{eq:r-scaled}, $r$ is still a real parameter but no longer equals to the physical radial coordinate of the system. The overall argument $kr e^{\ii \phi}$ of the Riccati-Hankel function in Eq.~\eqref{eq:pobc} satisfies $\Ip(kr e^{\ii \phi})>0$, and therefore square-integrability of the wave function as a function of $r$ is recovered. An example of such a scaled wave function is illustrated in Fig.~\ref{fig:csm_illustration}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}
    \centering
    \includegraphics[width=0.8\textwidth]{Chapter-4/csm_illustration.pdf}
    \caption{Effect on the reduced radial wave function of a typical $S$-wave resonance (Gamow state), due to complex-scaling of the $r$ contour. The solid (dotted) line corresponds to the real (imaginary) parts. After complex-scaling, it asymptotically converges to an exponentially decaying Riccati-Hankel function, $\hat{h}^+_0(pr e^{\ii\phi})=\hat{h}^+_0(\Tilde{p}r)=\exp(\ii \Tilde{p} r)$, where we define $\Tilde{p} = p e^{\ii \phi}$, the \emph{effective} wave number with $\Ip(\Tilde{p})>0$, so that it is normalizable just like bound-state wave functions.}
    \label{fig:csm_illustration}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

It was shown in Ref.~\cite{Afnan:1991kb} that the scaling of the radial coordinate $r$ is equivalent to a rotation in momentum representation that goes in the opposite (clockwise) direction with the same angle $\phi$. That is, if we consider the wave function of the state as a function of a momentum coordinate $q$, then complex scaling is implemented via
%
\begin{equation}
 q \to q e^{{-}\ii \phi} \,.
 \label{eq:q-scaled}
\end{equation}
%
This procedure then makes it possible to alternatively calculate resonance wave functions in momentum space. Note furthermore that scaling in momentum space can also be understood as a rotation of the branch cut in the complex-energy plane by an angle $2\phi$ clockwise, thereby exposing a section of the second Riemann sheet where resonances are located.

After this transformation, we can absorb the $e^{\ii \phi}$ phase into the wave number $k$ and define the \emph{effective} wave number as $\Tilde{k} = k e^{\ii \phi}$,
so that the asymptotic form in Eq.~\eqref{eq:pobc} is preserved as
\begin{equation}
 \psi_{l,k}(r e^{\ii \phi}) \xrightarrow[r \rightarrow \infty]{} N \, \hat{h}^+_l(\Tilde{k} r) \,.
\end{equation}
The scaling technique can thus be interpreted as mapping resonances from the \nth{4} quadrant in the complex-$k$ plane to the \nth{1} quadrant in the $\Tilde{k}$ plane.
For future reference, we note that, at the same time, it will effectively map bound states from the positive imaginary $k$ axis to the \nth{2} quadrant in the complex $\Tilde{k}$ plane.

\section{Complex scaling in 3D Cartesian coordinates}

While the method of complex scaling is easier to explain in a partial-wave framework~\citep{Afnan:1991kb,ho_method_1983}, the equivalent 3D Cartesian formulation, which is most appropriate for the finite-volume studies in Sec.\ref{TBA}, needs to be carefully stated. To that end we note that complex scaling of each individual component of $\vec{r}=(x,y,z)$ is equivalent to complex scaling of the radial coordinate $r$,
%
\begin{equation}
 r = \sqrt{x^2+y^2+z^2} \to \sqrt{(x e^{\ii \phi})^2+(y e^{\ii \phi})^2+(z e^{\ii \phi})^2} = r e^{\ii \phi} \,,
\label{eq:radial_coordinate_scaling}    
\end{equation}
%
but it leaves the angles $\theta$ and $\varphi$ in spherical coordinates $\vecr = (r,\theta,\phi)$ unaffected:
%
\begin{subalign}
 \cos{\theta} &= \frac{z}{r} = \frac{z e^{\ii \phi}}{r e^{\ii \phi}} \,, \\
 \tan{\varphi} &= \frac{y}{x} = \frac{y e^{\ii \phi}}{x e^{\ii \phi}} \,.
\end{subalign}
%
Therefore, we can apply complex scaling directly to Cartesian coordinates simply via $\vec{r}_i \to \vec{r}_i \ee^{\ii \phi}$. Note that we have implicitly altered usual Euclidean norm in Eq.~\eqref{eq:radial_coordinate_scaling} so it that preserves complex scaling, \viz
%
\begin{equation}
 |\vec{r}| = \sqrt{x^2+y^2+z^2}
 \mathtext{for} \vec{r}=(x,y,z) \mathtext{where} x,y,z \in \CC \,,
\end{equation}
instead of,
\begin{equation}
 |\vec{r}| \neq \sqrt{|x|^2+|y|^2+|z|^2} \, .
 \label{eq:norm-incorrect}
\end{equation}
This can be justified in the spirit of analytic continuation, because the standard Euclidean norm showing in Eq.~\eqref{eq:norm-incorrect} would not be analytic.

Finally, letting $\vec{q}$ be the momenta conjugate to $\vec{r}$, we see that the transformation
%
\begin{equation}
    \vec{q} \to \vec{q} \ee^{{-}\ii \phi}
\end{equation}
%
is necessary to preserve the canonical commutation relation,
%
\begin{equation}
    \left[ \vec{r} \ee^{\ii \phi},\vec{q} \ee^{{-}\ii \phi} \right] = \left[ \vec{r},\vec{q} \right] =  \ii \, .
\end{equation}
%
This agrees nicely with the claim made in Eq.~\eqref{eq:q-scaled}.

\subsection{Non-Hermiticity and the c-product}

In traditional quantum mechanics, one requires the Hamiltonian $H$ to be Hermitian ($H^\dagger=H$) to ensure that the energy spectrum, being a physical observable, is real and that time evolution is strictly unitary, \ie, the norm of quantum states are preserved under the time evolution operator $\ee^{{-}\ii H t / \hbar}$. However, when considering decay, an inherently time-dependent phenomenon, in a time-independent framework such as the complex-scaling method, the Hamiltonian is no longer Hermitian. Instead, in the present case, it becomes complex symmetric~\citep{Moiseyev:1998aa}, \ie,
\begin{equation}
H^\intercal=H \, .
\end{equation}
This permits the energy spectrum to include complex eigenvalues, which is precisely what is needed to describe resonances. In fact, the non-Hermiticity and the corresponding non-unitary time evolution of Gamow states are well aligned with the physical interpretation of resonances as metastable systems that ultimately decay.

Similarly to how non-degenerate eigenvectors of a Hermitian operator are orthogonal under the inner product defined on the Hilbert space, the non-degenerate eigenvectors of a complex symmetric operator are  orthogonal under the so-called ``c-product''~\cite{Moiseyev:1978aa,Moiseyev:2011}.
For eigenstates $\ket{\psi_1}$ and $\ket{\psi_2}$ with equal angular-momentum quantum numbers $(l,m)$, we define the c-product in coordinate representation as
%
\begin{equation}
 \braket{\psi_1 | \psi_2} = \int \dd r \, \psi_1(r)\psi_2(r) \,,
\label{eq:CP-vecr}
\end{equation}
%
and similarly in momentum space. Note that $\psi_1(r)$ appears without complex conjugation under the integral. This is precisely the c-product introduced in Ref.~\cite{Moiseyev:1978aa} with the notation $(\psi_1|\psi_2)$. In this manuscript, we use the standard notation $\braket{\psi_1 | \psi_2}$ with the implicit
understanding that for complex-scaled systems this is meant to denote the c-product.

Equivalently, one can change the definition of bra states so that no complex conjugation is involved when they are associated with a complex-scaled system. This is so even for bound states calculated with complex scaling. Although the energies of such states remain real, wave functions become complex when defined along the rotated contour and the orthogonality of states with different binding energies is ensured only if no complex conjugation is performed for bras, leading again to the c-product~\cite{Moiseyev:2011}. Ultimately, these concepts can be understood by properly distinguishing bra and ket states as, respectively, left and right
eigenvectors of the non-Hermitian complex-scaled Hamiltonian~\cite{Afnan:1991kb}. Even more rigorously, a comprehensive theory for Gamow bras and kets can be developed within the RHS formalism mentioned previously~\cite{delaMadrid:2012aa}. However, in practice we find it convenient and sufficient to employ complex scaling along with the c-product.

To prove c-orthogonality for complex symmetric $H$, we start by considering two non-degenerate eigenvectors (written as vectors instead of bras or kets for clarity),
%
\begin{subalign}
    H \vec{v_1} = E_1 \vec{v_1} \, ,
    \label{eq:eig_a}\\
    H \vec{v_2} = E_2 \vec{v_2} \, ,
    \label{eq:eig_b}        
\end{subalign}
%
where $E_1 \neq E_2$. Taking the transpose of \ref{eq:eig_b} gives,
%
\begin{equation}
    \label{eq:eig_trans}
    \vec{v_2}^\intercal H = E_2 \vec{v_2}^\intercal \, ,
\end{equation}
%
since $H^\intercal=H$. Note that complex conjugation was not involved. Now left-multiply \ref{eq:eig_a} by $v_2^\intercal$, right-multiply \ref{eq:eig_trans} by $v_1$, and finally subtract them to get,
%
\begin{spliteq}
    \vec{v_2}^\intercal H \vec{v_1} - \vec{v_2}^\intercal H \vec{v_1}
    & = (E_1-E_2) \vec{v_2}^\intercal \vec{v_1}\\
    0 & = (E_1-E_2) \vec{v_2}^\intercal \vec{v_1}\\
    0 & = \vec{v_2}^\intercal \vec{v_1} \, . \qed
\end{spliteq}

\ny{TODO: Self-orthogonal states}

\section{Purely outgoing boundary conditions}
\label{sec:POBC}

One can calculate Gamow states (and bound states) by solving the Schr\=odinger equation using an differential equation solver, subject to \emph{purely outgoing boundary conditions} (POBC).

The asymptotic behavior for the regular function is
\begin{equation}
    \phi_p(r) \xrightarrow[r \rightarrow \infty]{} \frac{\ii}{2} \left[ \Jostf_l(p) e^{-\ii pr} - \Jostf_l(-p) e^{\ii pr}  \right] \, .
\end{equation}
Then the first derivative is
\begin{equation}
    \phi'_p(r) \xrightarrow[r \rightarrow \infty]{} \frac{p}{2} \left[ \Jostf_l(p) e^{-\ii pr} + \Jostf_l(-p) e^{\ii pr}  \right] \, .
\end{equation}

Solving for $\Jostf_l(\pm p)$ gives
\begin{equation}
    \left( \frac{\phi'_p}{p} \pm \ii \phi_p \right) e^{\mp \ii pr} \xrightarrow[r \rightarrow \infty]{} \Jostf_l(\pm p) \, .
\end{equation}
But we know that $\Jostf_l(p) = 0$ for bound states and Gamow states. Therefore, for an integration interval $(0,R]$,
the boundary conditions
\begin{spliteq}
    \phi_p(0) &= 0 \\
    \frac{\phi'_p(R)}{p} + \ii \phi_p(R) &= 0 \\
\end{spliteq}
would yield the discrete bound state and Gamow state energies $E = \frac{p^2}{2\mu}$.

\section{Zeldovich regularization}
\label{sec:Zeldovich}

\ny{Only worked out for $S$-wave. Need to extend.}

Here we discuss Zeldovich regularization, as a method to calculate inner products
of configuration space wave functions determined via POBC as discussed in~\ref{sec:POBC}.
We will exploit the knowledge of asymptotic forms of wave functions to define a regularized inner-product.
Using this, we can implement CA-EC in configuration space, without the complications
of complex-scaled contours.

As $\Ip(p)<0$ for resonance states,
the wave function is ever increasing and not square integrable.
They are said to belong in a ``Rigged Hilbert Space''. Inner products and
normalization factors of Gamow states can be calculated by inserting a regularization
factor inside integrals which are otherwise divergent \cite{delaMadrid:2008sr},
allowing us to work directly in the RHS.

Let,
%
\begin{equation}
    \label{eq:zeldovich}
    \braket{\phi_1 | \phi_2}
    =\lim_{\mu\to 0}\int_{0}^{\infty}e^{-\mu r^2}\phi_1(r)\phi_2(r)\,dr
\end{equation}
%
Now assume that after some $r>R$, the wave functions becomes approximately equal to
their asymptotic forms, $N_1 e^{ip_1r}$ and $N_2 e^{ip_2r}$ respectively. Then the
integral can be broken into two.
%
\begin{spliteq}
    \braket{\phi_1 | \phi_2} 
    & \approx \lim_{\mu\to 0}\int_{0}^{R}e^{-\mu r^2}\phi_1(r)\phi_2(r)\,dr \\
    & \quad +N_1N_2\lim_{\mu\to 0}
    \int_{R}^{\infty}e^{-\mu r^2}e^{ip_1r}e^{ip_2r}\,dr \\
\end{spliteq}
%
where $N_{1,2}$ can be determined from $N=\phi(r)/e^{ipr}$ for any $r>R$. Using
results from \cite{kukulin1989theory} for the second term, and taking the limit in
the first term, this simplifies as follows.
%
\begin{equation}
    \label{eq:easy}
    \braket{\phi_1 | \phi_2} \approx \int_{0}^{R}\phi_1(r)\phi_2(r)\,dr
    +\frac{\ii N_1 N_2 e^{\ii(p_1+p_2)R}}{p_1+p_2}
\end{equation}
%
If $\phi_1(r)$ and $\phi_2(r)$ are solutions to finite-ranged potentials with ranges
$a_1$ and $a_2$, this result can be made exact by taking $R=\max\{a_1,a_2\}$.
For other potentials, a sufficiently large $R$ has to be considered such that
$V(r)\approx 0$ for all $r>R$.

Now, we have all the necessary tools for performing Galerkin projection on to the
reduced basis, to arrive at a smaller generalized eigenvalue problem. Rest of
the procedure would continue similarly.

\subsection{RHA-EC for Zeldovich method}

RHA-EC is straightforward in this configuration space and we can directly construct
free wave functions of the form $e^{\ii p r}$ and include them in the basis.

\subsection{CA-EC for Zeldovich method}

\ny{I need to check this math again. For now, please ignore this section.}

There will be no saving of additional memory usage by implementing CA-EC in
un-complex-scaled bases. Since RHA-EC already works better than CA-EC and is much
easier to implement, there is no reason to go with CA-EC,
other than for sanity checks, for which this section is dedicated.

\begin{lemma}
    For bound states, complex conjugation is equivalent to analytically rotating
    the function by an angle $2\phi$.
\end{lemma}

Consider a bound state training point $\ket{\psi}$.
%
\begin{equation}
    \ket{\psi}=\int_{0}^{\infty} dp\,p^2\, u(p)\ket{p}
\end{equation}
%
Now complex-scale the integration contour by an angle $-\phi$ (i.e., $\phi$ in the
clockwise direction).
%
\begin{equation}
    \ket{\psi}=\int_{-\phi} dp\,p^2\, u(p)\ket{p}
\end{equation}
%
This $u(p)$ is the reduced radial wave function that we obtain when we solve for
eigenstates with complex scaling. Now, consider the complex-conjugated wave function
on this contour.
%
\begin{equation}
    \ket{\Tilde{\psi}}=\int_{-\phi} dp\,p^2\, u^*(p)\ket{p}
\end{equation}
%
Since $\ket{\psi}$ is a bound state, $u(p)$ is real on the positive real axis.
Therefore, Schwarz reflection principle can be invoked to say $u^*(p)=u(p^*)$.
%
\begin{spliteq}
    \ket{\Tilde{\psi}}&=\int_{-\phi} dp\,p^2\, u(p^*)\ket{p}\\
    &=\int_{-\phi} dp\,p^2\, u(pe^{2i\phi})\ket{p}
\end{spliteq}
%
Since $u(pe^{2\ii\phi})$ is an analytic function, we can complex-scale the
contour back to the positive real axis.
%
\begin{equation}
    \ket{\Tilde{\psi}}=\int_{0}^{\infty} dp\,p^2\, u(pe^{2\ii\phi})\ket{p}
\end{equation}
%
The end result $\ket{\Tilde{\psi}}$ is simply the original state $\ket{\psi}$
rotated by an angle $2\phi$ on the complex $p$-plane.

Therefore, the equivalent of complex-conjugated vectors for the Zeldovich method
can be obtained by solving the Schr\"odinger's equation along a complex-scaled contour
($re^{-2\ii\phi}$) for some arbitrary angle $\phi$. As mentioned above, this method
is very inefficient in the Zeldovich implementation compared to RHA-EC.

\subsection{Analytic continuation of bound states}
\label{sec:analytic_cont}

In bound-to-resonance extrapolation, we need to work in a complex-scaled
basis to accommodate resonances.
We need the bound state wave functions expanded on this complex-scaled contour to use as
training data. These can be simply obtained by diagonalizing the Hamiltonian, analytically
continued onto the complex-scaled basis.

However, if one already has access to a set of pre-calculated bound state wave
functions on the real (un-rotated) contour, they
can be analytically continued onto the complex-scaled contour, without having to solve the Schr\"odinger's equation again, using the homogeneous equation,
%
\begin{equation}
    \ket{\psi}=\frac{1}{E-H_0}V\ket{\psi} \, ,
    \label{eq:homogenous}
\end{equation}
%
where $G(E)$ is the Green's operator defined as
%
\begin{equation}
    G(E)=\frac{1}{E-H_0} \, .
\end{equation}
%

Rewriting the above equation as in partial-wave notation gives
%
\begin{equation}
 u(\Tilde{q})=\int_{0}^{\infty} q\,dq\, \frac{1}{E-\frac{q^2}{2\mu}}V(\Tilde{q},q) u(q) \, .
\end{equation}
%
Now the reduced radial wave function $u(\Tilde{q})$ can be calculated on the complex-scaled contour in momentum space.

Equivalently, in configuration space, the Eq.~\eqref{eq:homogenous} reads,
%
\begin{equation}
 \psi(\Tilde{\vec{r}})={-}\frac{\mu}{2 \pi} \int_{0}^{\infty} d^3 \vec{r} \, \frac{e^{\ii p |\Tilde{\vec{r}}-\vec{r}|}}{|\Tilde{\vec{r}}-\vec{r}|} V(\vec{r}) \psi(\vec{r}) \, .
\end{equation}
%
However, as we work in a partial-wave basis, we should be able to find a simpler version of this equation that does not involve a 3-dimensional integration over $\vec{r}$.

To that end, let us start from Eq.~\eqref{eq:homogenous} and insert few complete basis sets of ``radial'' states introduced in Sec.~\ref{sec:partial_wave}.
\begin{spliteq}
    \braket{\Tilde{r},\Tilde{l},\Tilde{m} | \psi} &= \iint dr \, dr' \, \sum_{l,m} \sum_{l',m'} \bra{\Tilde{r},l,m} \frac{1}{E-H_0} \ket{r,l,m} \bra{r,l,m}V\ket{r',l',m'} \braket{r',l',m' | \psi} \\
    &= \int dr \, \sum_{l,m} \braketmatrix{\Tilde{r},\Tilde{l},\Tilde{m}}{\frac{1}{E-H_0}}{r,l,m} V(r) \braket{r,l,m | \psi}
\end{spliteq}

It is clear that we need an expression for $\braketmatrix{\Tilde{r},\Tilde{l},\Tilde{m}}{\frac{1}{E-H_0}}{r,l,m}$ which is the Green's function expanded in a partial-wave basis. This can be carried out as follows.
\begin{multline}
    \braketmatrix{\Tilde{r},\Tilde{l},\Tilde{m}}{\frac{1}{E-H_0}}{r,l,m} \\
     = \iint dE' \, dE'' \, \sum_{l',m'} \sum_{l'',m''} \braket{\Tilde{r},\Tilde{l},\Tilde{m} | E',l',m'} \braketmatrix{E',l',m'}{\frac{1}{E-H_0}}{E'',l'',m''} \braket{E'',l'',m'' | r,l,m} \\
     = \int dE' \, \sum_{l',m'} \braket{\Tilde{r},\Tilde{l},\Tilde{m} | E',l',m'} \frac{1}{E-E'} \braket{E',l',m' | r,l,m} \\
     = \delta_{\Tilde{l},l} \, \delta_{\Tilde{m},m} \int \frac{p'}{\mu} \, dp' \, \left[ \ii^l \sqrt{\frac{2\mu}{\pi p'}} \hat{j}_l(p' \Tilde{r}) \, \right]^* \frac{2\mu}{p^2-{p'}^2} \left[ \ii^l \sqrt{\frac{2\mu}{\pi p'}} \hat{j}_l(p' r) \, \right] \\
     = \delta_{\Tilde{l},l} \, \delta_{\Tilde{m},m} \frac{4\mu}{\pi} \int_0^\infty dp' \, \frac{\hat{j}_l(p' \Tilde{r}) \hat{j}_l(p' r)}{p^2-{p'}^2}
\end{multline}

This integral can be evaluated to obtain the partial-wave representation of the Green's function~\citep{Fuda:1973zz},
\begin{equation}
    \braketmatrix{\Tilde{r},\Tilde{l},\Tilde{m}}{\frac{1}{E-H_0}}{r,l,m} = {-} \frac{2 \mu}{p} \, \hat{j}_l(p r_{<}) \, \hat{h}^+_l(p r_{>}) \, \delta_{\Tilde{l},l} \, \delta_{\Tilde{m},m} \, ,
\end{equation}
where $r_<=\min(r,\Tilde{r})$ and $r_>=\max(r,\Tilde{r})$.
Finally, we have,
\begin{equation}
    u(\Tilde{r}) = {-}\frac{2 \mu}{p} \int_0^\infty dr \, \hat{j}_l(p r_{<}) \, \hat{h}^+_l(p r_{>}) V(r) u(r) \, .
\end{equation}