#  Complex eigenenergy of GDR for 16O by the Jost function within RPA
framework

K. Mizuyama1,2, N. Nhu Le3,4, T. V. Nhan Hao3,4, N. Hoang Tung5, and T. Dieu
Thuy3 corresponding author<EMAIL_ADDRESS>1 Institute of Research and
Development, Duy Tan University, Da Nang 550000, Vietnam
2 Faculty of Natural Sciences, Duy Tan University, Da Nang 550000, Vietnam
3 Faculty of Physics, University of Education, Hue University, 34 Le Loi
Street, Hue City, Vietnam
4 Center for Theoretical and Computational Physics, College of Education, Hue
University, 34 Le Loi Street, Hue City, Vietnam
5 Faculty of Physics, The University of Danang, University of Science and
Education, Da Nang 550000, Vietnam

###### Abstract

The Jost function method is extended within the framework of the RPA theory to
find poles on the complex energy plane that exhibit complex RPA eigenenergies.
Poles corresponding to the RPA excited states such as the giant resonance of
16O electric dipole excitations were successfully found on the complex energy
plane. Although the giant resonance has been known as a single resonance with
large strength and width, it is found that, at least withing the RPA
framework, the 16O electric dipole giant resonance is formed by multiple
poles, each of which is an independent pole with different widths, origins,
response properties to residual interactions, and components structures of the
density fluctuation..

##  I Introduction

In nuclear physics, the giant resonance is known as a resonance with a very
broad and large strength that appears in the relatively higher continuum
energy region above threshold in the cross section or strength function
expressed as a function of the excitation energy baldwin ; berman ; harakeh ;
bortignon . The random-phase-approximation (RPA) bohm  is a powerful tool for
the description of the collective excitation of nuclei, such as a giant
resonance liu ; ring-speth . Based on the understanding of the schematic model
of RPA ring , the RPA solution that is formed by the superposition of a large
number of $p$-$h$ excited configurations and gives a large strength due to the
effect of residual interactions is interpreted as a collective excitation
mode, and the giant resonance is considered to be such a collective excitation
mode.

There are two main methods for solving the equation in RPA theory: one that
uses a discrete basis and diagonalises the RPA Hamiltonian to obtain a
solution, and the continuum RPA (cRPA) shlomo  that takes into account the
boundary conditions of the continuum. In the former method, all RPA solutions
are obtained as discrete energy eigenstates, with giant resonances appearing
either as solutions giving large strength at the resonance energy or as many
spread discrete solutions near the resonance energy. The width is often
evaluated from the empirical reproducibility of the experimental data using
the Lorentz distribution function, assuming that the giant resonance has a
peak structure with a single width. Since the latter method (the cRPA method)
takes into account the boundary conditions of the continuum, it can represent
the behavior of the strength function as a smooth function of energy above
threshold. Without the assumption of a Lorentz distribution, as a result of
numerical calculations, the giant resonance is also represented as a smooth
function of energy, which has a peak structure with a width. However, the cRPA
method is not able to calculate the width itself.

There is a phenomenological model that describes the width of the giant
resonance in terms of mass number and dissipative term. It is well known that
the mass number dependence of experimental values of the giant resonance width
estimated by the Lorenz distribution function roughly follows this
phenomenological model speth ; aue . This phenomenological model is derived
from the Euler equations plus a dissipation term representing viscosity,
assuming that the giant resonance is a harmonic oscillation with damping due
to viscosity. In calculations such as the second RPA, which include higher-
order effects, the discrete strength (of resonances such as the giant
resonance) obtained from calculations using a discrete basis is fragmented and
spread in energy by the higher-order effects. This spreading of the strength
is interpreted as the broadening of the width of the resonance due to higher-
order effects (spreading width), however, this is not an argument made by
directly calculating the width. bortignon . Despite the fact that width is one
of the most important characters in the giant resonance, there is still no
method that can calculate the width of resonances obtained as the complex
eigenenergy states of the nucleus within the microscopic theoretical framework
such as the RPA theory.

The Jost function is a function that gives the energy eigenvalues of the
fundamental differential equations of a quantum system based on the
Hamiltonian of a system such as the Schrödinger equation as zeros on the
complex energy plane. The Jost function is given in differential and integral
forms, and its differential form is equivalent to the Wronskian. The energy
eigenvalues given by the zeros on the complex energy plane give not only bound
states but also resonance states. The imaginary part of the zeros on the
complex energy plane of the Jost function gives the width of the resonance.
However, the original Jost function jost  did not take into account channel
coupling (Note that the term “channel coupling” here includes a wide range of
meanings, such as the coupling between $p$-$h$ configurations). In order to
apply the Jost function to solve various physics problems, it is necessary to
extend the Jost function to be able to consider coupled channels. As a first
step, we have extended the Jost function to the Hartree-Fock-Bogoliubov theory
framework in Refs.jost-hfb ; jost-fano ; jost-class .

In this paper, we extend the Jost function method within the framework of RPA
theory in order to enable the Jost function method to find the complex
eigenenergies of the RPA solutions on the complex energy plane. The electric
dipole excitation of 16O is then chosen as the first application of the Jost-
RPA method and analyzed by calculating the poles of the RPA strength function
on the complex energy plane, adopting the Woods-Saxon potential for the mean
field, and the simple density dependent contact interaction for the residual
interaction.

##  II Formalism of Jost function for RPA

Since the Jost function is defined as a coefficient function relating regular
and irregular solutions of second-order differential equations such as the
Schrödinger equation, it is necessary to represent the RPA equation in the
form of a second-order differential equation in order to define and calculate
the Jost function within the framework of RPA theory.

In this section, we present the derivation of the Jost function within the
framework of RPA theory and the perturbed Green’s function, RPA response
function, and strength function using the Jost function.

###  II.1 Derivation of the Jost function

Defining the functions $X_{h}(\mbox{\boldmath$r$})$ and
$Y_{h}(\mbox{\boldmath$r$})$ as

$\displaystyle X_{h}(\mbox{\boldmath$r$};\omega)$ $\displaystyle=$
$\displaystyle\sum_{p}X_{ph}(\omega)\varphi_{p}(\mbox{\boldmath$r$})$ (1)
$\displaystyle Y_{h}(\mbox{\boldmath$r$};\omega)$ $\displaystyle=$
$\displaystyle\sum_{p}Y_{ph}(\omega)\varphi_{p}^{*}(\mbox{\boldmath$r$})$ (2)

using the $X$ and $Y$ amplitudes in the ordinary RPA equation known to be
expressed in the form $\begin{pmatrix}A&B\\\
B&A\end{pmatrix}\begin{pmatrix}X\\\
Y\end{pmatrix}=\hbar\omega\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}X\\\ Y\end{pmatrix}$ and the single particle
wave function $\varphi_{p}$, the differential equation for the partial wave
component of $X_{h}(\mbox{\boldmath$r$})$ and $Y_{h}(\mbox{\boldmath$r$})$ in
a spherically symmetric system is given by

$\displaystyle\begin{pmatrix}h_{lj}^{(q)}-\epsilon_{h}^{(q)}-\omega&0\\\
0&h_{lj}^{(q)}-\epsilon_{h}^{(q)}+\omega\end{pmatrix}\begin{pmatrix}X_{lj;h}^{L(q)}\\\
Y_{lj;h}^{L(q)}\end{pmatrix}$
$\displaystyle+\sum_{q^{\prime}h^{\prime}l^{\prime}j^{\prime}}\frac{\kappa_{qq^{\prime}}(r)}{r^{2}}\begin{pmatrix}\tilde{\varphi}_{lj;h}^{L(q)}\tilde{\varphi}_{l^{\prime}j^{\prime};h^{\prime}}^{L(q^{\prime})}&\tilde{\varphi}_{lj;h}^{L(q)}\tilde{\varphi}_{l^{\prime}j^{\prime};h^{\prime}}^{L(q^{\prime})}\\\
\tilde{\varphi}_{lj;h}^{L(q)}\tilde{\varphi}_{l^{\prime}j^{\prime};h^{\prime}}^{L(q^{\prime})}&\tilde{\varphi}_{lj;h}^{L(q)}\tilde{\varphi}_{l^{\prime}j^{\prime};h^{\prime}}^{L(q^{\prime})}\end{pmatrix}\begin{pmatrix}X_{l^{\prime}j^{\prime};h^{\prime}}^{L(q^{\prime})}\\\
Y_{l^{\prime}j^{\prime};h^{\prime}}^{L(q^{\prime})}\end{pmatrix}$
$\displaystyle=0$ (3)

where $h_{lj}^{(q)}$ is the mean field Hamiltonian given by

$\displaystyle h_{lj}^{(q)}=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}+U_{lj}^{(q)}(r)$ (4)

(the centrifugal potential $\frac{\hbar^{2}l(l+1)}{2mr^{2}}$ is included in
$U_{lj}^{(q)}(r)$), and $\kappa_{qq^{\prime}}(r)$ is the residual interaction,
and $\tilde{\varphi}_{lj;h}^{L(q)}$ is the function which is defined by

$\displaystyle\tilde{\varphi}_{lj;h}^{L}(r)\equiv\frac{\langle
lj||Y_{L}||l_{h}j_{h}\rangle}{\sqrt{2L+1}}\varphi_{h}^{(q)}(r)$ (5)

with use of the hole state wave function $\varphi_{h}^{(q)}(r)$ satisfies
$h_{l_{h}j_{h}}^{(q)}\varphi_{h}^{(q)}=e_{h}^{(q)}\varphi_{h}^{(q)}$. The
subscription $h$ describes a hole state quantum numbers as
$h\in(n_{h},l_{h},j_{h})$, and $q$ denotes neutron or proton by $q=n$ or $p$.

By introducing a subscription $\alpha$ which expresses the particle-hole
transition configuration (as shown in Table.2), and the momentum
$k_{1,\alpha}^{(q)}$ and $k_{2,\alpha}^{(q)}$ which are defined by

$\displaystyle k_{1,\alpha}^{(q)}(\omega)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2m}{\hbar^{2}}(e^{(q)}_{\alpha}+\omega)}$ (6)
$\displaystyle k_{2,\alpha}^{(q)}(\omega)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2m}{\hbar^{2}}(e^{(q)}_{\alpha}-\omega)},$ (7)

Eq.(3) can be rewritten as

$\displaystyle\sum_{q^{\prime}\alpha^{\prime}}\left[\frac{\hbar^{2}}{2m}\begin{pmatrix}k_{1,\alpha}^{(q)2}&0\\\
0&k_{2,\alpha}^{(q)2}\end{pmatrix}\delta_{qq^{\prime}}\delta_{\alpha\alpha^{\prime}}\right.$
$\displaystyle\left.-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\delta_{qq^{\prime}}\delta_{\alpha\alpha^{\prime}}\right.\right.$
$\displaystyle\left.\left.+U_{\alpha}^{(q)}\begin{pmatrix}1&0\\\
0&1\end{pmatrix}\delta_{qq^{\prime}}\delta_{\alpha\alpha^{\prime}}\right.\right.$
$\displaystyle\left.\left.+\frac{\kappa_{qq^{\prime}}(r)}{r^{2}}\begin{pmatrix}\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}&\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}\\\
\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}&\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}\end{pmatrix}\right\\}\right]\begin{pmatrix}X_{\alpha^{\prime}}^{(q^{\prime})}\\\
Y_{\alpha^{\prime}}^{(q^{\prime})}\end{pmatrix}$ $\displaystyle=0.$ (8)

