format: text

text:
# Gamow Shell Model Description of Weakly Bound Nuclei and Unbound Nuclear States

## I Introduction

The major theoretical challenge in the microscopic description of weakly bound nuclei is the rigorous treatment of both the many-body correlations and the continuum of positive-energy states and decay channels. A fully symmetric description of the interplay between scattering states, resonances, and bound states in the many-body wave function requires a close interplay between methods of nuclear structure and nuclear reactions. This mutual cross-fertilization, which cannot be accomplished without overcoming a traditional separation between nuclear structure and nuclear reaction methods, is a splendid opportunity for opening a new era in the nuclear theory of loosely bound systems.

In many respects, weakly bound nuclei are much more difficult to treat theoretically than well-bound systems . The major theoretical difficulty and challenge is the treatment of the particle continuum. For weakly bound nuclei (or for nuclear states above the particle threshold), the continuum of positive-energy states and resulting decay channels must be taken into account explicitly. As a result, many cherished approaches of nuclear theory such as the conventional shell model (based on a single-particle basis of bound states) and the pairing theory must be modified.

There are many factors which make the coupling to the particle continuum important. Firstly, even for a bound nucleus, there appears a virtual scattering into the phase space of unbound states. Although this process involves intermediate scattering states, the correlated bound states must be particle stable, i.e., they must have zero width. Secondly, the properties of unbound states, i.e., above the particle (or cluster) threshold, directly reflect the continuum structure. In addition, continuum coupling directly affects the effective nucleon-nucleon interaction.

The impact of the particle continuum was discussed in the early days of the multiconfigurational shell model (SM) in the middle of the last century. However, thanks to the tremendous success of the large-scale SM in terms of interacting nucleons assumed to be perfectly isolated from an external environment of scattering states , the continuum-related matters had been swept under the rug. An example of an impact of the continuum that goes beyond the standard SM physics is the so-called Thomas-Ehrman shift  appearing in, e.g., the mirror nuclei 13C, 13N, which is a salient effect of a coupling to the continuum depending on the position of the respective particle emission thresholds. The mathematical formulation of the problem of nuclear states embedded in the continuum of decay channels goes back to Feshbach , who introduced the two subspaces containing the discrete and scattering states. Feshbach succeeded in formulating a unified description of nuclear reactions for both direct processes in the short-time scale and compound nucleus processes in the long-time scale. As far as nuclear structure is concerned, the treatment of excited states near or above the decay threshold has been a playground of the continuum shell model (CSM) . Unfortunately, a unified description of nuclear structure and nuclear reaction aspects is much more complicated and became possible in realistic situations only at the end of the last century (see Ref.  for a recent review).

In the CSM, including the recently developed Shell Model Embedded in the Continuum (SMEC) , the scattering states and bound states are treated on an equal footing. So far, most applications of the CSM, including SMEC, have been used to describe limiting situations in which there is coupling to one-nucleon decay channels only. However, by allowing only one particle to be present in the continuum, it is impossible to apply the CSM to ‘Borromean systems’ for which $`A`$\- and $`(A`$-2)-nucleon systems are particle-stable but the intermediate $`(A`$-1)-system is not. Various approaches, including the hyperspherical harmonic method or the coupled-channel approach, have been developed to study structure and reaction aspects of three-body weakly bound nuclei . However, most of these models utilize the particle-core coupling which does not allow for the exact treatment of core excitations and the antisymmetrization between the core nucleons and the valence particles.

The reason for limiting oneself to only one particle in the continuum in the CSM has been two-fold. First, the number of scattering states needed to properly describe the underlying dynamics can easily go beyond the limit of what present computers can handle. Second, treating the continuum-continuum coupling, which is always present when two or more particles are scattered to unbound levels, is difficult. There have been only a few attempts to treat the multi-particle case  and, unfortunately, the proposed numerical schemes, due to their complexity, have never been adopted in microscopic calculations involving multiconfiguration mixing. Consequently, an entirely different approach is called for.

Recently, we formulated and tested the multiconfigurational shell model in the complete Berggren basis , the so-called Gamow Shell Model (GSM). (For application to two-particle resonant states, see also Ref. .) In this paper, GSM is applied to systems containing several valence neutrons. The single-particle (s.p.) basis of GSM is given by the Berggren ensemble, which contains Gamow states and the non-resonant continuum. The Gamow states  (sometimes called Siegert  or resonant states) were introduced for the first time in 1928 to study the $`\alpha `$resonances. Gamow defined complex-energy eigenstates $`E=E_0i\mathrm{\Gamma }/2`$ in order to describe the particle emission in the quasi-stationary formalism. Indeed, if one looks at the temporal part of such a state, which is $`e^{iE_0t/\mathrm{}}e^{\mathrm{\Gamma }t/(2\mathrm{})}`$, one notices that the squared modulus of the wave function has the time-dependence $`e^{\mathrm{\Gamma }t}`$, and one can identify $`\mathrm{}/(\mathrm{\Gamma }\mathrm{log}2)`$ with the half-life of the system.

Formally, the resonant states are generalized eigenstates of the time-independent Schrödinger equation with purely outgoing boundary conditions. They correspond to the poles of the $`S`$-matrix in the complex energy plane lying on or below the positive real axis; they are regular in origin and satisfy purely outgoing asymptotics. In the quasi-stationary approach with Gamow states, each observable $`O`$ is complex. An interpretation of these complex values has been given by Berggren : the real part of the matrix element gives the average value, while the imaginary part represents the uncertainty of the mean value. This is due to the finite lifetime of the Gamow state which implies that none of the measurements in this state can have a well-defined probability.

In the previous pilot work  we showed first applications of the GSM. In this work, we give the details of calculations and demonstrate first applications of GSM to particle distributions and transition matrix elements.

The paper is organized as follows. Section II discusses how to calculate the matrix elements in the Berggren basis. The completeness relations valid for the single-particle resonant states are briefly reviewed in Sec. III, and some numerical examples involving the Berggren set of the Woods-Saxon potential are presented. Section IV describes the GSM Hamiltonian used in our work. The extension of the completeness relations to the many-body case is described in Sec. V. Sections VI and VII contain the GCM analysis of 18-22O and 6-10He, respectively. Finally, Sec. VIII contains the main conclusions of the paper.

## II Matrix elements in the Berggren basis

Gamow functions are solutions of the Schrödinger equation which are regular at the origin and have the outgoing wave asymptotics; i.e., the radial part behaves as $`e^{ikr}`$ at large distances. In the case of the spherical one-body potential, the resonant wave function $`\varphi _{nj\mathrm{}}`$ carrying the s.p. angular momentum $`\mathrm{}j`$ can be written as a product of the usual angular part and the radial wave function $`u_{nj\mathrm{}}(r)/r`$. It is customary to introduce the notation

$`\stackrel{~}{u}_{nj\mathrm{}}(r)`$ $`=`$ $`u_{nj\mathrm{}}(r)^{}`$ (1)
$`\stackrel{~}{\varphi }_{nj\mathrm{}}`$ $`=`$ $`\varphi _{nj\mathrm{}}(u\stackrel{~}{u}).`$ (2)
For bound states, one can always introduce a phase convention which makes the radial wave function real. That is, for bound states $`\stackrel{~}{\varphi }_{nj\mathrm{}}=\varphi _{nj\mathrm{}}`$. The following discussion concerns the specific properties of the radial Gamow wave functions $`u_{nj\mathrm{}}(r)`$. Of course, one should always remember that the angular part is always present, but its treatment is standard. Consequently, in the GSM, only radial matrix elements require special attention.

### A Normalization of Gamow states and one-body matrix elements

The norm $`N_i`$ of a resonant state,

$$N_i^2=_0^+\mathrm{}u_i^2(r)𝑑r,$$
(3)
and the radial matrix elements calculated in the Berggren basis,

$$O_{if}=_0^+\mathrm{}u_f(r)O(r)u_i(r)𝑑r,$$
(4)
are diverging, but this difficulty can be avoided by means of a regularization procedure . Zel’dovich proposed to multiply the integrand of a radial matrix element by a Gaussian convergence factor :

$`u_f|O|u_i=\underset{ϵ0}{lim}{\displaystyle _0^+\mathrm{}}e^{ϵr^2}u_f(r)O(r)u_i(r)𝑑r.`$ (5)
In this expression, $`|u_f`$ and $`|u_i`$ stand for single-particle (s.p.) states, and $`O(r)`$ is a radial part of a one-body operator. Using the Zel’dovich regularization method, Berggren has shown  that the Gamow states, together with scattering states, form a complete basis. In particular, using definition (5) with $`O(r)`$=1, one can demonstrate that all Gamow states can be orthonormalized, and such orthonormalized functions can be used to calculate matrix elements. Unfortunately, the method of Zel’dovich, even though important on formal grounds, cannot be used in numerical applications due to the difficulty in approaching the limit in (5) for diverging integrals.

An equivalent and more practical procedure, justified by the apparatus of the analytical continuation, was proposed by Gyarmati and Vertse . For that, let us define the following functional on $`(\mathrm{};V_{lim})`$ :

$`F(V_o)={\displaystyle \frac{O_{ij}}{N_iN_f}},`$ (6)
where $`V_o`$ is the depth of the potential generating s.p. wave functions $`u_f`$ and $`u_i`$, $`V_{lim}`$ is the depth of the potential for which one of these functions is bound and the other one is at zero energy, and $`O(r)`$ is some analytical operator. This functional is defined in such a way because the integral converges in the domain $`(\mathrm{};V_{lim})`$. It represents the radial matrix element $`u_f|O|u_i/u_f/u_i`$ between two not necessarily normalized wave functions. Since $`u_f`$, $`u_i`$ are bound, one can make them real, $`u_f^{}=u_f`$, $`u_i^{}=u_i`$. In Ref. , the analytical continuation of $`F`$ is made using the Padé approximants. In the present work we shall refer to the technique of the complex rotation  which allows calculation of $`F`$ with $`V_o>V_{lim}`$. To see that, let us call $`f(r)`$ one of three integrands $`u_f(r)O(r)u_i(r)`$, $`u_f^2(r)`$, or $`u_i^2(r)`$ and let us take $`V_o<V_{lim}`$. Since $`f`$ is analytical on $`𝖢\text{ }`$ (see Fig. 1), then, following the Cauchy theorem, one has:

