Recent Progress in Gamow Shell Model Calculations of Drip Line Nuclei

The Gamow shell model (GSM) is a powerful method for the description of the exotic properties of drip line nuclei. Internucleon correlations are included via a configuration interaction framework. Continuum coupling is directly included at basis level by using the Berggren basis, in which, bound, resonance, and continuum single-particle states are treated on an equal footing in the complex momentum plane. Two different types of Gamow shell models have been developed: its first embodiment is that of the GSM defined with phenomenological nuclear interactions, whereas the GSM using realistic nuclear interactions, called the realistic Gamow shell model, was introduced later. The present review focuses on the recent applications of the GSM to drip line nuclei.


Introduction
Exotic nuclei have been studied for many years using a new generation of accelerators, which are now able to reach nuclear drip lines [1][2][3][4]. Contrary to well-bound nuclei, which are closed quantum systems, drip line nuclei can be seen as open quantum systems, as they are either weakly bound or unbound with respect to the particle emission threshold [5]. Many interesting phenomena appear at drip lines, such as a halo structure [1,6,7] and particle emission in resonance states [4,8]. Continuum coupling plays an important role in these loosely bound and unbound nuclear systems [5]. The proper description of nuclei at drip lines is one of the main challenges of nuclear theory, which was mostly developed to account for the structure of well-bound nuclei [5,6].
A clear consequence of the strong intertwinings of the continuum degrees of freedom and internucleon correlations at drip lines consists of the odd-even staggering found in the helium chain [9,10]. Indeed, odd helium isotopes (except 3 He) are all resonances and bear widths of several hundreds of keV [10][11][12]. Conversely, the even-even helium isotopes 4,6,8 He are bound, with 6,8 He both exhibiting halo properties [13][14][15]. To accurately reproduce nuclear halos, many-body wave functions in asymptotic regions must be treated properly, which demands to take into account continuum coupling [1,6,7,[16][17][18][19]. Adding to that, these weakly bound and unbound drip line nuclei also provide good laboratories to understand the single-particle structure, continuum coupling, internucleon correlations, and nucleon-nucleon (NN) interactions, which are not well understood in these regions.

Method
GSM is built within a configuration interaction framework based on the one-body Berggren basis [5,[36][37][38]. The Berggren basis [40,49] is generated by a finite-range potential, which can be written as the solutions of the one-body Schröndinger equation in the complex momentum space, which reads where l is the orbital angular momentum of the nucleon motion, m is the mass of the nucleon, r stands for the radius, andh is the reduced Planck constant. The momentum k and wave function u(k, r) can be complex. U(r) is the finite-range potential, which is, in practice, a Woods-Saxon (WS) [50] or GHF potential [33,44]. When considering protons, the Coulomb potential must be included in U(r). Bound, resonance, and scattering states can then be generated. The eigenenergy of single-particle states in the above equation is complex in general, and readsẽ n = k 2 /2m = e n − iγ n /2, where n denotes the state [40,49]. e n stands for the energy, whereas γ n represents the particle decay width, so that γ n = 0 for bound states and γ n > 0 for resonance states. A schematic Berggren basis set of states in the complex k-plane is illustrated in Figure 1. The wave function of a resonance state is not square-integrable, as its exponential increase in modulus implies that the wave function of a resonance state cannot be normalized with conventional techniques [40,49]. Consequently, one has to rely on the complex scaling method, which has been seen to properly account for the normalization of resonance states [51]. The completeness relation borne by Berggren basis states [40,49] reads where |n states are bound states and resonance states inside the L + contour of Figure 1. These states are called pole states, as they are the S-matrix poles of the finite-range potential. |k states are scattering states and follow the L + contours in the complex k-plane, starting from k = 0 and going to k → +∞, as shown in Figure 1. Scattering states initially form a continuum. Hence, in order to be used in numerical applications, the scattering states along the L + contour must be discretized with a Gauss-Legendre quadrature [5,49]. It has been checked that 10-45 states per contour are necessary to have converged results [5,38]. Once discretized, the Berggren basis is, in effect, the same as that of the harmonic-oscillator (HO) states within the standard SM [5,49]. Concerning resonance states, only narrow resonance states contribute to the physical states, and thus are included in the real calculations, whereas broad resonance states are not included, as they lie below the L + contour. In fact, the Berggren basis is the complex extension of the real-energy completeness relation of Newton [52], which consists of bound states and of a continuum of real-energy scattering states. Contrary to the Newton completeness relation [52], with which, only localized states can be expanded, the Berggren basis can expand unbound resonance states [40,49]. The many-body completeness relation is obtained by constructing Slater determinants from the one-body Berggren basis, which contains bound, resonance, and scattering states [5,49]. In the GSM, the Hamiltonian is represented by a complex symmetric matrix when using the one-body Berggren basis, which has to be diagonalized [5,49]. This process can be handled efficiently by using the complex symmetric extension of the Jacobi-Davidson method [49,53], where one can take advantage of the relatively small coupling to continuum states in order to have a fast convergence of calculations [5]. The full configuration space is extremely large due to the many scattering states within the model space. In practical calculations, however, we often truncate basis model spaces so that only two particles can occupy scattering states. It has been checked that this is sufficient to obtain converged results for both the energy and decay width of many-body states [5,38,54].
In GSM calculations, an effective Hamiltonian must be constructed. There are two main methods to build the effective Hamiltonian in GSM calculations. One is to construct an effective Hamiltonian based on realistic nuclear force [38,39,41], and hence in the frame of the realistic GSM, whereas the other one consists of using an effective phenomenological nuclear potential [5,36,37,49], in which, the parameters of the potential are optimized to reproduce experimental data. In the following, we give details about these two versions of GSM.

Realistic Gamow Shell Model Calculations
Within realistic GSM, one starts from the intrinsic Hamiltonian of an A-body system, which reads where p i is the nucleon momentum in laboratory frame, P = ∑ A i=1 p i is the center-ofmass (CoM) momentum of the system, and V (ij) NN is the realistic NN interaction, such as CD-Bonn [55] and N 3 LO [56] interaction. In the above Hamiltonian, the CoM energy is removed. In order to construct the effective Hamiltonian to be used in GSM calculations, an auxiliary potential is usually introduced [38,57,58], so that the Hamiltonian can be rewritten as, where has a one-body form, and Am ) is the residual two-body interaction, including the correction issued from the CoM motion. For the auxiliary potential U, we usually take a WS finite-range potential. To speed up the convergence of many-body calculations, the bare force is often softened by a similarity renormalization group method [59], such as the similarity renormalization group (SRG) and V low-k ,in order to remove the strong short-range repulsive core.
The realistic NN interaction is firstly defined in a relative momentum space. However, the many-body problem is usually solved in the laboratory frame (with, e.g., the HO basis), so that a transformation from relative and CoM coordinates to the laboratory frame is necessary. This procedure can be conveniently carried out in the HO basis via Brody-Moshinsky brackets [60]. In the HO basis, the two-body completeness relation reads where |αβ is the two-particle state of the HO basis. The two-body interaction in the HO basis is given by where N shell = 2n + l, indicates that a finite summation is performed up to N shell . The GSM calculations are carried out in the Berggren basis, so that the transformation of the interaction matrix elements from the HO basis to the Berggren basis needs to be carried out. This is achieved, in practice, by computing the overlaps between the Berggren and HO basis wave functions, where |ab (|cd ) is a two-particle state of the Berggren basis. For identical particles (protonproton or neutron-neutron), the overlap of the two-body state reads where J is the total angular momentum of the two-particle state, and j is the angular momentum of a single-particle basis state. The a|α ( b|β ) is the overlap of the one-body basis state, and δ αβ is the Kronecker delta. For the proton-neutron coupling, the overlap of the two-body state is simply given by The overlaps of one-body basis states are directly obtained from an integration in r-space where u(r) and R(r) are the radial parts of the single-particle Berggren and HO basis wave functions, with l, j, and t being the orbital, total angular momentum, and isospin quantum number, respectively. The single-particle wave functions of resonance and continuum states are not localized and hence are not square-integrable. The transformation defined by Equation (7), in fact, utilizes the short-range nature of nuclear force. Indeed, the Gaussian fall-off of the HO wave function renders the overlaps integrable, even when one considers resonances or scattering states of complex energy. For long-range operators, such as the one-body kinetic energy and Coulomb potential, using Equation (7) is not suitable in practice. In this case, we use the exterior complex scaling technique [51] to treat the kinetic and Coulomb operator, i.e., terms proportional to p 2 and 1/r, respectively, with the Berggren basis. The obtained interaction matrix elements in the Berggren basis are complex, and associated operators are non-Hermitian. The many-body perturbation theory (MBPT), named the fullQ-box folded-diagram method [61], is employed to construct the realistic complex effective Hamiltonian in the defined model space for GSM calculations. The complex-k Berggren basis states in the model space are non-degenerate; therefore a non-degenerateQ-box folded-diagram perturbation, i.e., the extended Kuo-Krenciglowa (EKK) method [62], has been used. For this, we first calculate theQ-box using MBPT in the Berggren complex-k basis, where E is starting energy, P and Q represent the model space and the excluded space, respectively, with P + Q = 1. TheQ-box is composed of irreducible valence-linked diagrams [57,58], which can be built order-by-order. V and H are the two-body interaction and two-body Hamiltonian, respectively, and H 0 is the unperturbed one-body Hamiltonian. The derivatives of theQ-box are defined as where s denotes the s-th derivative.
The effective Hamiltonian H eff can then be constructed in operator form [63], written as where H eff stands for Hamiltonian shifted by an energy E, with H eff is obtained by performing iterations of Equation (13), which is equivalent to calculate folded-diagrams, where one considers high-order contributions by summing up the subsets of diagrams to finite order. When convergence is obtained, the effective Hamiltonian is given by H eff = H eff − E, and the effective interaction reads V eff = H eff − PH 0 P. The extendedQ-box folded-diagram calculations provide a useful approach for including effects from the continuum coupling and core polarization [57,58,61,62].

Gamow Shell Model with Phenomenological Nuclear Potential
Within the Gamow shell model with phenomenological nuclear potential, the nucleus is assumed to be a system of N v valence particles outside a frozen closed core, from which, core polarization is absent [5,36,37,49]. The GSM Hamiltonian, expressed with intrinsic nucleon-core cluster-orbital shell model coordinates [64], writes where p i is the nucleon momentum in cluster-orbital shell model frame, U core is the singleparticle nucleon-core potential, and V is the phenomenological NN interaction between valence nucleons. µ i and M core stand for the reduced mass of the nucleon and the mass of the core, respectively. The p i p j M core term accounts for the two-body recoil term. As seen in Equation (15), the GSM has two components: the one-body part The core-valence particle potential U core is usually a WS potential, in which a spin-orbit term is included. The NN interaction V takes the form of an effective phenomenological NN interaction, such as Minnesota [65], Furutani-Horiuchi-Tamagaki (FHT) [66,67], or effective field theory (EFT) [18,56] interactions. The parameters within the effective Hamiltonian in Equation (15), both one-and two-body interactions, need to be optimized to reproduce experimental data. For optimizations, the χ 2 optimization method is employed, where one uses the Gauss-Newton algorithm augmented by the singular value decomposition technique to calculate the Jacobian pseudo-inverse [43].

Neutron-Rich Oxygen and Fluorine Isotopes
Neutron-rich oxygen isotopes form a particularly interesting chain for experimental and theoretical research. Firstly, the proton number Z = 8 shows magical properties for the neutron-rich oxygen isotopes, which provide a good laboratory to perform configuration interaction (shell-model) calculations [22,38,68,69]. Secondly, the nuclei 22 O and 24 O exhibit doubly magicity at the neutron number N = 14 and 16, respectively, [70][71][72][73]. Thirdly, experiments have shown that the 25 O and 26 O are unbound and decay by one-and twoneutron emission, respectively, [8,74]. Experimental studies suggest that 24 O is the heaviest bound isotope of the oxygen chain [8,74]. However, the loosely unbound property of 26 O, which is only −18 keV unbound [8], is a strong incentive to investigate the bound or unbound character of 28 O, which should have a magicity of N = 20. Consequently, the neutron-rich oxygen isotopes provide an ideal laboratory to study many-body correlation, continuum coupling, and single-particle structure. By adding one valence proton to the neutron-rich oxygen isotopes, one obtains the fluorine isotopes at the neutron drip line, which can sustain six additional neutrons after 25 F, hence, up to 31 F, which is suspected to be at the neutron drip line of the fluorine chain [75]. This dramatic change is called an "oxygen anomaly". Moreover, many exotic properties develop at the neutron drip line for fluorine isotopes, such as halos in 29 F [76] and 28 F within the island of inversion around N = 20 [77], and thus fluorine isotopes provide a very interesting ground for theoretical studies.

Realistic Gamow Shell Model Calculations
We have developed realistic GSM with the Berggren basis using a WS potential, while the realistic effective Hamiltonian is constructed within the model space using a nondegenerateQ-box folded-diagram method [38]. We first employed it to investigate the neutron-rich oxygen isotopes up to and beyond the neutron drip line [38]. In our calculations, the realistic CD-Bonn potential [55] was used. To speed up the convergence of many-body calculations, the bare force is usually softened to remove the strong short-range repulsive core. The V low−k method [59] is used for that matter in Ref. [38].  24,25,26 O, compared with available experimental data [8,10,74]. The resonances are indicated by shades, and their widths (in MeV) are given by the number below or above the levels. The light blue shade indicates the 3/2 + many-body scattering states (with permissions from Ref. [38]). The "CGSM" stays for "core Gamow shell model". Figure 2 shows the calculated low-lying states of 24−26 O, along with experimental data [8,10,74]. Our realistic GSM calculations [38] reproduce the experimental excitedstate spectrum well, including the observed resonance widths. The ground state energies and one-neutron separation energies S n of the neutron-rich oxygen isotopes are also calculated [38] (see Figure 3) and compared to the experimental data [8,10,74,78]. The WS parameters used, taken from Ref. [38], reproduce the experimental 1s 1/2 and 0d 3/2 single-particle energies well, including the decay width of the 0d 3/2 state, but give the 0d 5/2 energy as lower than the experimental data, at about 1.17 MeV [10]. The results presented in Figure 3 show that adopting the experimental 0d 5/2 energy can dramatically improve calculations. Overbinding in the GSM calculations of oxygen isotopes after 24 O is obtained in Ref. [38], which is caused by the absence of the three-nucleon force (3NF).  [8,74,78]. "GSM with WS SPE" indicates that the calculations were performed with Woods-Saxon (WS) single-particle energies (SPE), whereas "GSM with optimized SPE" means that the calculations were performed with the 0d 5/2 SPE replaced by its experimental value (with permissions from Ref. [38]).

Ab-initio Realistic GSM Calculations within GHF Basis
GSM is usually performed using a basis generated by a WS potential [5,19,[36][37][38][39], whose parameters must be determined by fitting experimental single-particle energies and resonance widths. However, the single-particle energies and resonance widths in the multishell case are sometimes difficult to assess due to the lack of experimental data for that matter [10]. We then developed an ab initio realistic GSM approach by introducing the GHF basis as the Berggren basis [41]. The GHF basis is obtained by using the same interaction as the one used in the construction of the effective SM Hamiltonian [41], and thus there is no parameter introduced in the GHF Berggren basis. Starting from the chiral next-to-next-toleading-order (NNLO opt ) force [79], we perform a nondegenerateQ-box folded-diagram calculation [38,62] in the GHF basis in order to construct a complex effective Hamiltonian. The energies and widths of single-particle orbitals can also be obtained self-consistently using the nondegenerateŜ-box folded-diagram method [41]. The neutron-rich fluorine isotopes have been extended to the p f -shell, using a cross-shell effective Hamiltonian with the following model space : {1s 1/2 , 0d 5/2 , 0d 3/2 } for the valence proton, and {1s 1/2 , 0d 3/2 + d 3/2 scattering states, 1p 3/2 + p 3/2 scattering states, 1p 1/2 + p 1/2 scattering states, f 7/2 scattering states} for valence neutrons. More details can be found in Ref. [41]. The constructed effective Hamiltonian was employed to study neutron-rich oxygen and fluorine drip line nuclei. Figure 4 shows the calculated ground-state energies and neutron separation energies S n of oxygen and fluorine isotopic chains, as well as comparisons with experimental data [78] and other theoretical calculations [31,68,[79][80][81][82][83]. The GSM calculations using a GHF basis and based on the NNLO opt [79] provide the correct location of the neutron drip line of oxygen isotopes and a good description of the unbound nuclei 25,26 O, which lie beyond the neutron drip line (see the left panel of Figure 4). Note that, when using the standard SM calculations with the USDB interaction [68], conventional SM calculations with NN + 3NF [82], or valence-space IMSRG (VS-IMSRG) calculations with NN + 3NF [81], the resonance and continuum couplings are absent. Complex CC [31] and GSM calculations [80] are displayed in Figure 4 for comparison.  [78] (the AME2016 extrapolated values are taken for 27,28 O and 30,31 F) and theoretical calculations from other groups: complex coupled-cluster (CC) with next-to-next-to-next-to-leading-order nucleon-nucleon CC with NNLO opt interaction [79], GSM [80], -space in-medium similarity renormalization group (VS-IMSRG) [81], SM with NN+3NF [82], SM with USDB [68], and SM with SDPF-M [83] (with permissions from Ref. [41]).
The results of fluorine isotopes are shown in the right panels of Figure 4. For comparison, standard SM calculations using USDB [68] and SDPF-M [83] effective interactions are also presented. All ground-state energies in Figure 4 are given with respect to the ground state of 22 O. Experiments revealed that 31 F is a neutron drip line nucleus [75]. Although our GSM calculations provide a lower energy of 31 F compared to that of 30 F, 31 F is still unbound compared to 29 F. However, our GSM calculations provide good descriptions of ground-state energies for 23−29 F.

GSM Calculations with Phenomenological Nuclear Potential
Many ab initio calculations, such as SM [82], VS-IMSRG [81], complex CC [31], and realistic GSM [38,39,41] calculations, have been employed for the description of neutronrich oxygen and fluorine isotopes. However, these calculations bear a large theoretical uncertainty. Furthermore, results arising from ab initio calculations depend on the realistic nuclear forces used (a short summary of the VS-IMSRG calculations based on different chiral nuclear forces can be found in Ref. [80]). Moreover, continuum coupling is absent in the VS-IMSRG [81] and SM [82] calculations. Similar situations also occur for neutron-rich fluorine isotopes, where few calculations have been performed and most of the calculations are absent for the continuum coupling [83,84]. Based on these grounds, we performed the GSM calculations with a phenomenological nuclear interaction for neutron-rich oxygen and fluorine drip line nuclei.   [42]). Results are compared with the experimental data available, represented by a star. The data for 25,26 O and 27,28 O are taken from experiment (see Refs. [8,74]) and evaluations given in AME2016 [78] (with permissions from Ref. [42]).
For the considered neutron-rich oxygen isotopes, the closed-shell nucleus 22 O is selected to be the inner core. The one-body potential is mimicked by a WS potential, whose parameters are adjusted to reproduce the single-particle spectrum of 23 O [10]. We use the pionless EFT interaction [85,86] as the two-body interaction. Owing to the few available data related to the oxygen drip line nuclei [10], only the leading order (LO) NN interaction of the EFT force is fitted to reproduce selected experimental data. The effect of the 3NF at LO is then effectively taken into account in the fitted parameters. Details can be found in Ref. [42]. We calculated the energies of the ground states of 24−28 O with GSM within sdp f (s 1/2 , p 1/2,3/2 , d 3/2,5/2 , f 5/2,7/2 partial waves) and sd (s 1/2 , d 3/2,5/2 partial waves) active spaces, using different EFT interactions (see details in Ref. [42]). The calculated groundstate energies of 24−28 O are shown in Figure 5. The calculations within the sd space show that the 25−28 O isotopes are unbound and that their binding energies are close to the experimental data [8,74,78] and to calculations performed within the sdp f space. However, the calculations obtained in the sd space provide an unbound 26 O ground state, by about 300 keV relative to the ground state of 24 O, which is a little higher than its experimental value, which is about 20 keV unbound [8]. Though the energy difference obtained using the two different model spaces is small, the calculation performed within the spd f space seems to be more reasonable. The GSM calculations performed within the sdp f space provide good agreements of the 23−26 O ground states with experimental data [8,74,78]; in particular, the two-neutron separation energy (S 2n ) of 26 O is about 20 keV [8]. The calculated ground state of 28 O is unbound in all three cases and located about 700 keV above the ground state of 24 Figure 6. Energies of 25−31 F with respect to the 24 O core, calculated within different theoretical frameworks and compared to experimental data [78]. Besides the GSM calculations using FHT and EFT interactions, calculations in the Hartree-Fock many-body perturbation theory (HF-MBPT) [88] and VS-IMSRG [84] frameworks, utilizing the harmonic-oscillator (HO) basis, hence being without continuum coupling, are presented (with permissions from Ref. [18]). Figure 6 shows the GSM calculations of the binding energies of fluorine isotopes using FHT and EFT interactions (see details in Ref. [18]), wich are compared with experimental data [78] and Hartree-Fock MBPT (HF-MBPT) [88] and VS-IMSRG [84] calculations, which are both performed in the HO basis. The energy of 25 F has been fixed to its experimental datum in all used models in Figure 6. We can see that all calculations reproduce the ground state energies of 25−28 F isotopes well, situated in the well-bound region, whereas differences start after 29 F, i.e., when one reaches the neutron drip line. Due to the absence of both multi-shell and continuum couplings, the VS-IMSRG calculations [84] provide visible differences, which are about 4-to 5-MeV in magnitude for 30,31 F. When applying the HF-MBPT method [88], the cross-shell couplings generated by the sd and p f shells are included, so that proper binding energies of up to 29 F are predicted. However, due to the lack of continuum coupling, the binding energies of 30,31 F are about 1 MeV away from experimental data. The GSM calculation performed with FHT and EFT interactions correctly provides binding energies of up to 31 F. Moreover, the odd-even staggering encountered from the 28 F isotope, typical of the presence of a strong proton-neutron interaction, is well reproduced, with 30 F being unbound and 31 F being loosely bound in our calculations.
Recent realistic shell model calculations [89] have pointed out that nuclear deformation plays an important role in the neutron drip line nuclei. Within GSM, deformation can be accounted for by configuration mixing using a cross-shell valence space [5]. Besides deformation, continuum coupling also gives important contributions in drip line nuclei. They are strongly coupled with continuum states near the particle-emission threshold, which provides additional binding energy [5]. This situation is unlike that occurring in well-bound systems, where one only has strong coupling with nearby deeply-bound single-particle states. ) r ( f m ) 3 1 F 2 9 F 2 7 F 3 1 F 2 9 F Figure 7. One-nucleon densities of the bound 27,29,31 F isotopes calculated with the GSM using the EFT interaction in the valence space as a function of radii r, respectively, depicted by short-dashed, long-dashed, and solid lines. The rms radii of these isotopes are shown in the inset (with permissions from Ref. [18]).
The two-neutron separation energy S 2n of 31 F is about 170 keV [78], which is sufficiently small for sustaining a two-neutron halo. We calculated the one-nucleon densities and neutron rms radii of the neutron-bound 27,29,31 F isotopes using GSM with the EFT interaction (see Figure 7). From our calculations, a halo clearly develops in the asymptotic region of 31 F. Indeed, on the one hand, the one-nucleon density of 31 F slowly decreases on the real axis and is about one to two orders of magnitude larger than those of 27,29 F in the asymptotic region. Added to that, on the other hand, the neutron rms radius of 31 F does not follow the trend noticed in 27,29 F, as the rms radius sharply increases compared to the associated values in 27 F and 29 F. Consequently, one can assume from these GSM calculations [18] that 31 F is a two-neutron halo state.

Realistic Gamow Shell Model Calculations of Neutron-Rich Calcium Isotopes
The long chain of calcium isotopes provides an ideal laboratory for both theoretical and experimental investigations of unstable isotopes. With two typical doubly magic isotopes, 40 Ca and 48 Ca, and two new magic isotopes discovered in the neutron-rich region, 52 Ca [90] and 54 Ca [91], the calcium chain is speculated to end the 70 Ca isotope. Its rich nuclear structure data [10] attract continued theoretical interest, especially using methods that include continuum coupling. The realistic GSM based on the realistic CD-Bonn [55] interaction has also been performed to investigate the properties of neutron-rich calcium isotopes up to the drip line.  Figure 8. Calculated one-neutron separation energies S n (a) and two-neutron separation energies S 2n (b), compared with data [78,92], and calculations obtained with SV-min density-functional theory (DFT) [93] and multireference IM-SRG (S 2n only) [94]. The S n calculations end at 60 Ca because odd isotopes heavier than 60 Ca become unbound in our GSM calculations (with permissions from Ref. [47]).
The calculated one-neutron separation energies S n and two-neutron separation energies S 2n are shown in Figure 8 and compared with experimental data [78,92], DFT [93], and IM-SRG [94] calculations. The calculated one-neutron separation energies S n show that 57 Ca is the heaviest odd-mass bound calcium isotope, which is consistent with MBPT calculations [95]. 59 Ca is weakly unbound with a small one-neutron separation energy S n = −326 keV in our GSM calculations. For the two-neutron separation energy S 2n , the GSM calculations are performed with two different cores, 48 Ca and 54 Ca. For 56,58,60 Ca, the two calculations give similar results. The calculated two-neutron separation energy S 2n is in good agreement with experimental data [78,92] and other theoretical calculations, e.g., with DFT [93] and IM-SRG [94] calculations. The large decrease in S 2n at neutron number N = 32 and 34 indicate that subshell closures occur therein, which has also been suggested from experiments [90,91] and theoretical calculations [94][95][96]. Moreover, the calculated two-neutron separation energy S 2n using GSM predicts that the two-neutron drip line of the calcium isotopes should be located at 70 Ca. This is consistent with the recent mean-field calculations of Ref. [97]. . Neutron effective single-particle energies (ESPE) with respect to the 48 Ca core, as a function of neutron number. The V low−k Λ = 2.6 fm −1 CD-Bonn interaction is utilized (with permissions from Ref. [47]).
In order to see the shell evolution of the calcium isotopes around the neutron number N = 32, 34, 40, and 50, we calculated effective single-particle energies (ESPE) based on the GSM effective Hamiltonian. Figure 9 shows the evolution of the valence neutron ESPEs when increasing the neutron number. The calculations show that large shell gaps between 1p 3/2 and 1p 1/2 and between 1p 1/2 and 0 f 5/2 exist, indicating that shell closures occur at N = 32 and 34, respectively. These results are consistent with experimental observations [90,91] and theoretical calculations [94][95][96]98]. The shell gap above the 0 f 5/2 orbit is reduced at around N = 40, implying a weakening of the N = 40 shell closure in the calcium chain. The 0g 9/2 shell becomes bound at N ≥ 40, which can enhance the stability of the heavy calcium isotopes. The observed 60 Ca isotope in experiments [99] may be an indication of this enhanced stability. Moreover, the calculated ESPEs show a significant shell gap at N = 50, implying a shell closure at 70 Ca.