Since $\alpha$ is defined for each multipolarity $L$ (i.e., the multipolarity
$L$ is fixed when $\alpha$ is defined), thereafter we will not explicitly show
$L$ in the formula.

Furthermore, Eq.(8) can be represented in the matrix form as

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}+\mbox{\boldmath$\mathcal{U}$}+\mbox{\boldmath$\mathcal{V}$}\right\\}\right]\vec{\mbox{\boldmath$\phi$}}=0$
(9)

where $\mathcal{K}$ and $\mathcal{U}$ are defined as the $N_{n}+N_{p}$
dimensional diagonal matrix defined as

$\displaystyle\mbox{\boldmath$\mathcal{K}$}=\begin{pmatrix}\begin{matrix}\mathcal{K}_{1}^{(n)}&0&\cdots&0\\\
0&\mathcal{K}_{2}^{(n)}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathcal{K}_{N_{n}}^{(n)}\end{matrix}&\begin{matrix}&&&\\\
&\text{\huge{0}}&&\\\ &&&\\\ &&&\\\ \end{matrix}\\\ \begin{matrix}&&&\\\
&\text{\huge{0}}&&\\\ &&&\\\ &&&\\\
\end{matrix}&\begin{matrix}\mathcal{K}_{1}^{(p)}&0&\cdots&0\\\
0&\mathcal{K}_{2}^{(p)}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathcal{K}_{N_{p}}^{(p)}\end{matrix}\end{pmatrix}$ (10)

and

$\displaystyle\mbox{\boldmath$\mathcal{U}$}=\begin{pmatrix}\begin{matrix}\mathcal{U}_{1}^{(n)}&0&\cdots&0\\\
0&\mathcal{U}_{2}^{(n)}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathcal{U}_{N_{n}}^{(n)}\end{matrix}&\begin{matrix}&&&\\\
&\text{\huge{0}}&&\\\ &&&\\\ &&&\\\ \end{matrix}\\\ \begin{matrix}&&&\\\
&\text{\huge{0}}&&\\\ &&&\\\ &&&\\\
\end{matrix}&\begin{matrix}\mathcal{U}_{1}^{(p)}&0&\cdots&0\\\
0&\mathcal{U}_{2}^{(p)}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathcal{U}_{N_{p}}^{(p)}\end{matrix}\end{pmatrix}$ (11)

by using the $2\times 2$ matrices

$\displaystyle\mathcal{K}_{\alpha}^{(q)}=\begin{pmatrix}k_{1,\alpha}^{(q)}&0\\\
0&k_{2,\alpha}^{(q)}\end{pmatrix}$ (12)

and

$\displaystyle\mathcal{U}_{\alpha}^{(q)}=U_{\alpha}^{(q)}\begin{pmatrix}1&0\\\
0&1\end{pmatrix},$ (13)

therefore $\mathcal{K}$ and $\mathcal{U}$ are totally $2(N_{n}+N_{p})\times
2(N_{n}+N_{p})$ dimensional diagonal matrix. $N_{q}$(for $q=n$ and $p$) is the
number of the p-h configurations. (For example, $N_{n}=N_{p}=7$ in the case of
the electric dipole excitation of 16O as shown in Table.2).

$\mathcal{V}$ is the $2(N_{n}+N_{p})\times 2(N_{n}+N_{p})$ matrix for the
residual interaction which is represented as

$\displaystyle\mbox{\boldmath$\mathcal{V}$}=$
$\displaystyle\begin{pmatrix}\begin{matrix}\mathcal{V}_{11}^{(nn)}&\mathcal{V}_{12}^{(nn)}&\cdots&\mathcal{V}_{1N_{n}}^{(nn)}\\\
\mathcal{V}_{21}^{(nn)}&\mathcal{V}_{22}^{(nn)}&\cdots&\mathcal{V}_{2N_{n}}^{(nn)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathcal{V}_{N_{n}1}^{(nn)}&\mathcal{V}_{N_{n}2}^{(nn)}&\cdots&\mathcal{V}_{N_{n}N_{n}}^{(nn)}\end{matrix}&\begin{matrix}\mathcal{V}_{11}^{(np)}&\mathcal{V}_{12}^{(np)}&\cdots&\mathcal{V}_{1N_{p}}^{(np)}\\\
\mathcal{V}_{21}^{(np)}&\mathcal{V}_{22}^{(np)}&\cdots&\mathcal{V}_{2N_{p}}^{(np)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathcal{V}_{N_{n}1}^{(np)}&\mathcal{V}_{N_{n}2}^{(np)}&\cdots&\mathcal{V}_{N_{n}N_{p}}^{(np)}\end{matrix}\\\
\begin{matrix}\mathcal{V}_{11}^{(pn)}&\mathcal{V}_{12}^{(pn)}&\cdots&\mathcal{V}_{1N_{n}}^{(pn)}\\\
\mathcal{V}_{21}^{(pn)}&\mathcal{V}_{22}^{(pn)}&\cdots&\mathcal{V}_{2N_{n}}^{(pn)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathcal{V}_{N_{p}1}^{(pn)}&\mathcal{V}_{N_{p}2}^{(pn)}&\cdots&\mathcal{V}_{N_{p}N_{n}}^{(pn)}\end{matrix}&\begin{matrix}\mathcal{V}_{11}^{(pp)}&\mathcal{V}_{12}^{(pp)}&\cdots&\mathcal{V}_{1N_{p}}^{(pp)}\\\
\mathcal{V}_{21}^{(pp)}&\mathcal{V}_{22}^{(pp)}&\cdots&\mathcal{V}_{2N_{p}}^{(pp)}\\\
\vdots&\vdots&\ddots&\vdots\\\
\mathcal{V}_{N_{p}1}^{(pp)}&\mathcal{V}_{N_{p}2}^{(pp)}&\cdots&\mathcal{V}_{N_{p}N_{p}}^{(pp)}\end{matrix}\end{pmatrix}$
(14)

with

$\displaystyle\mathcal{V}_{\alpha\alpha^{\prime}}^{(qq^{\prime})}=\frac{\kappa_{qq^{\prime}}(r)}{r^{2}}\begin{pmatrix}\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}&\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}\\\
\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}&\tilde{\varphi}_{\alpha}^{(q)}\tilde{\varphi}_{\alpha^{\prime}}^{(q^{\prime})}.\end{pmatrix}$
(15)

Note that $\mathcal{V}$ is a symmetric matrix. If we define
$2(N_{n}+N_{p})$-dimensional vector for the hole state wave functions
$\tilde{\varphi}_{\alpha}^{(q)}$ as

$\displaystyle\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{n}=\begin{pmatrix}\begin{pmatrix}\tilde{\varphi}_{1}^{(n)}\\\
\tilde{\varphi}_{1}^{(n)}\end{pmatrix}\\\
\begin{pmatrix}\tilde{\varphi}_{2}^{(n)}\\\
\tilde{\varphi}_{2}^{(n)}\end{pmatrix}\\\ \vdots\\\
\begin{pmatrix}\tilde{\varphi}_{N_{n}}^{(n)}\\\
\tilde{\varphi}_{N_{n}}^{(n)}\end{pmatrix}\\\ \begin{pmatrix}0\\\
0\end{pmatrix}\\\ \begin{pmatrix}0\\\ 0\end{pmatrix}\\\ \vdots\\\
\begin{pmatrix}0\\\ 0\end{pmatrix}\end{pmatrix}\mbox{, and
}\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{p}=\begin{pmatrix}\begin{pmatrix}0\\\
0\end{pmatrix}\\\ \begin{pmatrix}0\\\ 0\end{pmatrix}\\\ \vdots\\\
\begin{pmatrix}0\\\ 0\end{pmatrix}\\\
\begin{pmatrix}\tilde{\varphi}_{1}^{(p)}\\\
\tilde{\varphi}_{1}^{(p)}\end{pmatrix}\\\
\begin{pmatrix}\tilde{\varphi}_{2}^{(p)}\\\
\tilde{\varphi}_{2}^{(p)}\end{pmatrix}\\\ \vdots\\\
\begin{pmatrix}\tilde{\varphi}_{N_{p}}^{(p)}\\\
\tilde{\varphi}_{N_{p}}^{(p)}\end{pmatrix}\end{pmatrix},$ (16)

Eq.(14) can be represented as

$\displaystyle\mbox{\boldmath$\mathcal{V}$}=\sum_{q,q^{\prime}}\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{q}\frac{\kappa_{qq^{\prime}}(r)}{r^{2}}\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{q^{\prime}}^{\mathsf{T}}.$
(17)

$\vec{\mbox{\boldmath$\phi$}}$ is the RPA wave function which is defined by

$\displaystyle\vec{\mbox{\boldmath$\phi$}}=\begin{pmatrix}\phi_{1}^{(n)}\\\
\phi_{2}^{(n)}\\\ \vdots\\\ \phi_{N_{n}}^{(n)}\\\ \phi_{1}^{(p)}\\\
\phi_{2}^{(p)}\\\ \vdots\\\ \phi_{N_{p}}^{(p)}\end{pmatrix}$ (18)

as a $2(N_{n}+N_{p})$-dimensional vector using

$\displaystyle\phi_{\alpha}^{(q)}=\begin{pmatrix}X_{\alpha}^{(q)}\\\
Y_{\alpha}^{(q)}\end{pmatrix}.$ (19)

Eq.(9) is an $2(N_{n}+N_{p})$-dimensional simultaneous second-order
differential equation, which can be easily solved numerically by providing
appropriate boundary conditions for the given energy using the Numerov or
Runge-Kutta methods (the Numerov method is used in this paper). Since the
$2(N_{n}+N_{p})$-dimensional simultaneous second-order differential equations
have $2(N_{n}+N_{p})$ types of regular and non-regular solutions, the boundary
conditions for each are given as follows.

The regular solutions $\vec{\mbox{\boldmath$\phi$}}^{(r1;q\alpha)}$ (for
$\alpha\in(1,\cdots N_{q})$ and $q=n$ and $p$) are given as the solution
satisfying the boundary conditions at $r=0$ given by

$\displaystyle\lim_{r\to
0}\vec{\mbox{\boldmath$\phi$}}^{(r1;q\alpha)}=\begin{pmatrix}0\\\ \vdots\\\
0\\\
\phi_{\alpha}^{(q)}\to\begin{pmatrix}rj_{l_{\alpha}}(k_{1,\alpha}^{(q)}r)\\\
0\end{pmatrix}\\\ 0\\\ \vdots\\\ 0\end{pmatrix}$ (20)