$`{\displaystyle _{C_1}}f(z)𝑑z`$ $`+`$ $`{\displaystyle _{C_2}}f(z)𝑑z`$ (7)
$`+`$ $`{\displaystyle _{C_3}}f(z)𝑑z=0.`$ (8)
Since $`f`$ decreases exponentially for Re$`[z]>0`$, the integral $`_{C_2}f(z)𝑑z0`$ if $`R_f+\mathrm{}`$. For the same reason, the integrals $`_{C_1}f(z)𝑑z`$ and $`_{C_3}f(z)𝑑z`$ converge if $`R_f+\mathrm{}`$. Consequently, for $`R_f+\mathrm{}`$ one obtains:

$`{\displaystyle _R^+\mathrm{}}f(r)𝑑r={\displaystyle _0^+\mathrm{}}f(R+xe^{i\theta })e^{i\theta }𝑑x.`$ (9)
Hence, on the interval $`(\mathrm{};V_{lim})`$, one can define $`F`$ by Eq. (6), with the norm (3) given by

$`N_i`$ $`=`$ $`\sqrt{{\displaystyle _0^R}u_i^2(r)𝑑r+{\displaystyle _0^+\mathrm{}}u_i^2(R+xe^{i\theta })e^{i\theta }𝑑x}`$ (10)
and with the matrix element (4) of the form:

$`O_{if}`$ $`=`$ $`{\displaystyle _0^R}u_f(r)O(r)u_i(r)𝑑r`$ (11)
$`+`$ $`{\displaystyle _0^+\mathrm{}}[u_f(R+xe^{i\theta })O(R+xe^{i\theta })`$ (13)
$`\times u_i(R+xe^{i\theta })e^{i\theta }dx].`$
If $`u_f`$ and $`u_i`$ are bound- or decaying-state wave functions, one can write $`k_f=|k_f|e^{i\alpha _f}`$ and $`k_i=|k_i|e^{i\alpha _i}`$. As $`u_f(z)a_f(z)e^{ik_fz}`$ and $`u_ia_i(z)e^{ik_fz}`$ when Re$`[z]+\mathrm{}`$, with $`a_f`$ and $`a_i`$ the algebraic increasing functions, the integrals defining $`F`$ converge if one takes

$$\theta >\alpha _f+\alpha _i.$$
(14)
In addition, the expression for $`F`$ is analytical because $`F`$ is a function of converging integrals of analytical functions. Square roots in Eq. (10) cause no problems because $`N_i^2`$ and $`N_f^2`$ have always a positive real part. Consequently, following the theorem of analytic continuation, Eq. (6) defines also $`F`$ for $`V>V_{lim}`$. In this way, one may calculate the radial matrix elements of resonance states which a priori are not normalizable.

### B Scattering states

Scattering states represent the non-resonant continuum and explicitly enter the completeness relations discussed in Sec III. Their asymptotic behavior at $`r+\mathrm{}`$ is:

$$u(r)C^+H_{\mathrm{},\eta }^{(+)}(kr)+C^{}H_{\mathrm{},\eta }^{()}(kr),$$
(15)
where $`H_{\mathrm{},\eta }^{(\pm )}`$ denote Hankel (or Coulomb) functions. As usual, $`u(r)`$ is normalized to the Dirac $`\delta `$distribution:

$$_0^+\mathrm{}\stackrel{~}{u}(k,r)u(k^{},r)𝑑r=\delta (kk^{}),$$
(16)
which gives

$$C^+C^{}=\frac{1}{2\pi }.$$
(17)
Knowing $`u(r)`$, $`H_{\mathrm{},\eta }^{(+)}`$, $`H_{\mathrm{},\eta }^{()}`$, and their derivatives at point $`R`$, one may determine coefficients $`C^+`$ and $`C^{}`$ up to a normalization factor by solving the set of linear equations:

$`u(R)`$ $`=`$ $`C^+H_{\mathrm{},\eta }^{(+)}(kR)+C^{}H_{\mathrm{},\eta }^{()}(kR)`$ (18)
$`u^{}(R)`$ $`=`$ $`kC^+\left[{\displaystyle \frac{dH_{\mathrm{},\eta }^{(+)}}{dz}}\right]_{z=kR}+kC^{}\left[{\displaystyle \frac{dH_{\mathrm{},\eta }^{()}}{dz}}\right]_{z=kR}.`$ (20)
Finally, the scattering state is normalized to satisfy the condition (17).

### C Matrix elements involving scattering states

The calculation of matrix elements involving the scattering states is based on the complex rotation (10,11). However, the analytical continuation should be introduced differently than for resonant states. The one-body matrix element can be written as:

$`F(k_f)`$ $`=`$ $`{\displaystyle _0^R}u_f(r)V(r)u_i(r)𝑑r+A_fA_iF_{++}(k_f)`$ (21)
$`+`$ $`A_fB_iF_+(k_f)+B_fA_iF_+(k_f)`$ (22)
$`+`$ $`B_fB_iF_{}(k_f),`$ (23)
where

* $`u_f=A_fu_f^++B_fu_f^{}`$;
* $`u_i=A_iu_i^++B_iu_i^{}`$, where $`k_i`$ in $`F`$ is fixed and, in general, $`u_i`$ can be either a bound, resonant or scattering state;
* $`F_{s_fs_i}(k_f)={\displaystyle _0^+\mathrm{}}u_f^{s_f}(R+x)O(R+x)u_i^{s_i}(R+x)𝑑x`$ with $`s_f,s_i(+,)`$.

This separation is necessary because the presence of incoming and outgoing waves in the same integral does not allow one to find a unique path in the complex plane along which the integrand decreases exponentially. Consequently, for each $`F_{s_fs_i}`$ one has to consider the domain of the complex plane where it converges, and then one performs an analytical continuation with the appropriate angle $`\theta _{s_fs_i}`$.

Certain integrals cannot be regularized in the above sense. Those include $`F_+`$ and $`F_+`$ with $`u_i=u_f`$. For $`O(r)=1`$, the integrand tends toward a constant value at $`+\mathrm{}`$, independently of the value $`\theta _+`$. This can be immediately seen for neutrons, because with $`z=R+xe^{i\theta }`$ and with $`|z|+\mathrm{}`$, the product $`u^+(z)u^{}(z)\text{const}\times e^{ikz}\times e^{ikz}=\text{const}`$, and the corresponding integral diverges. In this case, however, it is easy to see that the integral is in fact a $`\delta `$ distribution, and it can be calculated by using a discrete representation of the Dirac $`\delta `$function,

$`\delta (kk_0){\displaystyle \frac{\delta _{k,k_0}}{\mathrm{\Delta }k}},`$ (24)
with $`\mathrm{\Delta }k`$ being the discretization step in $`k`$.

## III Completeness relation involving single-particle Gamow states

There exist several completeness relations involving resonant states. As shown by Lind , they all can be derived from Mittag-Leffler theory. In the following, we briefly discuss the Berggren completeness relation  which is used in our paper. The following discussion will concern the s.p. radial wave functions corresponding to a given partial wave $`(j,\mathrm{})`$.

We begin from the completeness relation of Newton :

$`{\displaystyle \underset{n}{}}|u_nu_n|+{\displaystyle _0^+\mathrm{}}|u_ku_k|dk=1,`$ (25)
where $`|u_n`$ are the normalized bound states and $`|u_k`$ are the scattering states along the real energy axis normalized according to (16). In the basis (25), one can expand any bound state or scattering state with real energy. Unfortunately, in the presence of narrow resonances, the discretization of the real energy continuum becomes cumbersome. The complex-energy formalism of Gamow states offers a simple remedy to this difficulty.

In order to derive the completeness relation with Gamow states, one has to deform the integration contour into the complex $`k`$plane, as shown in Fig. 2.

Following the residuum theorem, one obtains:

$`{\displaystyle _0^+\mathrm{}}|u_ku_k|𝑑k`$ $`+`$ $`{\displaystyle _{L_+}}|u_k\stackrel{~}{u_k}|𝑑k`$ (26)
$`=`$ $`2i\pi {\displaystyle \underset{k_n}{}}Res\left(|u_k\stackrel{~}{u_k}|\right)_{k=k_n},`$ (27)
where $`k_n`$ are the poles of $`|u_k\stackrel{~}{u_k}|`$ lying between the real axis and the complex contour.

In general, the scattering wave function $`u_k`$ can be written as :

$`u_k(r)=\sqrt{{\displaystyle \frac{𝒥^{}(k)}{2\pi 𝒥^+(k)}}}u_k^+(r)+\sqrt{{\displaystyle \frac{𝒥^+(k)}{2\pi 𝒥^{}(k)}}}u_k^{}(r).`$ (28)
In the above expression, $`𝒥^\pm `$ stands for the Jost function :

$`𝒥^\pm (k)=u^\pm {\displaystyle \frac{du}{dr}}u{\displaystyle \frac{du}{dr}}^\pm ,`$ (29)
where $`u^\pm H_{l,\eta }^\pm (kr)`$ when $`r+\mathrm{}`$, and $`uAu^++Bu^{}`$. For bound and decaying states, $`𝒥^{}(k)`$$``$0 ($`𝒥^{}(k)=0`$ only for capturing states with the incoming wave asymptotics). Consequently, the pole of $`u_k(r)`$ corresponds to zero of $`𝒥^+(k)`$. When $`kk_n`$, i.e., when $`u_k`$ approaches the resonant state, then:

$`|u_k\stackrel{~}{u_k}|{\displaystyle \frac{𝒥^{}(k_n)}{2\pi 𝒥^+(k)}}|u_{k_n}^+\stackrel{~}{u_{k_n}^+}|.`$ (30)
The derivative of the Jost function at $`k=k_n`$ is

$`\left[{\displaystyle \frac{d𝒥}{dk}}\right]_{k=k_n}=i𝒥^{}(k_n)Reg\left[{\displaystyle _0^+\mathrm{}}u_{k_n}^{+}{}_{}{}^{2}(r)𝑑r\right].`$ (31)
As this derivative is nonzero, we have as $`kk_n`$:

$`𝒥^+(k)(kk_n)\left[{\displaystyle \frac{d𝒥^+}{dk}}\right]_{k=k_n},`$ (32)
hence

$`|u_k\stackrel{~}{u_k}|{\displaystyle \frac{1}{2\pi i(kk_n)}}|u_n\stackrel{~}{u_n}|,`$ (33)
where the normalized resonant state is:

$`u_n(r)=u_{k_n}^+(r)\left(Reg\left[{\displaystyle _0^+\mathrm{}}u_{k_n}^{+}{}_{}{}^{2}(r)𝑑r\right]\right)^{\frac{1}{2}}.`$ (34)
Finally, the residuum at $`k=k_n`$ is