One-Proton and Two-Proton Decays in 16 Ne and 18 Mg Unbound Nuclei
Two-proton decay is one of the most important drip-line phenomena. It occurs in proton drip line nuclei, such as 48 Ni, 54 F, 54 Zn, 76 K, 16 Ne, and 19 Mg (see a review of this topic in Ref. [4] ). While 18 Mg has not been observed, it can decay in principle by proton and/or two-proton emissions. The GSM is then a suitable method to study these types of particle emissions. We carried out GSM calculations of the proton-rich carbon isotones of 14 O, which are all resonance [10], using 14 O as an inner core. The obtained energy spectra of carbon isotones are depicted in Figure 10 with respect to the ground state of 14 O. One can see that both the energies and widths of experimentally known eigenstates are well reproduced for the low-lying states in 15 F and 16 Ne [10]. We also provide predictions for the 17 Na and 18 Mg nuclear spectra, of which, there are no experimental data. Our calculations show that the 16 Ne and 18 Mg isotopes are unbound nuclei, where both oneproton separation energies S p and two-proton separation energies S 2p are negative, thereby indicating that two different particle-emission channels are open therein. To evaluate one-proton and two-proton decay widths, we changed the central potential depth V 0 of the WS core potential in order for the S p to become positive or very negative (see details in Ref. [48]). Consequently, it is possible to find a central potential depth for which only the two-proton decay channel is open, so that the obtained width is that of the two-proton emission. The obtained results are shown in Figure 11. As 15 F and 17 Ne are one-proton resonances, their width increases steadily with the Hamiltonian central potential depth. In contrast, one can see that the widths of 16 Ne and 18 Mg increase abruptly when the one-proton channel opens. The width of two-proton decay is almost constant with respect to the central potential depth below the one-proton emission threshold, and is also about 500 keV to 1 MeV above (see Figure 11). It is reasonable to assume that the two-proton decay width is almost independent of energy. Therefore, the GSM results shown in Figure 11, where only the two-proton channel is open, can be extrapolated to the physical case (indicated by an arrow in Figure 11). This two-proton decay width is about 10-15 keV for both 16 Ne and 18 Mg nuclei. The one-proton width can be assumed as the difference between the total width and the two-proton emission width of 10-15 keV. Then, our calculations show that one-proton emission is negligible for 16 Ne, whereas the one-proton decay width in 18 Mg is estimated to be about 85-90 keV. The obtained data for 16 Ne are also in agreement with experimental data [10, [100][101][102].  17 Na, and 18 Mg (lower panel) as a function of the difference ∆V 0 = V 0 − V (fit) 0 (fit) of the WS central potential depths (see details in Ref. [48]). Energies are depicted by blue disks and red lozenges for even and odd nuclei, respectively. Widths are represented by segments centered on disks and lozenges. The widths of 16 Ne and 18 Mg have been multiplied by 20 for readability. Energies are given with respect to the 14 O core. The physical GSM calculation, for which V 0 = V (fit) 0 , is indicated by an arrow (with permissions from Ref. [48]).

Summary
The Gamow shell model (GSM) is a powerful method for the description of the weakly bound and resonance properties of drip line nuclei. In the present review, we presented several recent applications of GSM dedicated to the study of drip line nuclei, including GSM calculations of neutron-rich oxygen and fluorine drip line nuclei, of the long chain of neutron-rich calcium isotopes, and of the unbound proton-rich 16 Ne and 18 Mg isotopes. For the neutron-rich oxygen and fluorine drip line nuclei, both the realistic GSM and GSM with phenomenological forces have been utilized. Our calculations have described the weakly-bound and unbound properties of drip line nuclei well. Furthermore, the unbound properties of the 28 O are obtained within the two both types of GSM calculations, and the two-neutron halo property of 31 F has been predicted in GSM calculations as well. The realistic GSM calculations provide good agreements of the neutron-rich calcium isotopes with experimental data, as GSM calculations predict that the one-and two-neutron drip line nuclei of calcium isotopes are 57 Ca and 70 Ca, respectively. For the unbound proton-rich 16 Ne and 18 Mg nuclei, GSM calculations provide calculations and predictions for their lowlying spectra. Added to that, the one-and two-proton emission widths could be estimated for 16 Ne and 18 Mg isotopes. Our calculations have shown that 16 Ne decays only by twoproton emission, whereas 18 Mg can decay through both one-and two-proton emission channels, whose widths are estimated to be about 85-90 and 10-15 keV, respectively.
GSM has thus been shown to be the tool of choice for the study of drip line nuclei. Many challenges remain to be overcome for the future applications of GSM: • Due to the large computational cost of the GSM calculations, the GSM has been applied for the neutron-rich nuclei with only one or two valence protons, and, for proton-rich nuclei, with only one or two valence neutrons in the non-resonant continuum. For example, the model space dimension of 31 F is about 10 7 with two valence particles in the continuum. It can easily reach 10 10 without truncations, which is untractable numerically. In the nuclear chart, most of the drip line nuclei need to be described with many valence particles (protons and neutrons). One can think of the neutron-rich Ne and Mg drip line isotopes, where both continuum coupling and strong internucleon correlations must be treated properly. These isotopes will provide a challenge for future GSM calculations, due to the large dimensions; • The diagonalization of the GSM Hamiltonian in order to obtain eigenstates of large resonance widths, such as the second 0 + 2 state in 8 He, is very difficult from a numerical point of view; • The dimensions of the GSM Hamiltonian matrices increase extremely quickly when one adds valence particles, and thus the treatment of the many-body Hamiltonian is difficult when using the configuration interaction framework. Other kinds of manybody methods are urgently needed. The two-particle reduced density matrix method is one of the promising methods to solve the dimensionality problem of the GSM many-body Hamiltonian [103]; • The unbound single-particle states of s waves in neutron-rich nuclei are anti-bound states, which are difficult to include in many-body GSM calculations. The consideration of many-body anti-bound states in GSM (the ground state of 10 Li is supposed to be anti-bound, for example [104]) is thus also a challenge for future applications of GSM. Data Availability Statement: All of the relevant data are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.