The boundary conditions for the regular solutions
$\vec{\mbox{\boldmath$\phi$}}^{(r2;q\alpha)}$ (for $\alpha\in(1,\cdots N_{q})$
and $q=n$ and $p$) are given by

$\displaystyle\lim_{r\to
0}\vec{\mbox{\boldmath$\phi$}}^{(r2;q\alpha)}=\begin{pmatrix}0\\\ \vdots\\\
0\\\ \phi_{\alpha}^{(q)}\to\begin{pmatrix}0\\\
rj_{l_{\alpha}}(k_{2,\alpha}^{(q)}r)\end{pmatrix}\\\ 0\\\ \vdots\\\
0\end{pmatrix},$ (21)

$\vec{\mbox{\boldmath$\phi$}}^{(r2;q\alpha)}$ corresponds to the negative
energy solution of $\vec{\mbox{\boldmath$\phi$}}^{(r1;q\alpha)}$ because there
is a relation between $k_{1,\alpha}^{(q)}$ and $k_{2,\alpha}^{(q)}$ as
$k_{2,\alpha}^{(q)}(\omega)=k_{1,\alpha}^{(q)}(-\omega)$.

The outgoing boundary conditions at the limit $r\to\infty$ for the irregular
solutions $\vec{\mbox{\boldmath$\phi$}}^{(\pm 1;q\alpha)}$ and
$\vec{\mbox{\boldmath$\phi$}}^{(\pm 2;q\alpha)}$ are given by

$\displaystyle\lim_{r\to\infty}\vec{\mbox{\boldmath$\phi$}}^{(\pm
1;q\alpha)}=\begin{pmatrix}0\\\ \vdots\\\ 0\\\
\phi_{\alpha}^{(q)}\to\begin{pmatrix}rh^{(\pm)}_{l_{\alpha}}(k_{1,\alpha}^{(q)}r)\\\
0\end{pmatrix}\\\ 0\\\ \vdots\\\ 0\end{pmatrix},$ (22)

and

$\displaystyle\lim_{r\to\infty}\vec{\mbox{\boldmath$\phi$}}^{(\pm
2;q\alpha)}=\begin{pmatrix}0\\\ \vdots\\\ 0\\\
\phi_{\alpha}^{(q)}\to\begin{pmatrix}0\\\
rh^{(\pm)}_{l_{\alpha}}(k_{2,\alpha}^{(q)}r)\end{pmatrix}\\\ 0\\\ \vdots\\\
0\end{pmatrix},$ (23)

respectively, where $h_{l}^{(\pm)}$ is the spherical Hankel function defined
by $h_{l}^{(\pm)}(kr)=j_{l}(kr)\pm in_{l}(kr)$.

We can define the regular and irregular solution matrix
$\mbox{\boldmath$\Phi$}^{(r)}$ and $\mbox{\boldmath$\Phi$}^{(\pm)}$ using the
regular and irregular solutions as

$\displaystyle\mbox{\boldmath$\Phi$}^{(r)}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\vec{\mbox{\boldmath$\phi$}}^{(r1;n1)},&\vec{\mbox{\boldmath$\phi$}}^{(r2;n1)},&\cdots&\vec{\mbox{\boldmath$\phi$}}^{(r1;pN_{p})},&\vec{\mbox{\boldmath$\phi$}}^{(r2;pN_{p})}\end{pmatrix}$

and

$\displaystyle\mbox{\boldmath$\Phi$}^{(\pm)}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\vec{\mbox{\boldmath$\phi$}}^{(\pm
1;n1)},&\vec{\mbox{\boldmath$\phi$}}^{(\pm
2;n1)},&\cdots&\vec{\mbox{\boldmath$\phi$}}^{(\pm
1;pN_{p})},&\vec{\mbox{\boldmath$\phi$}}^{(\pm 2;pN_{p})}\end{pmatrix}.$

These matrices are given as the $2(N_{n}+N_{p})\times 2(N_{n}+N_{p})$ matrix
which satisfy

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}+\mbox{\boldmath$\mathcal{U}$}+\mbox{\boldmath$\mathcal{V}$}\right\\}\right]\mbox{\boldmath$\Phi$}^{(r)}$
$\displaystyle=$ $0$ (26)

and

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}+\mbox{\boldmath$\mathcal{U}$}+\mbox{\boldmath$\mathcal{V}$}\right\\}\right]\mbox{\boldmath$\Phi$}^{(\pm)}$
$\displaystyle=$ $0$ (27)

respectively.

The Jost function $\mbox{\boldmath$\mathcal{J}$}^{(\pm)}$ is defined as the
coefficient matrix which connects $\mbox{\boldmath$\Phi$}^{(r)}$ and
$\mbox{\boldmath$\Phi$}^{(\pm)}$ as

$\displaystyle\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\mbox{\boldmath$\mathcal{J}$}^{(+)}\mbox{\boldmath$\Phi$}^{(-)\mathsf{T}}+\mbox{\boldmath$\mathcal{J}$}^{(-)}\mbox{\boldmath$\Phi$}^{(+)\mathsf{T}}\right].$
(28)

The Wronskian $\mbox{\boldmath$\mathcal{W}$}^{(\pm)}$ is defined by

$\displaystyle\mbox{\boldmath$\mathcal{W}$}^{(\pm)}$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}}{2m}\left[\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(\pm)}\right)-\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}\right)\mbox{\boldmath$\Phi$}^{(\pm)}\right],$

and it is very easy to confirm that this definition of Wronskian is constant
for the radial coordinate $r$, i.e. $\frac{\partial}{\partial
r}\mbox{\boldmath$\mathcal{W}$}^{(\pm)}=0$.

By inserting Eq.(28) into Eq.(LABEL:wrondef) and taking the limit
$r\to\infty$, it is very easy to obtain the relation formula between the Jost
function and Wronskian given by

$\displaystyle\mbox{\boldmath$\mathcal{J}$}^{(\pm)}$ $\displaystyle=$
$\displaystyle\pm\frac{2m}{i\hbar^{2}}\mbox{\boldmath$\mathcal{W}$}^{(\pm)}\mbox{\boldmath$\mathcal{K}$}.$
(30)

Applying the Green’s theorem to Eq.(26), we can obtain

$\displaystyle\left[\mbox{\boldmath$\chi$}^{(\pm)\mathsf{T}}(r)\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(r)}(r)-\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\chi$}^{(\pm)\mathsf{T}}(r)\right)\mbox{\boldmath$\Phi$}^{(r)}(r)\right]\pm
i\mbox{\boldmath$\mathcal{K}$}^{-1}$
$\displaystyle=\frac{2m}{\hbar^{2}}\int_{0}^{r}dr^{\prime}\mbox{\boldmath$\chi$}^{(\pm)\mathsf{T}}(r^{\prime})\left[\mbox{\boldmath$\mathcal{U}$}(r^{\prime})+\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\right]\mbox{\boldmath$\Phi$}^{(r)}(r^{\prime})$
(31)

where $\mbox{\boldmath$\chi$}^{(\pm)}$ is the free particle wave function
matrix which satisfies

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}+\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}\right]\mbox{\boldmath$\chi$}^{(\pm)}$
$\displaystyle=$ $0$ (32)

and the components are represented by the spherical Hankel function.

By inserting Eq.(28) into Eq.(31) and taking the limit $r\to\infty$, we can
obtain the integral form of the Jost function

$\displaystyle\mbox{\boldmath$\mathcal{J}$}^{(\pm)\mathsf{T}}=\mbox{\boldmath$1$}\mp\frac{2m}{\hbar^{2}}\frac{1}{i}\mbox{\boldmath$\mathcal{K}$}\int_{0}^{\infty}dr^{\prime}\mbox{\boldmath$\chi$}^{(\pm)\mathsf{T}}(r^{\prime})$
$\displaystyle\hskip
40.0pt\times\left[\mbox{\boldmath$\mathcal{U}$}(r^{\prime})+\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\right]\mbox{\boldmath$\Phi$}^{(r)}(r^{\prime}).$
(33)

###  II.2 Perturbed Green function and RPA response function

The Green’s function is defined as a function which satisfies

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}+\mbox{\boldmath$\mathcal{U}$}+\mbox{\boldmath$\mathcal{V}$}\right\\}\right]\mbox{\boldmath$\mathcal{G}$}^{(\pm)}(r,r^{\prime})$
$\displaystyle=\mbox{\boldmath$1$}\delta(r-r^{\prime}).$ (34)

can be represented as

$\displaystyle\mbox{\boldmath$\mathcal{G}$}^{(\pm)}(r,r^{\prime})$
$\displaystyle=\theta(r-r^{\prime})\mbox{\boldmath$\Phi$}^{(\pm)}(r)\left(\mbox{\boldmath$\mathcal{W}$}^{(\pm)-1}\right)\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}(r^{\prime})$
(35)
$\displaystyle+\theta(r^{\prime}-r)\mbox{\boldmath$\Phi$}^{(r)}(r)\left(\mbox{\boldmath$\mathcal{W}$}^{(\pm)-1}\right)^{\mathsf{T}}\mbox{\boldmath$\Phi$}^{(\pm)\mathsf{T}}(r^{\prime})$

by using the Wronskian. Note that this Green function is also given in the
form of $2(N_{n}+N_{p})\times 2(N_{n}+N_{p})$ matrix (The proof of Eq.(35) is
shown in Appendix..1).

The RPA equation when the external field exists is given by

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}+\mbox{\boldmath$\mathcal{U}$}+\mbox{\boldmath$\mathcal{V}$}\right\\}\right]\vec{\mbox{\boldmath$\phi$}}_{F}=\vec{\mbox{\boldmath$F$}}$
(36)

where $\vec{\mbox{\boldmath$F$}}$ is the external field which is expressed as

$\displaystyle\vec{\mbox{\boldmath$F$}}(r)$ $\displaystyle=$
$\displaystyle\sum_{q=n,p}\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{q}f_{q}(r)$
(37)

in the form of the $2(N_{n}+N_{p})$-dimensional vector. The solution of the
RPA equation with the external field Eq.(36)
$\vec{\mbox{\boldmath$\phi$}}_{F}$ is given by using the Green function
$\mbox{\boldmath$\mathcal{G}$}^{(\pm)}(r,r^{\prime})$ as

$\displaystyle\vec{\mbox{\boldmath$\phi$}}_{F}(r)=\int_{0}^{\infty}dr^{\prime}\mbox{\boldmath$\mathcal{G}$}^{(+)}(r,r^{\prime})\cdot\vec{\mbox{\boldmath$F$}}(r^{\prime})$
(38)

and the strength function is given by