$`Res\left(|u_k\stackrel{~}{u_k}|\right)_{k=k_n}={\displaystyle \frac{1}{2i\pi }}|u_n\stackrel{~}{u_n}|`$ (35)
and the completeness relation follows immediately:

$`{\displaystyle \underset{n}{}}|u_n\stackrel{~}{u_n}|+{\displaystyle _{L_+}}|u_k\stackrel{~}{u_k}|dk=1.`$ (36)
In the above equation $`|u_n`$ are the Gamow states (both bound states and the decaying resonant states between the real $`k`$axis and the complex contour). Relation (36) is the Berggren completeness relation which allows one to expand the states with complex $`k`$ inside the zone between real $`k`$axis and the complex contour. One may notice again that the resonances in Eq. (34) are normalized using the squared wave function and not the modulus of the squared wave function. This is a consequence of the analytical continuation which is used to introduce the normalization of Gamow states.

Figure 2 illustrates the ingredients entering Eq. (36). The resonant states, the poles of the $`S`$ matrix, are represented by the dots. They are divided into the bound, decaying, capturing, and antibound states (see, e.g., Refs. ). The relation (36) involves the bound and decaying states and the contour $`L_+`$ lying in the fourth quadrant of the complex-$`k`$ plane.

In practical applications, one has to discretize the integral in (36) :

$$_{L_+}|u_k\stackrel{~}{u_k}|dk\underset{i=1}{\overset{N_d}{}}|u_i\stackrel{~}{u_i}|,$$
(37)
where $`u_i(r)=\sqrt{\mathrm{\Delta }_{k_i}}u_{k_i}(r)`$ and $`\mathrm{\Delta }_{k_i}`$ is the discretization step. It follows from the definition of $`u_i(r)`$ that

$`u_i|\stackrel{~}{u_j}=\delta _{i,j},`$ (38)
and the discretized Berggren relation (36) takes the form:

$`{\displaystyle \underset{n}{}}|u_n\stackrel{~}{u_n}|+{\displaystyle \underset{i=1}{\overset{N_d}{}}}|u_i\stackrel{~}{u_i}|1.`$ (39)
This relation is formally identical to the standard completeness relation in a discrete basis and, in the same way, leads to the eigenvalue problem $`H|\mathrm{\Psi }=E|\mathrm{\Psi }`$. However, as the formalism of Gamow states is non-hermitian, the matrix $`H`$ is complex symmetric.

Up to this point, the choice of the contour in Eq. (2) has been completely arbitrary. In practice, however, one wants to minimize the number of discretization points $`N_d`$ along $`L_+`$. This can be achieved if the scattering functions on the contour (or, rather, their phase shifts) change smoothly from point to point. This condition can be met if the contour does not lie in the vicinity of a pole, especially the narrow resonant state. If this condition is met, the states appearing in the basis (39) can be naturally divided into:

* Bound states – lying on the imaginary $`k`$-axis (or negative real energy axis),
* Narrow decaying states – lying close to and below the real $`k`$-axis (or below the positive real energy axis). Those states can be interpreted as physical resonances of the system,
* Non-resonant continuum – represented by the scattering states along $`L_+`$. Physically, those are the building blocks of the non-resonant background.

These definitions can be extended to many-body states in the complex energy (or momentum) plane. In the following, we shall clearly distinguish between resonant, resonance, and non-resonant states.

The completeness relations derived above hold in every $`(j,\mathrm{})`$ channel. Consequently, in practical calculations, one has to take different contours for different partial waves. As discussed below, the choice of the contour depends on the distribution of resonant states in the complex $`k`$-plane.

In a number of papers (see, e.g., Refs. ), resonant states were applied to problems involving continuum in the so-called pole expansion, neglecting the contour integral in Eq. (36). The importance of the contour contribution was investigated in Refs.  where it was concluded that if one is aiming at a detailed description, the non-resonant contribution must be accounted for. This point will be clearly seen in several examples discussed below.

### A Completeness of the one-body Berggren basis: illustrative examples

In this section, we shall discuss examples of the Berggren completeness relation in the one-body case. The s.p. basis is generated by the spherical Woods-Saxon (WS) potential:

$`V(r)`$ $`=`$ $`V_0f(r)V_{\mathrm{so}}4\stackrel{}{l}\stackrel{}{s}{\displaystyle \frac{1}{r}}{\displaystyle \frac{df(r)}{dr}},`$ (41)
$`f(r)`$ $`=`$ $`\left[1+\mathrm{exp}\left({\displaystyle \frac{rR_0}{d}}\right)\right]^1.`$ (42)
In all examples of this section, the WS potential has the radius $`R_0`$ = 5.3 fm, diffuseness $`d=0.65`$ fm, and the spin-orbit strength $`V_{so}`$ = 5.0 MeV. The depth of the central part is varied to simulate different situations.

The complex contour corresponds to three straight segments in the complex $`k`$plane, joining the points: $`k_0=0.0i0.0`$, $`k_1=0.2i0.2`$, $`k_2=0.5i0.0`$ and $`k_3=2.0i0.0`$. The contour is discretized with a different number of points: $`n`$=60, 80, 100, 120, 160, and final results are obtained using the Richardson extrapolation method. In the examples considered in this section, we shall expand the $`2p_{3/2}`$ state, $`|u_{\mathrm{WS}}`$, either weakly bound or resonant, in the basis $`|u_{\mathrm{WS}^\mathrm{B}}(k)`$ generated by the WS potential of a different depth:

$`|u_{\mathrm{WS}}`$ $`=`$ $`{\displaystyle \underset{i}{}}c_{k_i}|u_{\mathrm{WS}^\mathrm{B}}(k_i)`$ (43)
$`+`$ $`{\displaystyle _{L_+}}c(k)|u_{\mathrm{WS}^\mathrm{B}}(k)𝑑k,`$ (44)
cf. Eq. (36). In the above equation, the first term in the expansion represents contributions from the resonant states while the second term is the non-resonant continuum contribution. Since the basis is properly normalized, the expansion amplitudes meet the condition:

$`{\displaystyle \underset{i}{}}c_{k_i}^2+{\displaystyle _{L_+}}c^2(k)𝑑k=1.`$ (45)
In all cases considered, the $`0p_{3/2}`$ and $`1p_{3/2}`$ orbitals are well bound (by $`40`$ MeV and $`18`$ MeV, respectively) and do not play any role in the expansion studied.

In the first example, we shall expand the $`2p_{3/2}`$ s.p. resonance (0.25–i0.20 MeV) of a WS potential of the depth $`V_0`$=58 MeV in the basis generated by the WS potential of the depth $`V_0^B`$=59 MeV (here the $`2p_{3/2}`$ s.p. resonance has an energy of 0.16–i0.09 MeV). The density of the expansion amplitudes is shown in Fig. 3. One can see that the contribution from the non-resonant continuum is essential even though the $`2p_{3/2}`$ basis state is a resonance. In this example, as one might expect, the contribution from the resonance state in the basis is dominant.

The second example shown in Fig. 4 deals with the case of a $`2p_{3/2}`$ state that is bound in both potentials. Here $`V_0`$=64 MeV and $`V_0^{(B)}`$=62 MeV, and the $`2p_{3/2}`$ state lies at $``$0.91 MeV and $``$0.33 MeV, respectively.

As in the previous case, the contribution from the resonant (here: bound) state in the Berggren basis dominates and the contribution of the non-resonant continuum is small although non-negligible. In this figure, one can notice the small cusp at Re\[$`k`$\]=0.2, even though the density $`c^2(k)`$ is an analytic function of $`k`$. This apparent paradox is due to the fact that $`c^2(k)`$ is plotted as a function of Re\[$`k`$\] (Im\[$`k`$\]=0). Moreover, the path in the complex plane is continuous but not derivable at $`k`$ = $`0.2i0.2`$. These two aspects contribute to the appearance of the ‘discontinuous feature’, which of course has no physical meaning.

An interesting situation is presented in Fig. 5. Here the unbound $`2p_{3/2}`$ state ($`V_0`$=59 MeV) is expanded in a WS basis containing the bound $`2p_{3/2}`$ level ($`V_0^{(B)}`$=62 MeV). Consequently, the non-resonant continuum has to supply the imaginary part of the resonance’s wave function.

The last example (Fig. 6) corresponds to $`V_0`$=64 MeV and $`V_0^{(B)}`$=59 MeV. This is the most intriguing case since one expresses a bound (real) state in the basis which contains only complex wave functions (the contribution from well bound $`0p_{3/2}`$ and $`1p_{3/2}`$ s.p. states is negligible). In this case, the non-resonant continuum annihilates the imaginary component of the $`2p_{3/2}`$ s.p. resonance contained in the basis. Indeed, in this example, the contribution from the contour is dominant.

To see the convergence of the wave function obtained by the expansion method, in Fig. 7 we show the root mean square (RMS) deviation of the calculated $`2p_{3/2}`$ wave function from the exact result as a function of the number of discretization points along the contour. While the wave functions converge fairly quickly, the convergence of complex energies is slightly slower; hence, one has to employ the Richardson extrapolation method to get an energy precision of the order of a keV.

## IV Gamow Shell Model Hamiltonian

The GSM Hamiltonian applied in this work consists of a one-body term and a zero-range two-body interaction. The spherical one-body potential was taken in a WS form (41). In our study, resonant states are determined using the generalized shooting method for bound states which requires an exterior complex scaling. The numerical algorithm for finding Gamow states for any finite-depth potential has been tested on the example of the Pőschl-Teller-Ginocchio (PTG) potential , for which the resonance energies and wave functions are known analytically. Energies of all PTG resonances with a width of up to $`90`$ MeV are reproduced with a precision of at least $`10^6`$ MeV.

Contrary to the traditional shell model, the effective interaction of CSM cannot be represented as a single matrix calculated for all nuclei in a given region. The GSM Hamiltonian contains a real effective two-body force expressed in terms of space, spin, and isospin coordinates. The matrix elements involving continuum states are strongly system-dependent, and they have to be determined for each case separately. This creates an additional difficulty, but there is also a pay-off. Namely, the resulting two-body matrix elements fully take into account the spatial extension of s.p. wave functions.

In this work, as a residual interaction we took the surface delta interaction