$\displaystyle S_{F}(\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{\pi}\mbox{ Im }\int
dr\vec{\mbox{\boldmath$F$}}^{\mathsf{T}}(r)\cdot\vec{\mbox{\boldmath$\phi$}}_{F}(r)$
$\displaystyle=$ $\displaystyle-\frac{1}{\pi}\mbox{ Im }\int\int
drdr^{\prime}\vec{\mbox{\boldmath$F$}}^{\mathsf{T}}(r)\cdot\mbox{\boldmath$\mathcal{G}$}^{(+)}(r,r^{\prime})\cdot\vec{\mbox{\boldmath$F$}}(r^{\prime})$

Since the inverse of the Wronskian is included in the Green function as shown
in Eq.(35) and the Wronskian is related to the Jost function as shown by
Eq.(30), the pole of the strength function may be found on the complex energy
$\omega$ plane as a solution of

$\displaystyle\det\mbox{\boldmath$\mathcal{J}$}^{(+)}(\omega)=0.$ (41)

The strength function can also be expressed as

$\displaystyle S_{F}(\omega)$ $\displaystyle=$
$\displaystyle-\frac{1}{\pi}\int dr\sum_{q=n,p}f_{q}(r)\mbox{ Im
}\delta\rho_{F,q}(r)$ $\displaystyle=$
$\displaystyle-\frac{1}{\pi}\sum_{qq^{\prime}}\int\int
drdr^{\prime}f_{q}(r)\mbox{ Im
}R_{qq^{\prime}}(r,r^{\prime})f_{q^{\prime}}(r^{\prime})$

by using the density fluctuation $\delta\rho_{F,q}(r)$ and RPA response
function $R_{qq^{\prime}}(r,r^{\prime})$ which were defined by

$\displaystyle\delta\rho_{F,q}(r)$ $\displaystyle=$
$\displaystyle\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{q}^{\mathsf{T}}\cdot\vec{\mbox{\boldmath$\phi$}}_{F}(r)$
(44) $\displaystyle R_{qq^{\prime}}(r,r^{\prime})$ $\displaystyle=$
$\displaystyle\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{q}^{\mathsf{T}}\cdot\mbox{\boldmath$\mathcal{G}$}^{(+)}(r,r^{\prime})\cdot\vec{\tilde{\mbox{\boldmath$\varphi$}}}_{q^{\prime}}$
(45)

Since $\vec{\mbox{\boldmath$\phi$}}_{F}(r)$ is $2(N_{n}+N_{p})$-dimensional
vector which is represented by

$\displaystyle\vec{\mbox{\boldmath$\phi$}}_{F}=\begin{pmatrix}\phi_{F,1}^{(n)}\\\
\phi_{F,2}^{(n)}\\\ \vdots\\\ \phi_{F,N_{n}}^{(n)}\\\ \phi_{F,1}^{(p)}\\\
\phi_{F,2}^{(p)}\\\ \vdots\\\ \phi_{F,N_{p}}^{(p)}\end{pmatrix}\hskip
5.0pt\mbox{ with }\hskip
5.0pt\phi_{F,\alpha}^{(q)}=\begin{pmatrix}X_{F,\alpha}^{(q)}\\\
Y_{F,\alpha}^{(q)}\end{pmatrix},$ (46)

the density fluctuation $\delta\rho_{F,q}(r)$ can be decomposed by each
transition component $\alpha$ as

$\displaystyle\delta\rho_{F,q}(r)$ $\displaystyle=$
$\displaystyle\sum_{\alpha=1,N_{q}}\tilde{\varphi}_{\alpha}^{(q)}\left(X^{(q)}_{F,\alpha}(r)+Y^{(q)}_{F,\alpha}(r)\right)$
(47) $\displaystyle=$
$\displaystyle\sum_{\alpha=1,N_{q}}\delta\rho_{F,q}^{(\alpha)}(r).$ (48)

Note that the capability of such a transition component decomposition is one
of the features of our RPA method using the Jost function. This is because the
existing cRPA method cannot perform a transition component decomposition of
density fluctuations because the density fluctuations or RPA response
functions are directly obtained.

##  III Numerical setup and check

The first Riemann sheet on which the bound state exists is analytically
connected to the Riemann sheet on which the pole of resonance exists by a
branch cut line extending from the branching point given by the threshold
energy on the real axis of complex energy. As will be discussed in more detail
later in Sec.III.3, there are as many Riemann sheets as there are sign
combinations of the imaginary part of the complex momentum defining the
complex energy plane, and the number of types of complex momentum is
determined by the number of transition configurations (see Table 2 for $E1$
dipole of 16O), so the heavier the nucleus, the more Riemann sheets that are
defined. Therefore, in this paper the electric dipole excitations of 16O
(known as relatively light spherical nuclei) are calculated and analyzed using
the Woods-Saxon potential for the mean field $U_{lj}^{(q)}(r)$ and simple
density-dependent interactions for the residual interactions
$\kappa_{qq^{\prime}}$.

In this section, the model and parameters used in this paper and the results
of a comparison with the cRPA method as a numerical check will be presented in
Sec.III.1 and III.2, respectively. In Sec.III.3, an explanation of the
definition of the Riemann sheet in 16O electric dipole excitations and a
numerical check of the analytic continuation is presented.

Table 1: The bound single-particle levels for neutron and proton obtained by using the Woods-Saxon potential model. The unit is MeV. |  | Neutron | Proton |
---|---|---|---|---
| $s_{1/2}$ | -36.17 | -31.16 |
| $p_{3/2}$ | -21.31 | -16.84 |
| $p_{1/2}$ | -16.38 | -11.95 |
| $d_{5/2}$ | -6.81 | -2.95 |
| $s_{1/2}$ | -4.90 | -1.43 |

###  III.1 Model and parameters

The Woods-Saxon potential model and the parameters are given by

$\displaystyle U_{lj}^{(n)}(r)$ $\displaystyle=$ $\displaystyle
V_{0}^{(n)}f_{WS}(r)+V_{1}^{(n)}\mbox{\boldmath$l$}\cdot\mbox{\boldmath$s$}\frac{1}{r}\frac{df_{WS}(r)}{dr}$
(49) $\displaystyle+\frac{\hbar^{2}l(l+1)}{2mr^{2}}$ $\displaystyle
V_{0}^{(n)}$ $\displaystyle=$
$\displaystyle-60\left(1-0.67\frac{N-Z}{A}\right)$ (50) $\displaystyle
V_{1}^{(n)}$ $\displaystyle=$ $\displaystyle
15\left(1-0.67\frac{N-Z}{A}\right)$ (51) $\displaystyle U_{lj}^{(p)}(r)$
$\displaystyle=$ $\displaystyle
V_{0}^{(p)}f_{WS}(r)+V_{1}^{(p)}\mbox{\boldmath$l$}\cdot\mbox{\boldmath$s$}\frac{1}{r}\frac{df_{WS}(r)}{dr}$
(52) $\displaystyle+\frac{\hbar^{2}l(l+1)}{2mr^{2}}+V_{C}(r)$ $\displaystyle
V_{0}^{(p)}$ $\displaystyle=$
$\displaystyle-60\left(1+0.67\frac{N-Z}{A}\right)$ (53) $\displaystyle
V_{1}^{(p)}$ $\displaystyle=$ $\displaystyle
15\left(1+0.67\frac{N-Z}{A}\right)$ (54) $\displaystyle f_{WS}(r)$
$\displaystyle=$ $\displaystyle\frac{1}{1+\exp((r-R)/a)}$ (55) $\displaystyle
R$ $\displaystyle=$ $\displaystyle r_{0}(A-1)^{\frac{1}{3}}$ (56)
$\displaystyle r_{0}$ $\displaystyle=$ $\displaystyle 1.2,\hskip 10.0pta=0.65$
(57) $\displaystyle
V_{C}(r)=\left\\{\begin{array}[]{cc}\frac{Ze^{2}}{2R}\left(3-\left(\frac{r}{R}\right)^{2}\right)&(r<R)\\\
\frac{Ze^{2}}{r}&(r\geq R)\end{array}\right.$ (60)

We adopt and use the following residual interaction model and parameters

$\displaystyle\kappa_{qq}$ $\displaystyle=$
$\displaystyle\frac{1}{2}t_{0}(1-x_{0})+\frac{t_{3}}{12}\left((5+x_{3})\rho-(2+4x_{3})\rho_{q}\right)$
$\displaystyle\kappa_{np}$ $\displaystyle=$ $\displaystyle\kappa_{pn}$ (62)
$\displaystyle=$ $\displaystyle
t_{0}\left(1+\frac{1}{2}x_{0}\right)+\frac{t_{3}}{12}(5+x_{3})\rho$

with $t_{0}=-1100$ MeV fm3, $t_{3}=16000$ MeV fm6, $x_{0}=0.5$ and
$x_{3}=1.0$.

For the external field $f_{q}(r)$ in Eq.(37), we adopt

$\displaystyle f_{n}(r)$ $\displaystyle=$ $\displaystyle e\frac{Z}{A}r$ (63)
$\displaystyle f_{p}(r)$ $\displaystyle=$ $\displaystyle-e\frac{N}{A}r$ (64)

as an $E1$ dipole operator where $A=N+Z$, with $N=Z=8$ for 16O.

Figure 1: (Color online) Comparison of the cRPA method (dotted curve) with the
method we developed in this paper (Jost-RPA, solid curve) by the $E1$ dipole
excitation strength of 16O.The strength function $S_{F}$ given as a function
of complex energy $\omega$ is plotted as a function of the real part
$\omega_{r}$ of complex energy. Note that in the upper and lower panels, the
imaginary part of complex energy $\omega_{i}=0$ and $0.02$ MeV, respectively.

In this paper, we solved the simultaneous second-order differential equation
Eq.(9) by taking into account the boundary condition of continuum in the
coordinate space up to $r\leq 20$ fm with an equidistant interval $\Delta
r=0.2$ fm.

The single particle levels of neutron and proton for both hole and particle
bound states which are obtained by using the Woods-Saxon potential are shown
in Table.1. The number assigned to the combination of angular momentum that
couples with the angular momentum of the hole state to make L=1 is the
subscription $\alpha$ shown in Table.2.

Table 2: The combination of angular momenta for $L=1$ associated with the configurations of the particle-hole excitations in the dipole excitations of 16O and the subscript $\alpha$ that designates it for neutron and proton, respectively. $q$ | $\alpha$ | $[lj]^{(q)}_{\alpha}$ | $[l_{h}j_{h}]^{(q)}_{\alpha}$ | $e^{(q)}_{\alpha}$
---|---|---|---|---
$(=n,p)$ |  |  |  | [MeV]
(neutron part)
| $1$ | $s_{1/2}$ | $p_{1/2}$ | $-16.38$
| $2$ | $d_{3/2}$ | $p_{1/2}$ | $-16.38$
| $3$ | $d_{5/2}$ | $p_{3/2}$ | $-21.31$
n | $4$ | $s_{1/2}$ | $p_{3/2}$ | $-21.31$
| $5$ | $d_{3/2}$ | $p_{3/2}$ | $-21.31$
| $6$ | $p_{3/2}$ | $s_{1/2}$ | $-36.17$
| $7$ | $p_{1/2}$ | $s_{1/2}$ | $-36.17$
(proton part)
| $1$ | $s_{1/2}$ | $p_{1/2}$ | $-11.95$
| $2$ | $d_{3/2}$ | $p_{1/2}$ | $-11.95$
| $3$ | $d_{5/2}$ | $p_{3/2}$ | $-16.84$
p | $4$ | $s_{1/2}$ | $p_{3/2}$ | $-16.84$
| $5$ | $d_{3/2}$ | $p_{3/2}$ | $-16.84$
| $6$ | $p_{3/2}$ | $s_{1/2}$ | $-31.16$
| $7$ | $p_{1/2}$ | $s_{1/2}$ | $-31.16$

###  III.2 Comparison with cRPA method

The cRPA method and the Jost function method extended within the framework of
the RPA theory developed in this paper (henceforth referred to as the Jost-RPA
method) are the same method in the sense that the RPA calculation is performed
by taking into account the boundary conditions of continuum. However, the cRPA
and Jost-RPA methods are completely different in the sense that the equations
to be solved are different.

The cRPA method calculates the unperturbed one-particle Green’s function that
satisfies the boundary conditions of the continuum and then calculates the
unperturbed response function using it. The unperturbed response function is
substituted into the Bethe-Salpeter equation, which is an integral equation
equivalent to the RPA equation, and the RPA response function is obtained
using the inverse matrix method in coordinate space. In contrast, the Jost-RPA
method calculates the perturbed one-particle Green’s function Eq.(35) using
the regular and non-regular solutions obtained by directly solving the RPA
equation given in the form of a differential equation Eq.(9), taking into
account the boundary conditions of the continuum, and the Jost function (which
is calculated using them). It is then used to calculate the RPA response
function, Eq.(45).

Since the two methods are thus different, we compared the two methods in the
$E1$ strength function. Both the cRPA and Jost-RPA methods allow the RPA
response function to be calculated as a function of continuous complex energy.
The strength function can be calculated using the RPA response function, as
shown in Eq.(LABEL:strength4). Fig.1 shows the $E1$ strength function for 16O
as a function of the real part of the complex energy ($\omega_{r}$) using the
cRPA and Jost-RPA methods. In the upper and lower panels, the imaginary part
of the complex energy is set to $\omega_{i}=0$ and $0.02$ MeV, respectively.
The comparison of the cRPA and Jost-RPA methods in the figure shows that the
results of the cRPA and Jost-RPA methods are in perfect agreement.

Note that the peaks appearing in the energy region below threshold in the
lower panel are bound states, but these are not seen in the upper panel. This
is because the poles of the bound states appear on the real axis of the
complex energy (i.e., they have no width) and cannot be plotted as a function
of $\omega_{r}$ with $\omega_{i}=0$. In contrast, the strength as background
is due to the contribution of the continuum, and the resonance peaks when the
width of the resonance (imaginary part of the pole) is very large compared to
$\omega_{i}$ are less affected by $\omega_{i}$ and the shape of the strength
function is almost unchanged. These are relatively well-known properties of
cRPA, and one can clearly and easily distinguish between the bound states and
the others.

Thus, the cRPA method can calculate the strength above the threshold,
including the width, but it cannot quantitatively calculate the individual
width of each resonance peak. In contrast, our Jost-RPA method can find the
corresponding pole (real part is the resonance energy, imaginary part is the
resonance width) for each resonance peak individually as the zeros of the Jost
function on the complex energy Riemann sheet. The results are shown in the
next section (Sec.IV).

Figure 2: (Color online) $|\det\mbox{\boldmath$\mathcal{J}$}|$ plotted as a
function of the complex energy $\omega$. The values of
$|\det\mbox{\boldmath$\mathcal{J}$}|$ with the complex energy defined as the
first Riemann sheet ($\omega_{0}$) is shown in the Im $\omega_{i}>0$ region.
$|\det\mbox{\boldmath$\mathcal{J}$}|$ with $\omega^{(1)}_{1}$,
$\omega^{(1)}_{3}$ and $\omega^{(1)}_{4}$ are shown in the Im $\omega_{i}<0$
region. And they are connected at the branch-cut lines which are existing on
the real axis continuously. The sky-blue solid, red dotted and blue dashed
curves are the branch cute lines which connect $\omega_{0}$ with
$\omega^{(1)}_{1}$, $\omega^{(1)}_{3}$ and $\omega^{(1)}_{4}$, respectively.

###  III.3 Riemann sheets of the complex energy planes

The complex energy $\omega$ can be defined as a function of the complex
momentum $k_{1,\alpha}^{(q)}$ and $k_{2,\alpha}^{(q)}$ (Eqs.(6) and (7)).
According to the basics of complex analysis, the types of on the complex
energy plane (Riemann sheets) are determined by the sign of the imaginary part
of the complex momentum and its combinations. In the case of 16O, three
momentum can be defined for each of $k_{1,\alpha}^{(q)}$ and
$k_{2,\alpha}^{(q)}$ for the proton and neutron, respectively, which means
that 4,096 types of Riemann sheets can be defined mathematically. However,
Riemann sheets in which the imaginary parts of $k_{1,\alpha}^{(q)}$ and
$k_{2,\alpha}^{(q)}$ are simultaneously negative are unphysical. Therefore,
there are a total of 128 types of physically meaningful Riemann sheets, 64
types each corresponding to the regions Re $\omega>0$ and Re $\omega<0$.
Furthermore, there are only 12 types of Riemann sheets in total, 6 each in the
regions Re $\omega>0$ and Re $\omega<0$, that are analytically connected on
the real axis of $\omega$ with the first Riemann sheet ($\omega_{0}$), where
the imaginary part of all momentum is defined as positive. The definitions of
$\omega_{0}$ and the 12 types of Riemann sheets are shown below.

$\displaystyle\omega_{0}$ $\displaystyle=$ $\displaystyle\omega(\mbox{ Im
}k_{1,\alpha}^{(q)}>0,\mbox{ Im }k_{2,\alpha}^{(q)}>0)$ (65) for all $\alpha$
and $q$

$\displaystyle\omega_{1}^{(s)}$ $\displaystyle=$ $\displaystyle\omega(\mbox{
Im }k_{s,\alpha}^{(q)}<0\mbox{ for }e_{\alpha}^{(q)}=-11.95)$ (66)
$\displaystyle\mbox{ other $\mbox{ Im }k_{1,\alpha}^{(q)}$ and $\mbox{ Im
}k_{2,\alpha}^{(q)}$ are positive},$

$\displaystyle\omega_{2}^{(s)}$ $\displaystyle=$ $\displaystyle\omega(\mbox{
Im }k_{s,\alpha}^{(q)}<0\mbox{ for }e_{\alpha}^{(q)}\geq-16.38)$ (67)
$\displaystyle\mbox{ other $\mbox{ Im }k_{1,\alpha}^{(q)}$ and $\mbox{ Im
}k_{2,\alpha}^{(q)}$ are positive},$

$\displaystyle\omega_{3}^{(s)}$ $\displaystyle=$ $\displaystyle\omega(\mbox{
Im }k_{s,\alpha}^{(q)}<0\mbox{ for }e_{\alpha}^{(q)}\geq-16.84)$ (68)
$\displaystyle\mbox{ other $\mbox{ Im }k_{1,\alpha}^{(q)}$ and $\mbox{ Im
}k_{2,\alpha}^{(q)}$ are positive},$

$\displaystyle\omega_{4}^{(s)}$ $\displaystyle=$ $\displaystyle\omega(\mbox{
Im }k_{s,\alpha}^{(q)}<0\mbox{ for }e_{\alpha}^{(q)}\geq-21.31)$ (69)
$\displaystyle\mbox{ other $\mbox{ Im }k_{1,\alpha}^{(q)}$ and $\mbox{ Im
}k_{2,\alpha}^{(q)}$ are positive},$

$\displaystyle\omega_{5}^{(s)}$ $\displaystyle=$ $\displaystyle\omega(\mbox{
Im }k_{s,\alpha}^{(q)}<0\mbox{ for }e_{\alpha}^{(q)}\geq-31.16)$ (70)
$\displaystyle\mbox{ other $\mbox{ Im }k_{1,\alpha}^{(q)}$ and $\mbox{ Im
}k_{2,\alpha}^{(q)}$ are positive},$ $\displaystyle\omega_{6}^{(s)}$
$\displaystyle=$ $\displaystyle\omega(\mbox{ Im }k_{s,\alpha}^{(q)}<0\mbox{
for }e_{\alpha}^{(q)}\geq-36.17)$ (71) $\displaystyle\mbox{ other $\mbox{ Im
}k_{1,\alpha}^{(q)}$ and $\mbox{ Im }k_{2,\alpha}^{(q)}$ are positive},$

where $s=1,2$. Riemann sheets $\omega_{1-6}^{(1)}$ are expected to be
connected with $\omega_{0}$ at $11.95\leq\mbox{ Re }\omega\leq 16.38$ MeV,
$16.38\leq\mbox{ Re }\omega\leq 16.84$ MeV, $16.84\leq\mbox{ Re }\omega\leq
21.31$ MeV, $21.31\leq\mbox{ Re }\omega\leq 31.16$ MeV, and $31.16\leq\mbox{
Re }\omega\leq 36.17$ MeV, respectively. And $\omega_{1-6}^{(2)}$ are expected
to be connected with $\omega_{0}$ in the negative energy region. We confirmed
these analytical continuation numerically in Fig.2.

In Fig.2 we show $|\det\mbox{\boldmath$\mathcal{J}$}|$ calculated using
$\omega_{0}$ in the region $\omega_{i}>0$, and in the region $\omega_{i}<0$ we
show $|\det\mbox{\boldmath$\mathcal{J}$}|$ calculated using the Riemann sheets
($\omega_{1}^{(1)}$, $\omega_{3}^{(1)}$ and $\omega_{4}^{(1)}$-planes)
connected to $\omega_{0}$ in the region $\omega_{r}<30$ MeV
($\omega_{2}^{(1)}$ is also connected with $\omega_{0}$ in a very narrow
region $16.38<\omega_{r}<16.84$ MeV but is not shown in this figure for ease
of seeing the figure).

Figure 3: (Color online) The $E1$ strength function of 16O plotted as a
function of $\omega_{r}$ with $\omega_{i}=0$ MeV in the upper panel, and
$|\det\mbox{\boldmath$\mathcal{J}$}(\omega_{i}^{(1)})|/|\det\mbox{\boldmath$\mathcal{J}$}(\omega_{0})|$
for $i\in 1,3$ and $4$ in the lower three panels. The poles labeled (a)-(g)
are poles which can contribute to the strength function. The sky-blue solid,
red dotted and blue dashed lines (with both side arrows) denote the branch cut
lines for three Riemann sheets ($\omega_{1}^{(1)}$, $\omega_{3}^{(1)}$ and
$\omega_{4}^{(1)}$) which connect with the first Riemann sheet ($\omega_{0}$).

##  IV Analysis

###  IV.1 Poles of $E1$ strength of 16O on the complex energy plane

The solution of Eq.(9) is given at the zero point of the Jost function (i.e.
the value of complex energy $\omega$ that satisfies Eq.(41)) and the bound
states appear on the real axis of complex energy below the threshold. It is
believed that the poles of the resonance state exists on the Riemann sheet
that is analytically connected to the first Riemann sheet ($\omega_{0}$-plane)
on the real axis above threshold. The real part of the resonance pole gives
the resonance energy, and the imaginary part the half-width of the resonance.
However, the poles exhibited by $|\det\mbox{\boldmath$\mathcal{J}$}|=0$ also
include unphysical ones due to hole-hole excited configurations, which are
known to cancel at the level of the unperturbed response function and do not
contribute to the strength function, and they appear on the real axis of the
complex energy. In Fig.2, such unphysical poles appear at $\omega_{r}=$ 19.21
MeV and 19.79 MeV. The presence of such unphysical poles makes it difficult to
see the physical poles that exist near the real axis when
$|\det\mbox{\boldmath$\mathcal{J}$}|$ is plotted on the complex energy plane.
Considering the property that the poles of resonances which have widths
existing above the threshold do not appear on the real axis,
$\frac{|\det\mbox{\boldmath$\mathcal{J}$}(\omega_{i}^{(1)})|}{|\det\mbox{\boldmath$\mathcal{J}$}(\omega_{0})|}$
(for $i\in 1,3,4$) is plotted instead of $|\det\mbox{\boldmath$\mathcal{J}$}|$
in Fig.3 in order to find the poles corresponding to the resonances.

The top panel of Fig. 3 shows the $E1$ strength function of 16O as a function
of $\omega_{r}$ with $\omega_{i}=0$ MeV, in order to see the correspondence
between the poles on the complex energy plane and the peaks of the strength
function. Fig.4 is an enlarged view of the region of $-0.01\leq\omega_{i}\leq
0.0$ MeV and $15.77\leq\omega_{r}\leq 15.81$ MeV of $\omega_{1}^{(1)}$-plane
near the low-lying peak in Fig.3.

Figure 4: (Color online) In Fig. 3, the poles (a), (b), and (c), which are
located very close to the real axis of the complex energy plane, are shown
enlarged.

The pole that can be moved to the first Riemann sheet by rotating the branch
cut line that analytically connects the first Riemann sheet
($\omega_{0}$-plane) to the other Riemann sheets on the complex energy around
the branch point, e.g. by complex scaling, is the pole of resonance, which
contributes as a peak to the strength function. The poles shown as (a)-(g) in
Fig.3 are the poles that are considered to contribute to the peak of the
strength function. Their values of (a)-(g) are given in the second column of
Table.3.

Each of these poles found in the complex plane is considered to be linearly
independent, and the giant electric $E1$ dipole resonance seems to be formed
by four independent poles (d), (e), (f), and (g) from Fig.3. This result shows
that, at least within the framework of RPA theory, one of the features of
giant resonances, the ’large width’, is not given by the imaginary part of a
single pole, but is formed by several independent poles. In the following
subsections, further analysis will be carried out to investigate the
properties of each of the poles.

Figure 5: (Color online) The $E1$ strength function which is calculated by RPA, RPA with $\kappa_{np}=0$, and the unperturbed response shown in the upper three panels. Trajectories of poles (a)-(g) which shows the connection between the RPA and unperturbed response, obtained by varying the constant parameters $f_{1}$ and $f_{2}$ multiplied to the residual interaction are shown in the lower panel. The filled circle ($\bullet$), square ($\square$) and triangle ($\triangle$) symbols show the position of poles which are calculated by the RPA (i.e. $f_{1}=1,f_{2}=1$), RPA with $\kappa_{np}=0$ ($f_{1}=1,f_{2}=0$) and unperturbed ($f_{1}=0,f_{2}=0$) solutions, respectively. (See text and Table 3) Table 3: The values of poles ((a)-(g)) on the complex energy planes, which are calculated by the RPA, RPA with $\kappa_{np}=0$ and the unperturbed response, which are corresponding to the filled circle ($\bullet$), square ($\square$) and triangle ($\triangle$) symbols shown in the lower panel of Fig.5. Pole | RPA | RPA with $\kappa_{np}=0$ | Unperturbed($\kappa_{qq^{\prime}}=0$) | Riemann surface |
---|---|---|---|---|---
| [MeV] | [MeV] | [MeV] |  |
(a) | $15.79-i3.35\times 10^{-3}$ | $15.48-i0.15$ | $15.47-i0.26$ |  $\omega_{1}^{(1)}$-plane |
(b) | $16.03-i0.07$ | $16.16-i0.18$ | $15.47-i0.26$ |  |
(c) | $17.15-i0.02$ | $16.59-i1.39\times 10^{-3}$ | $16.58-i1.33\times 10^{-3}$ |  $\omega_{3}^{(1)}$-plane |
(d) | $19.34-i0.53$ | $17.32-i3.64\times 10^{-2}$ | $16.58-i1.33\times 10^{-3}$ |  |
(e) | $19.28-i0.95$ | $19.32-i0.92$ | $19.32-i0.92$ |  |
(f) | $20.76-i0.34$ | $20.41-i0.27$ | $20.36-i0.26$ |  |
(g) | $21.71-i0.63$ | $21.58-i7.45\times 10^{-3}$ | $21.50-i1.33\times 10^{-3}$ |  $\omega_{4}^{(1)}$-plane |

###  IV.2 Trajectories of poles as a response to residual interactions

Multiplying the residual interactions $\kappa_{qq}$ (for both $\kappa_{nn}$
and $\kappa_{pp}$) and $\kappa_{np}$ by constants $f_{1}$ and $f_{2}$,
respectively, as $\kappa_{qq}\to f_{1}\times\kappa_{qq}$, $\kappa_{np}\to
f_{2}\times\kappa_{np}$, and varying them from 0 to 1, the trajectories of the
poles can be drawn as shown in the bottom panel of Fig.5. The top three panels
show, in order from top to bottom, the RPA strength function
($f_{1}=f_{2}=1$), the RPA strength function with $f_{1}=1$ and $f_{2}=0$
(i.e. $\kappa_{np}=0$) and the unperturbed strength function ($f_{1}=f_{2}=0$,
i.e. $\kappa_{qq^{\prime}}=0$). The positions of the poles corresponding to
the peaks of the RPA strength function, the RPA strength function with
$\kappa_{np}=0$ and the unperturbed strength function are shown by circle
($\bullet$), square ($\square$) and triangle ($\triangle$) symbols
respectively in the bottom panel. The specific values of the circle
($\bullet$), square ($\square$) and triangle ($\triangle$) for the poles in
(a)-(g) are given in columns 2, 3 and 4 of Table 3, respectively.

Focus on poles (d), (e), (f) and (g), which seem to be related to giant
resonances. It can be seen from Fig.5 that, first of all, poles (d), (e), (f)
and (g) have different origins. The pole (d) arises from the pole of the
unperturbed neutron resonance of
$\nu\left[d_{3/2}\otimes\left(p_{1/2}\right)^{-1}\right]_{1}$ at
$16.58-i1.33\times 10^{-3}$ MeV. This unperturbed neutron resonance appears as
a very sharp peak in the unperturbed strength function. The pole (e) is almost
unaffected by residual interactions and no corresponding peak in the
unperturbed strength function. This is thought to be a shape resonance created
by a mean field with a very wide width and only a small contribution as
background of the strength function. The pole (f) originates from the
unperturbed proton resonance of
$\pi\left[d_{3/2}\otimes\left(p_{3/2}\right)^{-1}\right]_{1}$ at $20.36-i0.26$
MeV and appears in the unperturbed strength function as a peak with a rather
broad width. The effect of residual interactions is small. The pole (g) arises
from the pole of the unperturbed neutron resonance of
$\nu\left[d_{3/2}\otimes\left(p_{3/2}\right)^{-1}\right]_{1}$ at
$21.50-i1.3\times 10^{-3}$ MeV. This unperturbed neutron resonance appears as
a very sharp peak in the unperturbed strength function.

The influence of the residual interactions is remarkable for poles (d) and
(g), both of which originate from the poles of the neutron unperturbed
resonance. As $\kappa_{np}$ increases, both the resonance energy (real part of
the pole) and the width (imaginary part of the pole) of the pole (d) increase
significantly. The significant shift of the resonance energy of the $E1$
dipole resonance to higher energy due to residual interactions is a typical
property of the collective excited state shown by the schematic model ring .
However, the response property of the pole (d) to the residual interaction is
not only the resonance energy (the real part of the pole) shift to higher
energy, but also the significant increase in the width (the imaginary part of
the pole). In the case of the pole (g), the response to the residual
interaction shows almost no energy shift, but a remarkable increase in the
width. This implies that this is mainly due to the dominance of coupling with
the non-resonant continuum. And this may be related to the fact that only the
pole (g) lies on a different Riemann sheet from the other resonances, and
there is no proton-originated pole on the same Riemann sheet (see Fig.3).

Poles (d) and (f) originate from the poles of the neutron and proton
unperturbed resonances, respectively, so they originally belong to different
Riemann sheets. However, the $\kappa_{np}$ effect mixes the neutron and proton
components so that (d) and (f) are finally poles belonging to the same Riemann
sheet. However, pole (f) has a very small energy shift, and its width does not
change much compared to pole (d). This may be due to the fact that the
unperturbed resonance (shape resonance) character tends to remain stronger at
pole (f) due to the Coulomb barrier.

###  IV.3 Component structure analysis of density fluctuations

In the Jost-RPA method, the density fluctuation $\delta\rho_{F,q}(r)$ is
defined by Eq.(44) and can be decomposed into components with subscription
$\alpha$, as shown in Eq.(48). Fig. 6 and 7 show the imaginary part of the RPA
density fluctuations and unperturbed density fluctuations for poles (d), (f)
and (g), respectively. The energy of the density fluctuations is the real part
of the poles, i.e., on the real axis of the complex energy plane
($\omega_{i}=0$). In each figure, the density fluctuations for protons and
neutrons are shown in the upper panel by solid curves, while the density
fluctuations for the $\alpha=2,3$ and $6$ components added together are shown
by dotted curves. In the middle and lower panels, the main contributing
components are shown when the proton and neutron density fluctuations are
decomposed into transition components, respectively. Note that the density
fluctuations shown in these figures are normalized by the value of the
strength function $S_{F}$ at a given energy ($\omega_{r}$).

Looking at the solid curves in the top panels of Fig. 6, the density
fluctuations at poles (d) and (f) both roughly show the neutron density and
proton density oscillating in antiphase, which is a typical feature of the
collective motion of isovector dipole excitations. The density fluctuations of
the pole (g) also show movement of the neutron and proton densities in
antiphase, but the proton density is likely to jump out of the nucleus.

Figure 6: (Color online) The density fluctuations $\delta\rho_{F}$ defined by
Eq.(44) for the poles (d), (f) and (g) which are expected to be related with
the giant dipole resonance are shown in the top panels by the red solid (for
proton) and blue solid (for neutron) curves. The dotted curves show the
summation of the $\alpha=2$, $3$ and $6$ components. Component decomposition
of the density fluctuation for proton and neutron are shown in the middle and
bottom panels, respectively.  Figure 7: (Color online) The same figure with
Fig.6 but calculated with the unperturbed response ($\kappa_{qq^{\prime}}=0$).

In the case of the pole (d) (originating from the neutron unperturbed
transition $\alpha=2$, i.e.
$\nu\left[d_{3/2}\otimes\left(p_{1/2}\right)^{-1}\right]_{1}$), both proton
and neutron components ($\delta\rho_{F,p}$ and $\delta\rho_{F,n}$) consist of
a superposition of in-phase $\alpha=2,3$ and $6$ components, plus $\alpha=4$
and $5$ components with antiphase. In the case of the pole (f) (originating
from the proton unperturbed transition $\alpha=5$, i.e.
$\pi\left[d_{3/2}\otimes\left(p_{3/2}\right)^{-1}\right]_{1}$), the $\alpha=5$
component in the proton component $\delta\rho_{F,p}$ is enhanced by the in-
phase superposition of the $\alpha=2,3$ and $6$ components, to which the
antiphase $\alpha=4$ component is added. The neutron components
$\delta\rho_{F,n}$) consists of a superposition of in-phase $\alpha=2,3$ and
$6$ components, plus $\alpha=4$ and $5$ components with antiphase. In the case
of the pole (g) (originating from the neutron unperturbed transition
$\alpha=5$, i.e.
$\nu\left[d_{3/2}\otimes\left(p_{3/2}\right)^{-1}\right]_{1}$), the proton
component $\delta\rho_{F,p}$ is in the form that the intranuclear part ($r<3$
fm region) of the amplitude created by the in-phase superposition of the
$\alpha=2,3$ and $6$ components is damped by the $\alpha=4$ component, leaving
mainly only the extranuclear component. The neutron component
$\delta\rho_{F,n}$ is dominated by the $\alpha=5$ component and the influence
of the other components is relatively small.

The component structure of the density fluctuations of pole (d) is roughly
antiphase to each other for neutrons and protons, whereas the density
fluctuations of pole (f) have a neutron component very similar to that of pole
(d), but the proton component does not have the antiphase structure of
neutrons. This may be because the proton component still keeps its unperturbed
resonance (shape resonance) character due to the Coulomb barrier, although
poles (d) and (f) influence each other by belonging to the same Riemann sheet.
Or it could be that pole (d) and pole (f) overlap each other due to their wide
widths, and the component structure of pole (d) just appears in the component
structure of the density fluctuations of pole (f), but actually pole (f)
itself is a shape resonance with little effect from residual interactions.

The component structure of the density fluctuations of the pole (g) is very
different from that of the poles (d) and (f). This may be due to the fact that
pole (g) belongs to a different Riemann sheet than (d) and (f). The difference
in component structure between protons and neutrons is assumed to be due to
the fact that only unperturbed resonances (shape resonances) of neutrons exist
in the same Riemann sheet. This may mean that, when the $E1$ excitation
occurs, in which the neutron and proton move in opposite phases, the neutron
tends to stay in the nucleus due to the characteristic of shape resonance,
while the proton mainly couples to the nonresonant continuum and tends to move
away from the nucleus.

These are inferences based on the analysis which we have done in this paper.
In order to clarify whether the inferences are correct, it is necessary to
analyze the pole-by-pole contribution to the strength function and density
fluctuations, excluding overlaps from other resonances, etc. However, this
requires changing the completeness of the system using the complex scaling
method aoyama  or Berggren’s method berggren , because the completeness of the
system is defined by the first Riemann sheet newton  (The contribution of
resonances is not explicitly included in the completeness). These issues are
for the future.

For the moment, the analysis at the current stage of this paper has shown that
poles (d), (f) and (g), which may be related to the $E1$ giant dipole
resonance, basically differ in their properties (resonance energy, width,
response to residual interactions, component structure of density
fluctuations) for each pole. Especially when a pole belongs to a different
Riemann sheet from the others (e.g. pole (g)), the properties are more
markedly different.

##  V Summary and perspective

In this paper, the Jost function method is extended within the framework of
RPA theory to find the complex eigenvalues of the RPA equation on the complex
energy plane (Jost-RPA method). As a first application of the Jost-RPA method,
we chose 16O electric dipole excitations and performed numerical calculations
using the Woods-Saxon potential as the mean field and simple density-dependent
interactions as the residual interactions. The cRPA and Jost-RPA methods are
identical in terms of solving the RPA problem by considering the boundary
conditions of the continuum, but the methods for solving the equations are
completely different. Therefore, we first compared both methods by calculating
the $E1$ strength function using the same potentials and residual
interactions. It was found that the result of the Jost-RPA method is in
perfect agreement with the result of the cRPA method.

In the original Jost function method, it is known that the complex energy
eigenvalues of the system are obtained as the zeros of the Jost function on
the complex energy plane that is analytically connected to the first Riemann
sheet on the branch cut line. In the Jost-RPA method, a very large number of
complex-energy Riemann sheets are defined, even for 16O $E1$ excitations.
However, only a very limited number of Riemann sheets are analytically
connected to the first Riemann sheet by a branch cut line with a branch point
at the threshold of each configuration on the real axis of complex energy. We
checked numerically whether the Jost functions computed on those Riemann
sheets actually have an analytical connection with the first Riemann sheet.
Then, the poles (complex energy eigenvalues) corresponding to the peaks of the
$E1$ strength function were successfully found numerically on the complex
energy plane in analytical connection with the first Riemann sheet.

The poles found on the complex energy plane revealed that the $E1$ giant
dipole resonance of 16O is formed by three independent poles with different
resonance energies and widths. Two of these three poles belong to the same
Riemann sheet, and only one belongs to a different Riemann sheet. Trajectory
analysis of the poles in terms of their response to residual interactions also
shows that these three poles originate from different poles of unperturbed
resonance and that the character of their response to residual interactions is
also different for each pole. It was found that response characteristics to
the residual interaction differed from pole to pole: a pole shifted towards
higher energies with increasing width as the residual interaction became
stronger, another pole was less affected by the residual interaction, and
another pole only increased in width. The component structure of the density
fluctuations is also characteristic for each pole, but the interpretation is
not clear because of the possibility of overlap due to the width of the poles
(especially for poles belonging to the same Riemann sheet). Only one thing is
clear: the poles belonging to different Riemann sheets seem to have little
influence on each other and to have quite different properties.

Based on the residue theorem in complex function theory, it is possible to
extract and show the contribution of a specific pole in the complex plane in
terms of physical quantities such as density fluctuations or strength
functions. However, the completeness of the system is defined by the first
Riemann sheet, and the resonance poles are not explicitly included in the
completeness. Therefore, in order to apply the residue theorem for more
detailed analysis, the completeness must be modified using the complex scaling
method or Berggren’s method. Further analysis by improving the Jost-RPA method
in this direction is expected to clarify points that were unclear in the
analysis of this paper.

In this paper, RPA calculations for $E1$ dipole excitations of 16O have been
performed using the Woods-Saxon potential for the mean field and simple
density-dependent interactions for the residual interactions, but other
excitation modes and other nuclei are naturally of interest. However, for a
more quantitatively reliable analysis, it is also important to perform self-
consistent calculations sagawa ; terasaki ; colo  based on effective two-body
nuclear forces, such as the Skyrme interactions SkII ; SIII ; SGII ; SkMs ;
SkP ; SLy4 ; SkO . This is because it is well known that calculations that
break the self-consistency do not satisfy the sum rule for dipole excitations
sil , and there is no phenomenological potential model for cases such as
neutron-rich nuclei. Removing spurious modes is another important issue ring ;
ripka ; lane ; rpa-fam . In this paper, an approximation to remove the
spurious mode from the external field was used, but this approximation may not
work well in some cases, so it is preferable to remove it directly from the
RPA solution. The cRPA has already developed a method to remove the spurious
mode from the RPA response function sup-crpa . A method to remove the spurious
mode directly from the solution should also be developed for the Jost-RPA
method.

We are planning to continue our research using the Jost-RPA method developed
in this paper, improving the above-mentioned issues step by step, in order to
further analyze and understand the excited structure of nuclei in more detail
from the viewpoint of the resonance poles in the near future.

## Acknowledgments

This work is funded by Vietnam National Foundation for Science and Technology
Development (NAFOSTED) under grant number “103.04-2019.329”. Tran Dieu Thuy
was funded by the Master, PhD Scholarship Programme of Vingroup Innovation
Foundation (VINIF), code VINIF.2022.TS128. This work was partially supported
by the Hue University under the Core Research Program, grant No.
NCM.DHH.2018.09.

###  .1 Proof of Eq.(35)

The Green’s theorem can lead the following equations

$\displaystyle\left[\mbox{\boldmath$\Phi$}_{0}^{(r)\mathsf{T}}(r)\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(r)}(r)-\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}_{0}^{(r)\mathsf{T}}(r)\right)\mbox{\boldmath$\Phi$}^{(r)}(r)\right]$
$\displaystyle=\frac{2m}{\hbar^{2}}\int_{0}^{r}dr^{\prime}\mbox{\boldmath$\Phi$}_{0}^{(r)\mathsf{T}}(r^{\prime})\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\mbox{\boldmath$\Phi$}^{(r)}(r^{\prime})$
(72) $\displaystyle\pm
i\mbox{\boldmath$\mathcal{K}$}^{-1}\mbox{\boldmath$\mathcal{J}$}^{(\pm)}_{0}$
$\displaystyle-\left[\mbox{\boldmath$\Phi$}_{0}^{(r)\mathsf{T}}(r)\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(\pm)}(r)-\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}_{0}^{(r)\mathsf{T}}(r)\right)\mbox{\boldmath$\Phi$}^{(\pm)}(r)\right]$
$\displaystyle=\frac{2m}{\hbar^{2}}\int_{r}^{\infty}dr^{\prime}\mbox{\boldmath$\Phi$}_{0}^{(r)\mathsf{T}}(r^{\prime})\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\mbox{\boldmath$\Phi$}^{(\pm)}(r^{\prime})$
(73)