$$V(1,2)=V_{\mathrm{SDI}}\delta \left(\stackrel{}{r}_1\stackrel{}{r}_2\right)\delta (r_1R)$$
(46)
with the the same value of $`R`$ as in the WS potential. In our exploratory GSM calculations, we consider two cases: (i) the chain of oxygen isotopes with the inert 16O core and active neutrons in the $`sd`$ shell, and (ii) the helium chain with the inert 4He core and active neutrons in the $`p`$ shell. The parameters of the GSM Hamiltonian are summarized in Table I. The WS potentials have been adjusted to s.p. states in one-neutron nuclei 17O and 5He. For 17O, the resulting WS $`0d_{5/2}`$ and $`1s_{1/2}`$ states are bound, with s.p. energies -4.142 MeV and -3.272 MeV, respectively, and $`0d_{3/2}`$ is a resonance with the s.p. energy 0.898–i0.485 MeV. The agreement with experimental data ($`e_{5/2_1^+}^{exp}`$=–4.143 MeV, $`e_{1/2_1^+}^{exp}`$ = -3.273 MeV, and $`e_{3/2_1^+}^{exp}`$ = 0.942 MeV, $`\gamma _{3/2_1^+}^{exp}`$ = 96 keV) is excellent. The strength of the SDI has been adjusted for a given configuration space to the experimental two-neutron separation energy of 18O. Since only $`0d_{3/2}`$ is a s.p. resonance, we shall include only a $`d_{3/2}`$ non-resonant continuum. In fact, the completeness relation requires taking the non-resonant continua corresponding to all partial waves ($`\mathrm{},j`$). However, if for a given partial wave no resonances are included in the basis, the corresponding non-resonant continua can be chosen along the real momentum axis. Since, to the first order, the inclusion of these continua should only result in the renormalization of the effective interaction, they can be ignored in most cases, except for Sec. V B.

The nucleus 5He, with one neutron in the $`p`$ shell, is unstable with respect to the neutron emission. Indeed, the $`J^\pi =3/2_1^{}`$ ground state of 5He lies 890 keV above the neutron emission threshold and its neutron width is large, $`\mathrm{\Gamma }`$=600 keV. The first excited state, $`1/2_1^{}`$, is a very broad resonance ($`\mathrm{\Gamma }`$=4 MeV) that lies 4.89 MeV above the threshold. Our WS potential yields single-neutron resonances $`p_{3/2}`$ and $`p_{1/2}`$ at $`E`$=0.745–i0.32 MeV and $`E`$=2.130–i2.936 MeV, respectively. In our model space we take resonances $`0p_{3/2}`$, $`0p_{1/2}`$, and the two associated complex continua $`p_{3/2}`$ and $`p_{1/2}`$. The strength of the SDI has been adjusted for a given configuration space to the experimental two-neutron separation energy of 6He.

For the $`N`$-body problem, the Hamiltonian matrix contains one- and two-body matrix elements. The one-body part corresponds to s.p. energies of basis states and contributes only to diagonal matrix elements. In general, the calculation of two-body matrix elements is performed by splitting the radial integral into 16 terms corresponding to all different possible asymptotic conditions of s.p. wave functions. Then, each term is regularized separately by an appropriate choice of angle of the external complex scaling, cf. Sec. II.

## V Many-body completeness relation with Gamow states

The discretized basis (39) can be a starting point for establishing the completeness relation in the many-body case, in a full analogy with the standard shell-model in a complete discrete basis, e.g., the harmonic oscillator basis. In this case one has:

$`{\displaystyle \underset{n}{}}|\mathrm{\Psi }_n\stackrel{~}{\mathrm{\Psi }_n}|1,`$ (47)
where the $`N`$-body Slater determinants $`|\mathrm{\Psi }_n`$ have the form $`|\varphi _1\mathrm{}\varphi _N`$, where $`|\varphi _k`$ are resonance (bound and decaying) and scattering (contour) s.p. states. The approximate equality in (47) is an obvious consequence of the continuum discretization, similarly as in (39). Like in the case of s.p. Gamow states, the normalization of the Gamow vectors in the configuration space is given by the squares of shell-model amplitudes:

$`{\displaystyle \underset{n}{}}c_n^2=1,`$ (48)
and not to the squares of their absolute values.

In the particular case of two-particle states, the completeness relation reads:

$$\underset{i_1,i_2}{}|\varphi _{i_1}\varphi _{i_2}_J_J\varphi _{i_1}\varphi _{i_2}|1.$$
(49)
This relation can be used to calculate the two-body matrix elements.

### A Determination of many-body bound and resonance states

Before discussing completeness relations in the many-body case, let us describe the method of selecting many-body resonances. In a standard shell model, one often uses the Lanczos method to find the low-energy eigenstates (bound states) in very large configuration spaces. This popular method is unfortunately useless for the determination of many-body resonances because of a huge number (continuum) of surrounding many-body scattering states, many of them having lower energy than the resonances. A practical solution to this problem is the two-step procedure proposed in Ref. :

* In the first step, one performs the pole approximation, i.e., the Hamiltonian is diagonalized in a smaller basis consisting of s.p. resonant states only. Here, some variant of the Lanczos method can be applied. The diagonalization yields the first-order approximation to many-body resonances $`|\mathrm{\Psi }_i^{(0)}`$, where index $`i`$ ($`i=1,\mathrm{},N`$) enumerates all eigenvectors in the restricted space. These eigenvectors serve as starting vectors (pivots) for the second step of the procedure.
* In the second step, one includes couplings to non-resonant continuum states in the Lanczos subspace generated by $`|\mathrm{\Psi }_j^{(0)}`$ ($`j[1,\mathrm{},N]`$).
* Finally, one searches among the $`M`$ solutions $`|\mathrm{\Psi }_{j;k}`$, ($`k=1,\mathrm{},M`$) for the eigenvector which has the largest overlap with $`|\mathrm{\Psi }_j^{(0)}`$.

This procedure is a variant of the Davidson method. Obviously, in the search for bound states, both Lanczos and Davidson methods can be used. In the following, we shall show that the procedure outlined above allows for an efficient determination of physical states within the set of all eigenvectors of a given Lanczos subspace.

As a representative example, let us consider the cases of 18O and 20O with the core of 16O, i.e., two- and four-particle systems, respectively. Here we employ the complex $`d_{3/2}`$ contour corresponding to two segments defined by the points: $`k_1=0+i0`$, $`k_2=0.2i0.05`$, and $`k_3=0.4i0.0`$. The contour is discretized with 17 and 9 points for 18O and 20O, respectively. From all the possible many-body configurations, we only keep the Slater determinants with an energy less than 5 MeV and with a width smaller than 1.7 MeV.

The results of calculations for the $`0^+`$ states in 18O are displayed in Fig. 8.

In this case, one obtains two bound states. An eigenstate lying just close to the real energy axis just above the two-neutron threshold is a candidate for a resonance. Even though the width of this state is very small, only the overlap with the states calculated in the pole approximation can give the answer. Figure 9 shows

the overlap of the ground-state wave function $`\mathrm{\Psi }_0`$ of 18O, calculated in the pole approximation, with all the $`0^+`$ eigenstates of the GSM Hamiltonian resulting from the full diagonalization. One can see that only one state (the GSM ground state) has a significant overlap with $`\mathrm{\Psi }_0`$; hence, the identification of the ground state wave function is unambiguous. Figure 10

illustrates a more challenging case of the third $`0_3^+`$ state of 18O. In spite of the fact that this state is embedded in the non-resonant $`0^+`$ continuum, its identification is straightforward. It is interesting to notice that those GSM states in Fig. 8 that represent the non-resonant background tend to align along regular trajectories. As discussed in Refs. , the shapes of these trajectories directly reflect the geometry of the contour in the complex $`k`$-plane. In the two-particle case, this information can be directly used to identify the resonance states.

Figure 11 shows the results of calculations for the $`0^+`$ states in 20O. As compared to the 18O case, the number of many-body states is much larger and the regular pattern of non-resonant states reflecting the structure of the contour is gone (the figure represents the projection of four-dimensional trajectories onto two dimensional space). While the two lowest (bound) states can be simply identified by inspection, for the higher-lying states it is practically impossible to separate the resonances from the non-resonant continuum. However, the procedure outlined above makes it possible to identify unambiguously the many-body resonance states. On the other hand, the method proposed in Refs.  cannot be easily applied.

The many-body resonances should be stable with respect to small deformations of the contour (as the physical solutions should not be dependent on the deformation of the basis). This observation offers an independent criterion for identifying resonance states. Figure 12 shows the effect of a small deformation of the contour on the stability of selected $`0^+`$ states in 20O. As expected, only the states which have previously been identified as resonances are stable with respect to small changes of the contour; the states belonging to the non-resonant continuum ‘walk’ in the complex energy plane following contour’s motion.

### B Completeness of the many-body Gamow basis: example of two interacting particles

In this section, we shall discuss the completeness of the many-body basis spanned by Gamow states. Since the number of configurations is growing extremely fast with the number of valence nucleons, and this is enhanced by including the non-resonant continuum, we shall restrict our discussion to the case of two valence particles. By analyzing the behavior of the wave function as the number of basis states grows, one can assess the impact of various truncations in the valence space. This is especially important for calculations with zero-range forces, such as the SDI interaction employed in this work, which require an energy cut-off. In the calculations contained in this section, the configuration space consists of all the Slater determinants of energy less than 35 MeV.

As a first example, we shall consider the convergence of the ground-state energy of 18O with the increasing size of the non-resonant phase space. The GSM s.p. space consists of the $`0d_{5/2}`$, $`1s_{1/2}`$ orbitals and the $`0d_{3/2}`$ Gamow resonance. This discrete basis is supplemented by adding successive continua: $`s_{1/2}`$, $`p_{1/2}`$, $`p_{3/2}`$, …, in the decreasing order of their importance. Since $`0d_{3/2}`$ is a resonance, a $`d_{3/2}`$ contour should be complex, and we take it according to Sec. V A. Other non-resonant continua are real and for their path we choose a straight segment with $`0<k<1.3`$ fm-1 . The discretization of each continuum, whether real or complex, is made with 10 points.

As one can see in Fig. 13 (top), the convergence is achieved with the $`d_{3/2}`$, $`f_{7/2}`$, $`d_{5/2}`$, $`s_{1/2}`$, and $`f_{5/2}`$ contours; the contributions from all the remaining partial waves with $`\mathrm{}>3`$ and $`\mathrm{}=1`$ are practically negligible.

As a rule of thumb, the continua which contribute most to the binding energy are those which are associated either with s.p. resonances in the basis or with weakly bound s.p. states. This is the case with the $`0d_{5/2}`$, $`1s_{1/2}`$, and $`0d_{3/2}`$ orbitals, and also with the $`0f_{7/2}`$ state which has an energy of +7.43 MeV and a width of 3.04 MeV. Even though this wide resonance is absent in the basis, it has nevertheless an indirect influence through the $`f_{7/2}`$ contour which ‘keeps memory’ of its presence. One may notice that the presence of the (complex) $`d_{3/2}`$ contour affects not only the real part of the ground-state energy but also its imaginary part; the ground-state energy has a spurious width of $``$ –130 keV in the pole approximation, and this width is reduced to $``$ –2 keV when the $`d_{3/2}`$ (10-point) contour is added. The spurious width of the ground state remains stable when other (real) contours are added.

The structure of the ground-state wave function of 18O calculated with the full non-resonant continuum discussed in the context of Fig. 13 is given in Table II. One can see that the configuration with 2 neutrons in the $`0d_{5/2}`$ shell dominates; the remaining configurations, including those with one and two neutrons in the non-resonant continuum, contribute with $``$ 15% to the wave function. Imaginary parts of squared amplitudes are generally very small. This is due to the small width of the $`0d_{3/2}`$ resonance included in the basis.

As a second example, we shall investigate the energy convergence for the weakly bound ground state of 6He. We shall assume that the structure of 6He can be described in the full $`0p`$ shell with two valence neutrons. Since the GSM s.p. space consists of the $`0p_{3/2}`$, $`0p_{1/2}`$ Gamow resonances, the associated $`p_{3/2}`$ and $`p_{1/2}`$ non-resonant continua should be complex. Here, the $`p_{3/2}`$ contour consists of three straight segments connected at points: $`k_0=0.0i0.0`$, $`k_1=0.3i0.35`$, $`k_2=0.5i0`$, and $`k_3=0.7i0`$. In the $`p_{1/2}`$ non-resonant channel, straight lines join the points: $`k_0=0.0i0.0`$, $`k_1=0.4i0.5`$, $`k_2=0.5i0`$, and $`k_3=0.7i0`$. For all other contours we take a segment \[0:0.7\] fm-1 of the real $`k`$axis. The $`p_{3/2}`$ contour is discretized with 20 points. All other contours, including the complex $`p_{1/2}`$contour, are discretized with 10 points.

As seen in Fig. 13 (bottom), the full energy convergence is attained with $`p_{3/2}`$ and $`p_{1/2}`$ contours; other scattering waves with $`\mathrm{}>1`$ are negligible. Coupling to the non-resonant $`p_{3/2}`$ and $`p_{1/2}`$ continua changes not only the ground-state energy but also its (spurious) width. In the pole approximation, the calculated ground state has a huge width of $``$ –2 MeV, which is reduced to $``$ –650 keV when the $`p_{3/2}`$ contour is added, and it reaches $`\mathrm{\Gamma }`$$``$–10 keV when both $`p_{3/2}`$ and $`p_{1/2}`$ non-resonant continua are included. Contrary to the case of 18O, the pole approximation is totally unreliable, and it does not even give a rough approximation to the energy and wave function of the ground state of 6He. The main reason is that in this case one attempts to describe the bound state using the basis which does not contain any bound state and, therefore, the non-resonant continuum is essential for compensating the resonance contribution. The analogous situation, discussed in the context of a one-body problem, can be found in Sect. III A (Fig. 6).

The structure of the ground-state wave function of 6He including all the ($`\mathrm{},j`$) continua shown in Fig. 13, is shown in Table II. One can see that the configuration with two neutrons in the $`0p_{3/2}`$ shell dominates, though the imaginary part of the corresponding squared amplitude is almost equal in magnitude to the real part. The amplitude of the $`0p_{1/2}^2`$ configuration is small, whereas the contributions from one and two particles in the non-resonant continuum, $`L_+^{(1)}`$ and $`L_+^{(2)}`$, are almost equally important.

## VI GSM Study of oxygen isotopes

In this section we shall discuss the GSM results for a chain of oxygen isotopes with several valence neutrons ($`N_{val}2`$). In the calculations, we assume the core of 16O. As discussed in Sec. IV, the valence neutrons are distributed over the $`1s_{1/2}`$ and $`0d_{5/2}`$ bound shells, the $`0d_{3/2}`$ Gamow resonance, and the discretized $`d_{3/2}`$ non-resonant continuum. The complex $`d_{3/2}`$ contour corresponds to two segments defined by the points: $`k_1=0+i0`$, $`k_2=0.2i0.05`$, and $`k_3=0.4i0.0`$. It is discretized with 9 points. Consequently, the discretized s.p. GSM space (39) consists of a total of 12 subshells on which valence neutrons are distributed. Let us reiterate that we also introduce the cut-off in the configuration space of the GSM. In this study, we shall include all Slater determinants with an energy (width) smaller than 5 MeV (1.7 MeV). Moreover, we shall take only Slater determinants with, at most, 2 neutrons in the non-resonant continuum. We have checked that configurations with a larger number of neutrons in the non-resonant continuum are of minor importance in all nuclei from 18O to 22O, and their contribution is $`<`$ 0.01% in all cases studied. Our aim is not to give the precise description of actual nuclei (for this, one would need a realistic Hamiltonian and a larger configuration space), but rather to illustrate the method, its basic ingredients, and underlying features.

#### 1 Ground states of oxygen isotopes

According to our calculations (see, e.g., Table III for 18O), the ground-state wave functions of 18-22O are dominated by the single shell model configuration $`0d_{5/2}^{N_{\mathrm{val}}}`$ having a weight of 80-95%; the non-resonant continuum contribution is relatively small. A one-neutron continuum provides between 0.1 and 1% of the wave function, whereas a two-neutron continuum yields a contribution which varies from $`0.01`$ to $`0.1`$% in these nuclei. This means that the main effect of the non-resonant continuum in these states is to cancel a spurious width induced by the $`0d_{3/2}`$ resonance included in the basis, and their influence on the real part of the energy can be neglected.

It is interesting to compare the squared amplitudes for 18O in Tables II and III. The values in Table II have been obtained by taking all the contours from $`s_{1/2}`$ to $`h_{11/2}`$; whereas the results in Table III have been calculated with only a $`d_{3/2}`$ contour. Increasing the non-resonant space gives rise to the depletion of the occupation of bound shells.

The predicted one-neutron separation energies $`S_n`$ for oxygen isotopes are displayed in Fig. 14 and compared with the data. The calculations reproduce experimental separation energies, including the magnitude of the odd-even staggering. The difference between GSM and the data are in the worst case $`500`$ keV.

#### 2 Level schemes of oxygen isotopes

Shown in Figs. 15 \- 17 are calculated and experimental spectra of 18-20O. Let us stress that the main purpose of this comparison is to demonstrate the internal coherence of GSM results and the numerical feasibility of GSM in realistic applications. Nevertheless, in spite of the simplicity of the Hamiltonian, the overall agreement between the GSM and experimental spectra is quite reasonable.

Figure 15 illustrates the case of 18O. One can see that the first three lowest states are relatively well described. Moreover, the one-neutron emission threshold is reproduced. In GSM, this threshold is obtained from the difference of calculated ground-state energies of 18O and 17O. It is important to note that all states above the one-neutron emission threshold are predicted to be resonances, i.e., their widths are positive, whereas all states below this threshold are calculated to have small spurious negative widths. This feature is certainly not imposed on the model and it proves both an internal consistency of the GSM as well as a correct numerical account of the completeness for many-body Gamow states. One should notice that in the calculations with a discretized continuum, it is impossible to assure that the bound states have a width exactly equal to zero. Instead, one expects that this numerical artifact of continuum discretization becomes numerically unimportant with increasing the number of points on the contour (see Table 1 of ). With 10-point discretization, the spurious width remains $`5`$ keV for all bound states of 18O, and we consider this a reasonable precision. Similarly, for resonances lying very close to the lowest particle emission threshold, hence having narrow widths, it may happen that their width results from the calculations as negative if the many-body completeness relation is numerically not well satisfied (see Fig. 11 and related discussion).

The calculated level scheme of 19O (see Fig. 16) reproduces the main experimental features. The ground state is $`5/2^+`$, as observed, but the higher-lying states $`3/2_1^+`$ and $`1/2_1^+`$ are given in the reversed order. Nevertheless, the next four states are rather well reproduced, in spite of the inversion of the $`5/2_2^+`$ and $`3/2_2^+`$ levels.

Finally, the GSM prediction for 20O is shown in Fig. 17. Here, the overall agreement between calculations and experiment is best for all the isotopes studied.

#### 3 Distribution of valence particles

Radial features of the density distribution of valence nucleons are always of great theoretical interest. They determine the nuclear size, electromagnetic transition rates, polarization charges, and many other nuclear properties. To assess the radial extent of the nucleonic distribution obtained in the GSM, we investigate the monopole form factor defined through the normalized radial Gamow wave functions $`u_i(r)`$ (34):

$$\rho (r)=\underset{i}{}n_i|u_i(r)\stackrel{~}{u}_i(r)|,$$
(50)
where $`n_ia_n^+a_n`$ is the GSM occupation coefficient of the s.p. Gamow orbital $`i`$. By definition, the form factor is normalized to the number of valence particles:

$$_{r=0}^{\mathrm{}^{}}\rho (r)𝑑r=N_{\mathrm{val}},$$
(51)
where the “prime” sign indicates that the integral is calculated by means of complex scaling.

Figure 18 shows the neutron form factors for the ground state and the highly excited resonance state $`2_5^+`$ of 18O.

As for all observables associated with bound states, the neutron form factor in the ground state of 18O is real even though the $`0d_{3/2}`$ resonant state is included in the basis. The imaginary part of $`\rho (r)`$ is very small; its value is less than $`3\times 10^5`$ for all values of $`r`$. On the other hand, for the $`2_5^+`$ resonance, having a width of $``$150 keV, the form factor is always complex. This feature is not an artifact of the discretization procedure. In fact, the imaginary part of the form factor has an interesting physical interpretation. Namely, it is associated with the uncertainty in the determination of the mean value (given by the real part) . In other words, in the decaying state the particle distribution cannot be given with an unlimited precision.

#### 4 Electromagnetic transition probabilities