where $\mbox{\boldmath$\Phi$}_{0}^{(r)}(r)$ and
$\mbox{\boldmath$\Phi$}_{0}^{(\pm)}(r)$ are regular and irregular solution
matrix of

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}+\mbox{\boldmath$\mathcal{U}$}\right\\}\right]\mbox{\boldmath$\Phi$}_{0}^{(r)}$
$\displaystyle=$ $0$ (74)

and

$\displaystyle\left[\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{K}$}^{2}-\left\\{-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
r^{2}}\mbox{\boldmath$1$}+\mbox{\boldmath$\mathcal{U}$}\right\\}\right]\mbox{\boldmath$\Phi$}_{0}^{(\pm)}$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$0$},$ (75)

respectively. And $\mbox{\boldmath$\mathcal{J}$}^{(\pm)}_{0}$ is the Jost
function which is defined by $\mbox{\boldmath$\Phi$}_{0}^{(r)}(r)$ and
$\mbox{\boldmath$\Phi$}_{0}^{(\pm)}(r)$. It should be noted that
$\mbox{\boldmath$\Phi$}_{0}^{(r)}(r)$, $\mbox{\boldmath$\Phi$}_{0}^{(\pm)}(r)$
and $\mbox{\boldmath$\mathcal{J}$}^{(\pm)}_{0}$ are given as the
$2(N_{n}+N_{p})\times 2(N_{n}+N_{p})$ diagonal matrix due to the absence of
the residual interaction $\mathcal{V}$ which has the off-diagonal components.