In the shell-model calculations with Gamow states, only radial matrix elements are treated differently as compared to the standard shell model. This means that the electromagnetic (EM) transition selection rules and the angular momentum and isospin algebra do not change. To calculate the EM transitions, one can no longer use the long wavelength approximation because of the presence of the non resonant continuum. Indeed, for the diagonal EM matrix elements $`n_il_ij_ik_i|O|n_fl_fj_fk_f`$ ($`k_i`$ = $`k_f`$) between the scattering states, the complex scaling cannot be carried out (see discussion around Eq. (24)). Furthermore, since in the long wavelength approximation the EM operators behave like $`r^\lambda `$, one has to deal with derivatives of delta functions, which is difficult to handle. Without the long wavelength approximation, however, these matrix elements become finite, because it is always possible to carry out a complex scaling with the Bessel function of the photon $`j_L(qr)`$, as $`q0`$. Moreover, as they represent a set of measure zero, the diagonal non-resonant EM matrix elements can be put to zero in the discretized calculation. As all the other matrix elements can be regularized, the EM matrix elements are all well defined.

Tables IV and V display the selected EM transition rates in 18O and 19O calculated in GSM with and without the contribution from the non-resonant $`d_{3/2}`$ continuum. In all calculations, we have taken the effective neutron charge $`\delta e_n=0.5e`$. While the radial form factor discussed in Sec. (VI 3) determines the structure of monopole operators, the EM probabilities probe the off-diagonal matrix elements; hence the continuum coupling manifests itself differently. As seen in Tables IV and V, the non-resonant continuum plays an important role both for transitions involving bound states only and for transitions involving an unbound state or states. In all cases, the real part of the transition rate is reasonably approximated in the pole approximation.

Similarly, as in the problem of spurious negative widths of bound states, the imaginary part of EM transitions between bound states must disappear in the limit of the non-discretized continuum. Indeed, in our calculations the imaginary part of such a transition rate tends to zero when the coupling to the non-resonant continuum is added. In the cases studied, taking only 10 points on the $`d_{3/2}`$ contour decreases the imaginary part by one to two orders of magnitude, depending on the transition. The real part of the transition probabilities changes by up to $``$ 10% when the contribution from the contour is considered. This is nevertheless an important change if one notices that the configurations with one and two neutrons in the non-resonant continuum amount to $``$ 0.1% of bound wave functions in these nuclei. For most transitions, the experimental values are reproduced within a factor of 3 by the GSM, and we find this agreement satisfactory in light of a very simple Hamiltonian employed.

Table V displays transition rates in 19O, indicated in Fig. 16 by arrows. The three lowest transitions are between bound states and the GSM calculations with contour predict a very small imaginary part (less than 0.001 W.u.). The next four transitions in Table V involve the unbound $`7/2_2^+`$, $`1/2_2^+`$, and $`5/2_3^+`$ levels lying above the one-neutron threshold, and the corresponding transition rates are, in general, complex. The effect is particularly pronounced for the $`1/2_2^+\stackrel{\text{E2}}{}5/2_3^+`$ transition between unbound states. As mentioned in Sec. VI 3, the imaginary part gives the uncertainty of the average value . In all cases, the real part of the matrix element is slightly influenced by the interference with the non-resonant background,

## VII GSM study of helium isotopes

A description of the neutron-rich helium isotopes, including Borromean nuclei 6,8He, is a challenge for the GSM. 4He is a well-bound system with the one-neutron emission threshold at 20.58 MeV. On the contrary, as discussed in Sec. IV, the nucleus 5He is a broad resonance. The two-neutron system 6He, on the contrary, is bound with the two-neutron emission threshold at 0.98 MeV and one-neutron emission threshold at 1.87 MeV. The first excited state $`2_1^+`$ at 1.8 MeV in 6He is neutron unstable with a width $`\mathrm{\Gamma }`$=113 keV.

The s.p. configuration space used includes both resonances $`0p_{3/2}`$, $`0p_{1/2}`$ and the two associated complex continua $`p_{3/2}`$ and $`p_{1/2}`$, which are discretized with 5 points each. The fact that the resonances in the basis are so broad requires particular attention when selecting points along the contour. They are given in Table VI.

In this case, we cannot significantly increase the density of points along the contour because the calculation for the heavier helium isotopes with more valence neutrons would not be feasible. On the other hand, we do not introduce any restriction on energies and widths of Slater determinants included. Also, no restriction on the number of neutrons in the non-resonant continuum is imposed, with the only exception being 9He, where we allow for at most 4 neutrons to occupy “contour shells”. Contrary to the discussion in Sect. V B, we neglect all remaining real continua in the present calculation. We shall see that this will have an impact on the relative weight of different configurations in the wave function.

#### 1 Spectra of helium isotopes

The energies of the lowest GSM states of helium isotopes are shown In Fig. 19, and the structure of their ground state wave functions is given in Table VII.

As seen in Table VII, the non-resonant continuum contributions $`L_+^{(n)}`$ are always essential, and, in some cases (e.g., 8,9He), they dominate the structure of the ground-state wave function. Moreover, as can be seen in the example of 8He, configurations with many neutrons in the non-resonant continuum are essential for fulfilling the completeness relation. In this particular case, the $`L_+^{(1)}`$ contribution is even more important than the contribution from the resonant states, and even the $`L_+^{(4)}`$ configuration gives a non-negligible contribution of the order of 1-2%. ($`L_+^{(n)}`$ amplitudes seem to decrease with $`n`$.) Without the contour, the predicted ground-state energy of 8He is +2.08 MeV and the spurious width is huge, $`\mathrm{\Gamma }=4.16`$ MeV. The inclusion of scattering states lowers the binding energy to $`1.6`$ MeV. Clearly, the pole approximation fails miserably in this case.

It is instructive to compare the wave function decomposition in Tables VII (only $`p_{3/2}`$ and $`p_{1/2}`$ contours are included) and II (all contours up to $`f_{7/2}`$ are included). As expected, the spread of the ground-state wave function over non-resonant continuum states is larger in the latter case, but this effect is relatively less important than in the case of 18O (cf. Fig. 13).

Our calculations reproduce the most important feature of 6He and 8He: the ground state is particle-bound, despite the fact that all the basis states lie in the continuum. In spite of a very crude Hamiltonian, rather limited configuration space, etc., the calculated ground state energies shown in Fig. 19 reproduce surprisingly well the experimental data. The odd-$`N`$ isotopes of 7,9He are calculated to be neutron resonances. The neutron separation energy anomaly (i.e., the increase of one-neutron separation energy when going from 6He to 8He) is reproduced. One can see that 8He is predicted to be better bound than 6He, though, contrary to the data, the two-neutron separation energy of 8He is smaller than that of 6He. Also, the energies of excited $`2_1^+`$ states are in fair agreement with the data.

The present calculations for helium isotopes have a systematic tendency to underbind. The spurious width of the ground state of 8He is $``$ 100 keV, largely due to the $`0p_{3/2}`$ resonance. For that reason, the removal of the spurious width (which reflects the accuracy of our calculations) cannot be done here as precisely as in the oxygen case. In order to check the stability of the results for 6,8He, we carried out additional calculations increasing the number of points along the complex contour. We have found that the value of the two-neutron separation energy difference, $`S_{2n}(^8`$He) - $`S_{2n}(^6`$He), which is negative with 5 states in each non-resonant continuum, depends sensitively on the number of discretization points. This is because $`S_{2n}(^8`$He) increases fast with the size of the non-resonant phase space. While with our schematic Hamiltonian we were unable to find a fully satisfactory result for the ground-state energy of 8He just by increasing the number of discretization points, we can see that there is a systematic improvement when the number of non-resonant shells increases. Therefore, one would be tempted to associate the so-called ‘helium anomaly’ in the position of two-neutron emission threshold when going from 6He to 8He with the strong enhancement in the occupation of non-resonant continuum states.

#### 2 Distribution of valence neutrons in 6He

Form factors of valence neutrons (50) in $`0_1^+`$ and $`2_1^+`$ states of 6He are shown in Figs. 20.

In this calculation, we have taken discretized $`p_{1/2}`$ and $`p_{3/2}`$ continua with 20 points each. Even with this very high number of continuum states, the imaginary part of the form factor in the ground state of 6He is still $``$ 2$`\times 10^3`$ fm-1. It is interesting to see that the radial distribution of neutrons in this case extends well above the range of the one-body potential, as expected for a weakly bound system. In the case of the first excited (resonance) $`2_1^+`$ state, $`\rho (r)`$ is complex; real and imaginary parts of the density are of a comparable size, which is obviously related to the large width ($``$ 500 keV) of this state.

## VIII Summary

This paper contains the detailed description of the Gamow Shell Model, which can be viewed as a straightforward extension of the standard diagonalization shell model that allows for a consistent treatment of bound states and the particle continuum, including both resonances and the non-resonant background. Our first application, based on a realistic, albeit simple, two-body SDI force, concerns many-neutron systems, including states and nuclei near the neutron-emission threshold.

In this work, we succeeded in overcoming several obstacles which traditionally plagued previous continuum shell model applications. In addition to the successful inclusion of the continuum-continuum coupling by means of the complex rotation technique (exterior complex scaling), we incorporated the non-resonant part of the continuum. This has been achieved by discretizing the contour in the complex $`k`$-plane for each partial wave. Another problem which has been solved in our study is the isolation of resonance states. As a result of the GSM diagonalization, one obtains a multitude of states corresponding to the many-body continuum, some being resonances and some representing the non-resonant background. Our work offers a simple prescription on how to identify the resonance states from the multitude of complex-energy eigenstates of the GSM Hamiltonian.

The GSM Hamiltonian is composed of a one-body term and a two-body residual interaction which is directly written in terms of space, spin, and isospin coordinates. The two-body matrix elements have to be determined separately for each case by means of radial integration. As a result, the resulting Hamiltonian matrix fully takes into account the continuum coupling, in particular the spatial extension of s.p. wave functions, determining the physics of halo nuclei.

We demonstrated that the contribution from the non-resonant continuum is essential, especially for unbound and near-threshold states. In some cases (e.g., 8,9He) non-resonant continuum components dominate the structure of the g.s. wave function. The so-called ‘pole approximation’ (resonant state expansion) breaks down in such cases. In addition, the inclusion of the non-resonant continuum also impacts results for bound states.

The results of our first calculations for binding energies, spectra, electromagnetic matrix elements, and nucleonic distributions are very encouraging. In particular, pairing correlations due to the continuum-continuum scattering can bind the ground states of 6,8He with a completely unbound basis provided by the s.p. resonances of 5He. The ‘Helium anomaly’ (an increase in one- and two-neutron separation energy when going from 6He to 8He) is explained in terms of the neutron scattering to the non-resonant background. In all cases considered, our calculations yield neutron resonances above the calculated neutron threshold – a property that is not guaranteed a priori by the formalism.