Introducing a $2(N_{n}+N_{p})\times 2(N_{n}+N_{p})$ dimensional matrix
$\mathcal{C}$ and calculating
$\mbox{\boldmath$\Phi$}^{(\pm)}(r)\mbox{\boldmath$\mathcal{C}$}\times$(Eq.(72))T
+$\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\times$(Eq.(73))T,
we can obtain

$\displaystyle\pm
i\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\mbox{\boldmath$\mathcal{J}$}^{(\pm)}_{0}\mbox{\boldmath$\mathcal{K}$}^{-1}$
$\displaystyle+\left[\mbox{\boldmath$\Phi$}^{(\pm)}(r)\mbox{\boldmath$\mathcal{C}$}\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}(r)\right)\right.$
$\displaystyle\hskip
40.0pt\left.-\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(\pm)\mathsf{T}}(r)\right)\right]\mbox{\boldmath$\Phi$}_{0}^{(r)}(r)$
$\displaystyle-\left[\mbox{\boldmath$\Phi$}^{(\pm)}(r)\mbox{\boldmath$\mathcal{C}$}\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}(r)\right.$
$\displaystyle\hskip
40.0pt\left.-\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\mbox{\boldmath$\Phi$}^{(\pm)\mathsf{T}}(r)\right]\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}_{0}^{(r)}(r)\right)$
$\displaystyle=\frac{2m}{\hbar^{2}}\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\int_{r}^{\infty}dr^{\prime}\mbox{\boldmath$\Phi$}^{(\pm)\mathsf{T}}(r^{\prime})\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\mbox{\boldmath$\Phi$}_{0}^{(r)}(r^{\prime})$
$\displaystyle+\frac{2m}{\hbar^{2}}\mbox{\boldmath$\Phi$}^{(\pm)}(r)\mbox{\boldmath$\mathcal{C}$}\int_{0}^{r}dr^{\prime}\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}(r^{\prime})\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\mbox{\boldmath$\Phi$}_{0}^{(r)}(r^{\prime})$
(76)

If we require

$\displaystyle\left[\mbox{\boldmath$\Phi$}^{(\pm)}(r)\mbox{\boldmath$\mathcal{C}$}\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}(r)\right)\right.$
$\displaystyle\hskip
40.0pt\left.-\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\left(\frac{\partial}{\partial
r}\mbox{\boldmath$\Phi$}^{(\pm)\mathsf{T}}(r)\right)\right]=-\mbox{\boldmath$1$}$
(77)
$\displaystyle\left[\mbox{\boldmath$\Phi$}^{(\pm)}(r)\mbox{\boldmath$\mathcal{C}$}\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}(r)\right.$
$\displaystyle\hskip
40.0pt\left.-\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\mbox{\boldmath$\Phi$}^{(\pm)\mathsf{T}}(r)\right]=\mbox{\boldmath$0$},$
(78)

then we can rewrite Eq.(76) as

$\displaystyle\pm
i\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\mbox{\boldmath$\mathcal{J}$}^{(\pm)}_{0}\mbox{\boldmath$\mathcal{K}$}^{-1}-\mbox{\boldmath$\Phi$}_{0}^{(r)}(r)$
$\displaystyle=\frac{2m}{\hbar^{2}}\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{C}$}^{\mathsf{T}}\int_{r}^{\infty}dr^{\prime}\mbox{\boldmath$\Phi$}^{(\pm)\mathsf{T}}(r^{\prime})\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\mbox{\boldmath$\Phi$}_{0}^{(r)}(r^{\prime})$
$\displaystyle+\frac{2m}{\hbar^{2}}\mbox{\boldmath$\Phi$}^{(\pm)}(r)\mbox{\boldmath$\mathcal{C}$}\int_{0}^{r}dr^{\prime}\mbox{\boldmath$\Phi$}^{(r)\mathsf{T}}(r^{\prime})\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\mbox{\boldmath$\Phi$}_{0}^{(r)}(r^{\prime}).$
(79)