Other applications of GSM, including the case of open proton and neutron shells, also employing more realistic effective interactions, are in progress. We are also working on the optimization of the non-resonant part of calculations (choice of the contours, distribution of discretization points, etc.). We are convinced that the Gamow Shell Model will become a very useful theoretical tool unifying structure and reaction aspects of weakly bound nuclei.

###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy under Contract Nos. DE-FG02-96ER40963 (University of Tennessee) and DE-AC05-00OR22725 with UT-Battelle, LLC (Oak Ridge National Laboratory).


text:
# Characterisation of the dynamical quantum state of a zero temperature Bose-Einstein condensate

## 1 Introduction

The realisation of Bose-Einstein condensation in dilute gases has made it possible to study the dynamics of quantum fields directly . In this paper, we address one of the most intriguing issues surrounding a Bose-Einstein condensate (BEC), namely the nature of its quantum state. Number states might seem to be the natural choice given that atoms cannot be created or destroyed at these energies. However, it is well known that, in an open environment, particles may be added or removed from the condensate, and this implies that the condensate is not in a pure Fock state. On the other hand, it has been suggested that there should be a highly entangled state of the condensate plus the environment . The coherence properties of the condensate have been demonstrated experimentally, indicating that a coherent state might be the most robust representation . This point turns out to be very important, as it is well-known that a coherent state propagating through an amplitude-dispersive medium evolves into a superposition of two coherent states . The generation of these Schrödinger cat-like states is considered to be a general property of nonlinear systems where dissipation effects are not very large. Recently, some schemes to create Schrödinger cat states in BECs have been proposed  and the effects of loss have also been considered.

It has also been pointed out that the state of the condensate might be squeezed, given that binary collisions between atoms play a very important role in the description of BECs  and are known to give rise to squeezing . We note that the question of describing the quantum state of a trapped BEC was addressed in reference  where a single mode approximation was considered, based on a symmetry-breaking picture. The authors showed that for repulsive interactions the state is number squeezed. In this paper, we use a multimode description in the Hartree approximation. In particular, we show how the state evolves under ballistic expansion, which is important as this is the basis of many experiments involving the condensates .

To study the evolution of the quantum state of a BEC, we have used the Heisenberg equation of motion obtained from the many-body Hamiltonian describing the condensate, and used the Hartree approximation, under which the $`n`$-atom wavefunction is written as a product of $`n`$ single atom wavefunctions. This results in a nonlinear Schrödinger equation for the condensate wavefunction, or the well-known Gross-Pitaevskii equation (GPE). The GPE has proved to be a very useful tool for studying the condensate dynamics, and its predictions are in very good agreement with experimental results at low temperatures . It is important to note that the GPE can be understood in terms of either coherent states or number states. Here, we make use of a number state description in order to obtain an analogue of the scenario presented in  for the case of light. We shall take “snapshots” of the wavefunction at various times and analyse them using the $`Q`$-function, which is a quantum-mechanical phase space function. We also want to mention that we are treating the condensate in a zero temperature regime; this means that we are basically referring to the predominant mechanisms occurring in the ground state of the condensate. At finite temperature, the mechanism is affected by the presence of excited states and therefore the dynamics are changed .

The paper is organised as follows: in Section 2, we present the Hamiltonian that describes the BEC and solve the equation of motion in the Schrödinger picture using the Hartree approximation. In Section 3, we calculate the $`Q`$-function for the condensate. The evolution of the state of the condensate in a ballistic expansion is then presented, and the generation of Schrödinger cat-like states is observed. This is relevant in studies of output couplers or atom lasers. Finally we conclude in Section 4.

## 2 Q-Function in the Hartree approximation

The phase space representation is useful for visualising the evolution of a quantum state, as the statistics of the state can easily be described using the quasiprobability distribution functions . In this case, we have made use of the $`Q`$-function which has the nice property of being always positive. This function can be considered to describe probability densities, and has the further advantage of being readily measurable by quantum tomographic techniques . The first experimental measurement of quantum states was reported by Smithey et. al. in 1993 . Since then, several proposals to reconstruct the quasiprobability distributions of a quantum system, such as the Wigner function, have been made. For instance, Kurtsiefer et. al.  reported a technique to measure the Wigner function of an ensemble of helium atoms in a double-slit experiment. Moreover, a proposal has been made to measure the quantum state of a BEC in an atomic interferometer based on Raman transitions , where the $`Q`$-function can be measured directly. The way of obtaining the quasiprobability function is analogous to doing unbalanced homodyne detection of an optical signal.

The $`Q`$-function of a pure quantum state $`|\mathrm{\Psi }`$ is defined as

$$Q(\alpha _r,\alpha _i)\frac{1}{\pi }|\alpha |\mathrm{\Psi }|^2,$$
(1)
where the two real numbers $`\alpha _r`$ and $`\alpha _i`$ describe the coherent state $`|\alpha `$. Therefore, the quasiprobability distribution depends directly on these two parameters, which span the phase space. An alternative form of equation (1), leads the $`Q`$-function to be defined by the diagonal matrix elements of a density operator $`\widehat{\rho }`$ in a pure coherent state $`|\alpha `$:

$$Q(\alpha _r,\alpha _i)\frac{1}{\pi }\alpha |\widehat{\rho }|\alpha .$$
(2)
Clearly, the $`Q`$-function is characterised by being always positive and normalised to unity. Thus, it can be regarded as describing probability densities, and the numbers $`\alpha _r`$ and $`\alpha _i`$ play the role of conjugate variables, eg number and phase.

We shall start by describing the system in terms of the many-body interacting Hamiltonian