By taking the “trace” of Eqs.(77) and (78), and applying the following basic
properties of the trace of the matrix

  1. (i)

$\mbox{Tr}[AB]=\mbox{Tr}[BA]$

  2. (ii)

$\mbox{Tr}[A+B]=\mbox{Tr}[A]+\mbox{Tr}[B]$

  3. (iii)

$\mbox{Tr}[A^{\mathsf{T}}B]=\mbox{Tr}[AB^{\mathsf{T}}]=\mbox{Tr}[B^{\mathsf{T}}A]=\mbox{Tr}[BA^{\mathsf{T}}]$,

we can notice that Eq.(77) becomes

$\displaystyle\frac{2m}{\hbar^{2}}\mbox{Tr}\left[\mbox{\boldmath$\mathcal{C}$}\mbox{\boldmath$\mathcal{W}^{(\pm)}$}\right]=\mbox{Tr}[\mbox{\boldmath$1$}],$
(80)

and Eq.(78) is a trivial identity. Therefore we find

$\displaystyle\mbox{\boldmath$\mathcal{C}$}=\frac{\hbar^{2}}{2m}\mbox{\boldmath$\mathcal{W}$}^{(\pm)-1}=\pm\frac{1}{i}\mbox{\boldmath$\mathcal{K}$}\mbox{\boldmath$\mathcal{J}$}^{(\pm)-1}.$
(81)

By inserting Eq.(81) into Eq.(79), finally we can derive

$\displaystyle\mbox{\boldmath$\Psi$}^{(\pm)}(r)$
$\displaystyle=\mbox{\boldmath$\Psi$}^{(\pm)}_{0}(r)+\int_{0}^{\infty}dr^{\prime}\mbox{\boldmath$\mathcal{G}$}^{(\pm)}(r,r^{\prime})\mbox{\boldmath$\mathcal{V}$}(r^{\prime})\mbox{\boldmath$\Psi$}^{(\pm)}_{0}(r^{\prime})$
(82)

with use of the definition of the Green’s function Eq.(35), where
$\mbox{\boldmath$\Psi$}^{(\pm)}(r)$ and
$\mbox{\boldmath$\Psi$}^{(\pm)}_{0}(r)$ are the RPA and unperturbed
“scattering” wave functions defined by

$\displaystyle\mbox{\boldmath$\Psi$}^{(\pm)}(r)$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\Phi$}^{(r)}(r)\mbox{\boldmath$\mathcal{J}$}^{(\pm)-1}$
(83) $\displaystyle\mbox{\boldmath$\Psi$}^{(\pm)}_{0}(r)$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\Phi$}^{(r)}_{0}(r)\mbox{\boldmath$\mathcal{J}$}_{0}^{(\pm)-1}.$
(84)

Since Eq.(82) is the Lippmann-Schwinger equation, it is proved that the RPA
Green function is given by Eq.(35).

## References

  * (1) G. C. Baldwin and G. S. Klaiber, Phys. Rev. 71, 3 (1947).
  * (2) B. L. Berman and S. C. Fultz, Rev. Mod. Phys. 47, 713 (1975).
  * (3) M. N. Harakeh and A. van der Woude, Giant resonances, Fundamental High-frequancy modes of nuclear excitation(Oxford University Press Inc., New York, 2001).
  * (4) P.F. Bortignon, A. Bracco and R.A. Broglia, Giant resonances, nuclear structure at finite temperature (Harwood Academic Publishers, 1998).
  * (5) D. Bohm and D. Pines, Phys.Rev. 92 609 (1953).
  * (6) K.F. Liu and N. Van Giai, Phys.Lett. 65B 23 (1976).
  * (7) P. Ring and J. Speth, Nucl. Phys. A235 315 (1974).
  * (8) P. Ring and P. Schuck, The nuclear many body problem (Springer, 1980).
  * (9) S. Shlomo and G. Bertsch, Nucl. Phys. A243, 507 (1975).
  * (10) Electric and Magnetic Giant Resonances in Nuclei, edited by J. Speth (World Scientific, Singapore, 1991).
  * (11) N. Auerbach and A. Yeverechyahu, Ann. Phys. 95, 35-52 (1975).
  * (12) R. Jost and A. Pais, Phys. Rev. 82, 840 (1951).
  * (13) K. Mizuyama, N. Nhu Le, T. Dieu Thuy, and T. V. Nhan Hao, Phys. Rev. C 99, 054607 (2019).
  * (14) K. Mizuyama, N. Nhu Le, and T. V. Nhan Hao, Phys. Rev. C 101, 034601 (2020).
  * (15) K. Mizuyama, H. Cong Quang, T. Dieu Thuy, and T. V. Nhan Hao, Phys. Rev. C 104, 034606 (2021).
  * (16) S. Aoyama, T. Myo, K. Kato and K. Ikeda, Prog. Theor. Phys, 116, 1 (2006).
  * (17) T. Berggren, Nucl. Phys. A109 (1968), 265.
  * (18) R. G. Newton, J. Math. Phys. 1 (1960), 319.
  * (19) H. Sagawa, Prog. Theor. Phys. Suppl. 142, 1 (2001).
  * (20) J. Terasaki, J. Engel, M. Bender, J. Dobaczewski, W. Nazarewicz, and M. Stoitsov, Phys. Rev. C 71, 034310 (2005).
  * (21) G. Col‘o, P. F. Bortignon, S. Fracasso, and N. Van Giai, Nucl. Phys. A788, 173c (2007).
  * (22) D. Vautherin and D. M. Brink, Phys.Rev. C 5, 626 (1972).
  * (23) M. Beiner, H. Flocard, N. Van Giai, and P. Quentin, Nucl. Phys. A238,29(1975).
  * (24) N. Van Giai, and H. Sagawa, Phys. Lett. 106B, 379 (1981).
  * (25) J. Bartel, P. Quentin, M. Brack, C. Guet, and H. B. Hakanson, Nucl. Phys. A386, 79 (1982).
  * (26) J. Dobaczewski, H. Flocard, and J. Treiner, Nucl. Phys. A422, 103 (1984).
  * (27) E. Chabanst, P. Bonche, P. Haensel, J. Meyer, and F. Shaeffer, Nucl. Phys. A635, 231 (1998).
  * (28) P. G. Reinhard, D. J. Dean, et al., Phys. Rev. C 60, 014316 (1999).
  * (29) T. Sil, S. Shlomo, B. K. Agrawal, and P.-G. Reinhard, Phys. Rev. C 73, 034316 (2006).
  * (30) J.-P. Blaizot and G. Ripka, Quantum Theory of Finite Systems (The MIT Press, Cambridge, 1986).
  * (31) A. M. Lane and J. Martorell, Ann. Phys. 129, 273 (1980).
  * (32) T. Nakatsukasa, T. Inakura and K. Yabana, Phys. Rev. C 76,024318 (2007).
  * (33) K. Mizuyama, G. Colo, Phys. Rev. C 85, 024307, 2012.