$`\widehat{H}`$ $`=`$ $`{\displaystyle \mathrm{d}^3𝐫\widehat{\varphi }^{}(𝐫,\tau )\left[\frac{\mathrm{}^2}{2m}_𝐫^2+V_{trap}(𝐫)\right]\widehat{\varphi }(𝐫,\tau )}`$ (4)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \mathrm{d}^3𝐫\mathrm{d}^3𝐫^{}\widehat{\varphi }^{}(𝐫,\tau )\widehat{\varphi }^{}(𝐫^{},\tau )\widehat{V}(𝐫𝐫^{})\widehat{\varphi }(𝐫,\tau )\widehat{\varphi }(𝐫^{},\tau )},`$
where $`\widehat{\varphi }^{}(𝐫,\tau )`$ and $`\widehat{\varphi }(𝐫,\tau )`$ are the boson field creation and annihilation operators respectively, and $`\widehat{V}(𝐫𝐫^{})`$ is the two-body interatomic potential. For a contact interaction potential, the normalised Heisenberg equation of motion for the operator $`\widehat{\varphi }`$ is given by

$$\mathrm{i}\frac{\widehat{\varphi }}{t}=\left[_𝐫^2+V_{trap}(𝐫)+\frac{nU_0}{\mathrm{}\omega _{trap}}\widehat{\varphi }^{}\widehat{\varphi }\right]\widehat{\varphi }.$$
(5)
In the Schrödinger picture, equation (5) takes the form $`\mathrm{i}\mathrm{}\frac{\mathrm{d}}{\mathrm{d}t}|\mathrm{\Phi }=\widehat{H}|\mathrm{\Phi }`$. The solution can then be found by expanding the state vector $`|\mathrm{\Phi }`$ in a Fock space

$$|\mathrm{\Phi }=\underset{n}{}A_n|n;t.$$
(6)
The solutions of the equation are then number states $`|n;t`$. In the Hartree approximation, it is possible to write the number states as

$$|n;t=\frac{1}{\sqrt{n!}}\left(dx\mathrm{\Psi }_n(x,t)\widehat{\varphi }^{}(x)\right)^n|0,$$
(7)
where $`\mathrm{\Psi }_n`$ is a single particle wavefunction and $`\widehat{\varphi }^{}`$ corresponds to the creation operator of the field.

This implies that each particle in the system experiences the same potential. Therefore, the single-particle function $`\mathrm{\Psi }_n`$ satisfies the GPE with a scaled nonlinearity proportional to the number of particles .

If we superimpose these states with a Poissonian distribution, the coefficient $`A_n`$ takes the form

$$A_n=\frac{\alpha _0^n}{\sqrt{n!}}\mathrm{exp}\frac{|\alpha _0|^2}{2}$$
(8)
so that the quantum state in the Hartree approximation is expressed as

$$|\mathrm{\Phi }=\underset{n}{}\frac{\alpha _0^n}{n!}\mathrm{exp}\frac{|\alpha _0|^2}{2}\left(dx\mathrm{\Psi }_n(x,t)\widehat{\varphi }^{}(x)\right)^n|0.$$
(9)
The $`Q`$-function given by equation (1) can be written in terms of a reference coherent state $`|\alpha ,\{\mathrm{\Psi }(x,t)\}`$ as

$$Q(\alpha _r,\alpha _i)=|\alpha ,\{\mathrm{\Psi }(x,t)\}|\mathrm{\Phi }|^2.$$
(10)
We can construct the reference coherent state $`|\alpha ,\{\mathrm{\Psi }(x,t)\}`$ using the $`n`$-particle eigenstate for the field with envelope $`\mathrm{\Psi }_{\overline{n}}(x,t)`$, where $`\overline{n}=|\alpha _0|^2`$ has the meaning of the average particle number

$$|n,\{\mathrm{\Psi }_{\overline{n}}(x,t)\}=\frac{1}{\sqrt{n!}}\left(dx\mathrm{\Psi }_{\overline{n}}(x,t)\widehat{\varphi }^{}(x)\right)^n|0.$$
(11)
Assuming Poissonian statistics, we obtain the many-body coherent state

$$|\alpha ,\{\mathrm{\Psi }_{\overline{n}}(x,t)\}=\mathrm{exp}\left(\frac{|\alpha |^2}{2}\right)\mathrm{exp}\left(\alpha dx\mathrm{\Psi }_{\overline{n}}(x,t)\widehat{\varphi }^{}(x)\right)|0.$$
(12)
Using equations (9), (10) and (12), the $`Q`$-function can be written

$$Q(\alpha _r,\alpha _i,t)=\mathrm{exp}\left(|\alpha _0|^2|\alpha |^2\right)\times \left|\underset{n=0}{\overset{\mathrm{}}{}}\frac{(\alpha ^{}\alpha _0)^n}{n!}\left(dx\mathrm{\Psi }_{\overline{n}}^{}(x,t)\mathrm{\Psi }_n(x,t)\right)^n\right|^2.$$
(13)
## 3 Quantum state of a Bose-Einstein Condensate at zero temperature

Let us consider Gaussian wavefunctions to describe the ground state of the trapped condensate

$$\mathrm{\Psi }_n=\mathrm{exp}\left[\left(\frac{x}{n1}\right)^2\right]$$
(14)
and assume that $`\overline{n}`$ is large so that the difference between $`\overline{n}1`$ and $`\overline{n}`$ can be neglected. The $`Q`$-function for such a case is

$$Q_G(\alpha _r,\alpha _i)=\mathrm{exp}\left(|\alpha _0|^2|\alpha |^2\right)\times \left|\underset{n=2}{\overset{\mathrm{}}{}}\frac{(\alpha ^{}\alpha _0)^n}{n!}\left[\sqrt{\pi }\left(\frac{\overline{n}^2n^2}{n^2+\overline{n}^2}\right)\right]^n\right|^2.$$
(15)
For an average number of particles $`\overline{n}=50`$ the $`Q`$-function (15) is shown in Figure 1.

A more realistic approximation to the ground state is given by the Thomas-Fermi solution

$$\mathrm{\Psi }_n=\frac{1}{\sqrt{n(n1)U_0}}\left(\mu _nV_{trap}\right)^{1/2},$$
(16)
where the chemical potential $`\mu _n`$ is determined by the normalisation of $`\mathrm{\Psi }_n`$. The $`Q`$-function is then expressed as

$`Q_{TF}(\alpha _r,\alpha _i)`$ $`=`$ $`\mathrm{exp}(|\alpha _0|^2|\alpha |^2)\times `$ (17)
$`|{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\alpha ^{}\alpha _0)^n}{n!}}\{{\displaystyle \frac{4}{3}}\mu _n[(\mu _n\mu _{\overline{n}})\mathrm{F}(\sqrt{{\displaystyle \frac{\mu _n}{\mu _{\overline{n}}}}},\sqrt{{\displaystyle \frac{\mu _{\overline{n}}}{\mu _n}}})`$
$`(\mu _n+\mu _{\overline{n}})\mathrm{E}(\sqrt{{\displaystyle \frac{\mu _n}{\mu _{\overline{n}}}}},\sqrt{{\displaystyle \frac{\mu _{\overline{n}}}{\mu _n}}})]\left\}^n\right|^2,`$
where $`\mathrm{F}(\phi ,k)`$ and $`\mathrm{E}(\phi ,k)`$ are incomplete elliptic integrals of the first and second kind respectively . For a harmonic trap potential $`V_1=\frac{1}{2}m\omega _{trap}^2r^2`$, the chemical potential is given by $`\mu _n=\left(\frac{3n(n1)U_0}{8}\right)^{2/3}`$. The $`Q`$-function for the same average number of particles as in the previous case is plotted in Figure 2.

From Figures 1 and 2, it is clear that the ground state of the condensate is very close to a squeezed state, since the $`Q`$-function has different widths in the $`\alpha _r`$ and $`\alpha _i`$-directions. On the one hand, the squeezing is generated by the nonlinearity inherent to the system. This can be explained in terms of a superposition of different numbers of particles with different phases, creating a deviation from the classical phase. On the other hand, we know that the coherent state interpretation arises from the robustness of these states in the presence of interactions between the condensate and its environment . It is well known that coherent states are a particular case of a more general class of minimum-uncertainty states, namely the squeezed states. It has been shown that the propagation of a coherent state in a nonlinear medium evolves into a superposition of two coherent states , while the quantum-mechanical superposition of two coherent states has been shown to lead to a transition from a Poisson distribution to a sub-Poissonian one and squeezing . In this sense, the robust description of the quantum state still applies, bearing in mind that the existence of the interactions between atoms yields squeezing of the quantum field and, as a result, entanglement between the atoms is also produced . We note that while a single coherent state does not evidence squeezing, a considerable amount is exhibited in a superposition of even just two coherent states.

This mechanism of achieving squeezing is the predominant one in the case of the ground state, without the presence of highly populated excited states. When we consider the presence of a thermal cloud, we can treat the system using the Bogoliubov procedure and analyse the behaviour of the collective excitations. It has been shown that squeezed quasiparticle excitations can be produced when considering the generation of correlated pairs by the interaction of an excited state with the ground state . Thus, we can talk about two different aspects of nonclassicality in the system, the predominant one for the ground state at zero temperature, and the prevalent one in the presence of excited states.

It is of particular interest to see how the $`Q`$-function behaves once the condensate is released from the trap as when the atomic cloud is allowed to expand ballistically , or in some output coupling schemes for atom lasers . The Thomas-Fermi solution was used as an initial condition to solve the GPE after the trapping potential is switched off and the $`Q`$-function (17) was then calculated from the dynamical evolution of the wavefunctions. We expect that once the trapping potential is switched off, the condensate will rethermalise. Comparing the quantum state of the condensate, for a trap frequency $`w_{trap}=2\pi \times 100`$ Hz, at $`2\times 10^3`$ s (shown in Figure 3) to the $`Q`$-function for a thermal state (Figure 4) makes it clear that rethermalisation is indeed taking place.

A very important aspect that has to be considered is the finite lifetime of the condensate, which is at most 10-20 s. The reason behind this, is that the thermal cloud is not the only source of decoherence. The condensate loses atoms because of Rayleigh scattering, external heating, and three-body decay. Information about the quantum state of the condensate is carried by the atoms that escape from it so that, due to interactions of the system with the environment, a quantum superposition is turned into a statistical mixture. From the processes mentioned above, three-body recombination is the most important one . In an experiment carried out by Stamper-Kurn et. al. , the measured loss rate per atom turned out to be $`4/s`$ for $`N=10^7`$ atoms. This rate scales as $`N^3`$, where $`N`$ is the number of atoms in the condensate. Considering a condensate with $`N=10^4`$, just one atom is lost per second; the decoherence time is 1 s. Bearing this fact in mind, it is perfectly sensible to ask whether, after opening the trap, it is still possible to see a Schrödinger cat-like state. When the states evolve further, a multi-component structure develops as shown in Figures 5 and 6, at $`5\times 10^3`$ s and $`1\times 10^2`$ s, respectively. In both figures, the quasiprobability function demonstrates that the quantum states obtained correspond to Schrödinger cat states like the ones generated in Kerr media . Moreover, we also see that the cat states have a different number of components at different times in the evolution, in a close analogy to the mechanism discussed by Yao for the case of a radiation field propagating in a nonlinear medium . It has been shown that the quantum-mechanical superposition of states can yield a state that exhibits squeezed fluctuations . This suggests that, even in the presence of losses due to three-body recombination, it would still be possible to observe Schrödinger cat-like states.

## 4 Discussion

In summary, we have described the ground state of a Bose condensed gas in the Hartree approximation in a phase-space representation using the $`Q`$-function. We have confirmed that the ground state of the condensate is in a squeezed state due to the presence of interactions between the atoms. We have also studied the dynamical evolution of the condensate in a ballistic expansion. The simulations show that immediately after the condensate is released from the trap the system rethermalises, with the quantum state being still close to a squeezed state. Continuing the evolution without the trapping potential, we obtain a multimode Schrödinger cat-like structure. In general, the “longevity” of these states is very small due to decoherence. However, the time scales of the calculations suggest that it might still be possible to observe them. The generation of nonclassical states by the time evolution of an initial coherent state under the influence of some nonlinear dynamics has been studied in the past for the case of radiation fields . In the case of a BEC the nonlinearity arises from binary correlations between the atoms, producing the effective interactions between them. We point out that this process basically refers to the ground state, as the presence of excited states changes the predominant mechanisms of achieving squeezing.

An issue that remains to be addressed is the one related to the experimental reconstruction of the quantum state of the condensate. We have mentioned that it is possible to use quantum tomographic techniques to measure the $`Q`$-function . As an alternative to these techniques, a scheme using projection synthesis was proposed by Baseia and co-workers . The scheme follows the idea of measuring the probability distribution of a state $`|f`$ to have an observable $`\theta `$

$$f|\theta \theta |f=P_N(\theta )=\frac{1}{2\pi }\left|\underset{n=0}{\overset{N}{}}c_n\mathrm{exp}(\mathrm{i}n\theta )\right|^2,$$
(18)
where the expectation value is replaced with the probability density by setting $`c_n`$ to zero for $`n>N`$, with large $`N`$. The probability density is proportional to the expectation value of the projector $`\widehat{\pi }=\kappa |\theta \theta |`$ where $`\kappa `$ is a positive constant. If the field is in a mixed state $`\widehat{\rho }`$, equation (18) can be written as

$$P_\rho (\theta )=\mathrm{Tr}(\widehat{\rho }\widehat{\pi })=\theta |\widehat{\rho }|\theta .$$
(19)
It has been shown  that by specifying a suitable reference state $`|B_b`$ in an experimental setup where a state $`|\psi _a`$ is coherently mixed in a beam splitter with the reference state $`|B_b`$, we can obtain the required probability (19).

In the case of the quasiprobability function $`Q`$, we would require a projector $`\widehat{\pi }=\kappa |\alpha \alpha |`$ so that its substitution in equation (19) yields

$$P_\rho (\alpha )=\mathrm{Tr}(\widehat{\rho }\widehat{\pi })=\alpha |\widehat{\rho }|\alpha =Q(\alpha _r,\alpha _i)$$
(20)
The form of the reference state that leads to the required projector is given by

$$|B_b=C\underset{k=0}{\overset{N}{}}b_k|k_b,$$
(21)
where $`C`$ is a normalisation constant and

$$b_k=\frac{2^N}{C}\left(\begin{array}{c}N\\ k\end{array}\right)^{1/2}\mathrm{exp}\left(\frac{\mathrm{i}k\pi }{2}\right)C_{Nk}^{}.$$
(22)
It is important to mention that this scheme involves projections onto truncated coherent states, which is a good approximation for large $`N`$. In other words, it works well for the case $`N\overline{n}`$, where $`\overline{n}`$ is the average particle number of the synthesised coherent state. The construction of the reference state (21) would be the main experimental challenge in this scheme. However, one might able to generate such a state using the scheme proposed in  where the authors show that it is possible to prepare an arbitrary (finite) superposition of Fock states by applying a well-defined succession of the displacement and creation operators. Experimentally a sequence of beam splitters at which conditional measurements are performed is required. In the case of BEC, this may be achieved by using repeated interference of condensates with different hyperfine states, similar to the Gedanken experiment proposed by Castin and Dalibard in reference  where atoms leak from two trapped condensates and are detected in the output channels of a 50-50 beam splitter. Exact procedure for such state preparation is the subject of our future work

This work has been supported by CONACyT, EPSRC, the Royal Commission for the Exhibition of 1851 and the EU.