The in-medium similarity renormalization group

The SRG is to evolve the Hamiltonian H(s) to be band-diagonal by using the continuous unitary transformation as [35, 34] $H\left(s\right)=U\left(s\right)H{U}^{†}\left(s\right)\equiv H^{\text{d}}\left(s\right)+H^{\text{od}}\left(s\right)\to H^{\text{d}}\left(\infty \right)\mathrm{,}$ (1) where s denotes the so-called flow parameter, and H^d(s) and H^od(s) are appropriately defined as the “diagonal” and “off-diagonal” parts of the Hamiltonian, respectively. Although the evolution should continue up to s→∞, a finite number of evolution steps is usually sufficient to make H(s) approach the band-diagonal form of H^d(∞).

Equation (1) expresses a general ideal. In practice, taking the derivative of Eq. (1), a flow equation is defined to evolve the Hamiltonian H(0), $\frac{\text{d}}{\text{d}s}H(s)=[\eta (s)\mathrm{,}H(s)]\mathrm{,}$ (2) where the anti-Hermitian generator $\eta (s)$ is related to the unitary transformation U(s) by $\eta (s)=\frac{\text{d}U(s)}{ds}{U}^{†}\left(s\right)=-{\eta }^{†}\left(s\right)\mathrm{.}$ (3) A commonly used generator is defined as $\eta (s)=\left[H^{\text{d}}\left(s\right)\mathrm{,}H\left(s\right)\right]=\left[H^{\text{d}}\left(s\right)\mathrm{,}{H}^{\text{od}}\left(s\right)\right]\mathrm{,}$ (4) which guarantees that the off-diagonal coupling of $H^{\text{od}}$ is driven exponentially to zero with increasing in the value of the flow parameter s [35]. In practice, the demand for strict diagonality is usually relaxed to band diagonality of the Hamiltonian matrix in a chosen basis, such as in relative momentum or harmonic oscillator (HO) spaces. In nuclear physics, the SRG is used to decouple the momentum or energy scales in free space to construct “soft” NN and 3N interactions, thereby rendering the nuclear Hamiltonian more suitable for ab initio many-body calculations [31, 33, 72-74].

The SRG is used to soften the nuclear force which has a hard core in free space. This renormalization can significantly accelerate ab initio calculations of nuclei. Another development of the SRG theory is the in-medium SRG (IMSRG) [36-38] which evolves the many-body Hamiltonian to block diagonal form. The decoupling between the lowest-energy ground state and excited states of the Hamiltonian directly provides the energy of the ground state of the nucleus. A distinct advantage of IMSRG, compared to the SRG free-space evolution, is its ability to approximately evolve 3, …, A-body operators using only two-body machinery. This simplification is primarily achieved through the use of normal ordering with respect to a reference state $|\text{Φ}\rangle$, usually the Hartree-Fock (HF) state.

The intrinsic Hamiltonian of A-body nuclear system is expressed as $H={\displaystyle \sum _{i=1}^A}\left(1-\frac{1}{A}\right)\frac{p_i^2}{2m}+{\displaystyle \sum _{i<j}^A}\left(v_{ij}^{N\text{}N\text{}}-\frac{p_i\cdot p_j}{mA}\right)+{\displaystyle \sum _{i<j<k}^A}{v}_{ijk}^{3N}\mathrm{,}$ (5) where p_i is the nucleon momentum in laboratory coordinates, and m is the nucleon mass, with $v^{N\text{}N\text{}}$ and $v^{3N}$ denoting the NN and 3N interactions, respectively. In order to generate the reference state in the IMSRG calculation, the HF equation for the intrinsic Hamiltonian Eq. (5) is first solved. The Wick’s theorem is applied to normal order all operators starting from a general second-quantized Hamiltonian with two- and three-body interactions, with respect to the HF ground state. $\begin{array}{l}H=E_0+{\displaystyle \sum _{ij}{f}_{ij}}:a_i^{†}{a}_j:+\frac{1}{{2!}^2}{\displaystyle \sum _{ijkl}{\text{Γ}}_{ijkl}}:a_i^{†}{a}_j^{†}{a}_l{a}_k:\\ +\frac{1}{{3!}^2}{\displaystyle \sum _{ijklmn}{W}_{ijklmn}}:a_i^{†}{a}_j^{†}{a}_k^{†}{a}_n{a}_m{a}_l:\mathrm{,}\end{array}$ (6) where E₀, f, Γ, and W correspond to the normal-ordered zero-, one-, two-, and three-body terms, respectively, given by $\begin{array}{c}{E}_0={\displaystyle \sum _i{T}_{ii}{n}_i}+\frac{1}{2}{\displaystyle \sum _{ij}\left(T_{ijij}+v_{ijij}^{N\text{}N\text{}}\right)}{n}_i{n}_j\\ +\frac{1}{6}{\displaystyle \sum _{ijk}{v}_{ijkijk}^{3N}}{n}_i{n}_j{n}_k\mathrm{,}\end{array}$ (7) $f_{ij}=T_{ij}+{\displaystyle \sum _k\left(T_{ikjk}+v_{ikjk}^{NN}\right)}{n}_k+\frac{1}{2}{\displaystyle \sum _{kl}{v}_{ikljkl}^{3N}}{n}_k{n}_l\mathrm{,}$ (8) ${\text{Γ}}_{ijkl}=T_{ijkl}+v_{ijkl}^{N\text{}N\text{}}+\frac{1}{4}{\displaystyle \sum _m{v}_{ijmklm}^{3N}}{n}_m\mathrm{,}$ (9) $W_{ijklmn}=v_{ijklmn}^{3N}\mathrm{,}$ (10) where $n_i=\theta \left({\epsilon }_{\text{F}}-{\epsilon }_i\right)$ represents the occupation numbers in the reference state $|\text{Φ}\rangle$, with ${\epsilon }_{\text{F}}$ denoting the Fermi energy of the reference state and T representing the kinetic part of the Hamiltonian. From Eqs. (7)-(10), it is evident that all the zero-, one- and two-body parts of the Hamiltonian contain the in-medium effects from the free-space 3N interactions.

The exact treatment of the 3NF is computationally expensive. Therefore, the residual 3NF is usually neglected, which provides a reasonably good approximation in nuclear structure calculations. The omission of the residual normal-ordered three-body component of the Hamiltonian has been shown to result in only 1%-2% discrepancy in ground-state and excited-states energies for light and medium-mass nuclei [16, 75]. The approximation of normal-ordered two-body (NO2B) for the Hamiltonian has been proved to be useful and beneficial in practical calculations, offering an efficient method to include 3NF effects in nuclear many-body calculations, thereby avoiding the computational burden of directly dealing with three-body operators.

Similar to the evolution of the Hamiltonian, the operators of other observables can also be evolved using the flow equation $\frac{\text{d}}{\text{d}s}O(s)=[\eta (s)\mathrm{,}O(s)]\mathrm{.}$ (11) The Magnus expansion was usually used in matrix differential equations [76], and was applied to reformulate the IMSRG [77] for more efficient calculations. In the Magnus approach, the IMSRG transformation can be written as an exponential expression [76], $U(s)\equiv e^{\text{Ω}\left(s\right)}\mathrm{.}$ (12) The Magnus evolution operator Ω(s) works for both the Hamiltonian and other observable operators, which allows the derivation of the flow equation for the anti-Hermitian Magnus operator Ω(s), $\frac{\text{d}\text{Ω}}{\text{d}s}={\displaystyle \sum _{k=0}^{\infty }}\frac{B_k}{k!}{\text{ad}}_{\text{Ω}}^k\left(\eta \right)\mathrm{,}$ (13) where Bk denote the Bernoulli numbers and ${\text{ad}}_{\text{Ω}}^0\left(\eta \right)=\eta$ (14) ${\text{ad}}_{\text{Ω}}^k\left(\eta \right)=\left[\text{Ω}\mathrm{,}{\text{ad}}_{\text{Ω}}^{k-1}\left(\eta \right)\right]\mathrm{.}$ (15) In practical calculations, η and Ω are truncated along with their commutators at the two-body level, called the Magnus(2) approximation. The series of nested commutators generated by ${\text{ad}}_{\text{Ω}}^k$ are recursively evaluated until a satisfactory convergence of the right-hand side of Eq. (13) is reached [77]. At each integration step, U(s) is used to construct the Hamiltonian H(s) via the Baker-Campbell-Hausdorff (BCH) formula $H(s)\equiv e^{\text{Ω}}\left(s\right)H(0)e^{-\text{Ω}}\left(s\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{1}{k!}{\text{ad}}_{\text{Ω}\left(s\right)}^k\left(H(0)\right)\mathrm{.}$ (16) The Magnus formulation offers a significant advantage, as it enables the evaluation of arbitrary observables by utilizing the final Magnus operator $\text{Ω}(\infty )$, $O(\infty )\equiv e^{\text{Ω}\left(\infty \right)}O(0)e^{-\text{Ω}\left(\infty \right)}\mathrm{.}$ (17) The computational effort for solving the IMSRG(2) flow equations is primarily dictated by the two-body flow equation, which exhibits polynomial complexity of $\mathcal{O}\left(N^6\right)$ based on the single-particle size N.

The deformed IMSRG

The use of deformations as degrees of freedom in nuclear many-body problems can make the calculations more efficient [78, 79]. The standard IMSRG conserves spherical symmetry with a single reference, which works for closed-shell nuclei. To calculate open-shell nuclei, symmetry-breaking schemes have been developed, including both single- and multi-reference approaches. The single-reference Hartree-Fock-Bogoliubov (HFB) IMSRG, which selects a single HFB state as the reference state, has been proposed [48]. The HFB quasiparticle state breaks the particle number conservation, necessitating that particle number projection be performed. To choose a reference state closer to the true solution, the multi-reference IMSRG with particle-number-projected spherical HFB [49, 50] has been suggested. Calculations based on the Bogoliubov quasiparticle states significantly complicate the formalism and increase computational costs. Using the m scheme, a single HF reference state can be constructed for any even-even nuclei, with the particle number conserved but rotational symmetry broken. This deformed reference state may better reflect the intrinsic structure of a deformed nucleus and capture more correlations through symmetry restoration, which would otherwise be many-particle-many-hole excitations in the spherical scheme. The expected symmetry preconsiderations, e.g., as in the symmetry-adapted approach [80, 81], provide an efficient way to capture the expected features of nuclear states of interest while simultaneously reducing the computational cost.

As indicated in Sec. 2.1, the standard IMSRG is limited to extracting the ground-state energy of a closed-shell nucleus. An extension of the IMSRG to the deformed scheme would be useful for the description of open-shell nuclei. Therefore, we developed the D-IMSRG method [53] within the deformed HF basis, i.e., the m-scheme HF basis. First, the axially deformed HF equation of the even-even nucleus is solved within the spherical HO basis. Under the j-scheme, the initial Hamiltonian is typically expressed in the spherical HO basis. Using the Wigner-Eckart theorem, the matrix elements of operators, including the Hamiltonian in the j-scheme, can be converted to the matrix elements in the m-scheme, $\begin{array}{c}\langle {\alpha }^{\text{'}}{j}^{\text{'}}{m}^{\text{'}}\left|O_{kq}\right|\alpha jm\rangle =\langle jmkq|j^{\text{'}}{m}^{\text{'}}\rangle \\ \times \frac{1}{\sqrt{2j^{\text{'}}+1}}\langle {\alpha }^{\text{'}}{j}^{\text{'}}\Vert O_k\Vert \alpha j\rangle \mathrm{,}\end{array}$ (18) where $|\alpha jm\rangle$ are quantum states with good angular momentum j and projection m. The quantity α contains all other quantum numbers needed to completely specify the quantum state. Okq is a spherical tensor of rank k; $\langle jmkq|j^{\text{'}}{m}^{\text{'}}\rangle$ is the Clebsch-Gordan coefficient; and $\langle {\alpha }^{\text{'}}{j}^{\text{'}}\Vert O_k\Vert \alpha j\rangle$ denotes the reduced matrix elements. The one-, two-, and three-body matrix elements can thus be converted from j-scheme to m-scheme through Eq. (18). The m-scheme HF single-particle levels obtained are twofold degenerate with respect to the angular momentum projection quantum number m of the orbital (i.e., the energies are the same with respect to ±m). Filling the deformed HF single-particle levels up to the Fermi surfaces of neutrons and protons in ±m pairing from low to high $\left|m\right|$, keeps the axial, parity, and time-reversal symmetries of the even-even ground state, thereby creating an oblate or prolate deformed HF reference state [82]. Subsequently, the intrinsic A-body Hamiltonian in Eq. (5) is normal ordered with respect to the deformed A-dependent reference state $|\text{Φ}\rangle$ (i.e., the m-scheme HF ground state of the target nucleus). The off-diagonal parts of the Hamiltonian are consistent with the standard IMSRG, and the flow equations are evolved using the Magnus expansion Eqs. (16)-(17).

Subsequently, the ground-state energy and other observables can be calculated using the D-IMSRG ground-state wave function $|\text{Ψ}\rangle =e^{-\text{Ω}}|\text{Φ}\rangle$ (here $|\text{Φ}\rangle$ is the deformed HF reference state of the nucleus) expressed as $E=\langle \text{Ψ}\left|H\right|\text{Ψ}\rangle =\langle \text{Φ}\left|e^{\text{Ω}}H{e}^{-\text{Ω}}\right|\text{Φ}\rangle =\langle \text{Φ}\left|\tilde{H}\right|\text{Φ}\rangle \mathrm{,}$ (19) $O=\langle \text{Ψ}\left|O\right|\text{Ψ}\rangle =\langle \text{Φ}\left|e^{\text{Ω}}O{e}^{-\text{Ω}}\right|\text{Φ}\rangle =\langle \text{Φ}\left|\tilde{O}\right|\text{Φ}\rangle \mathrm{.}$ (20) In the D-IMSRG, the reference state is just the ground state of the deformed even-even nucleus. However, performing exact symmetry restoration of the D-IMSRG wave function is mathematically cumbersome and computationally expensive due to the exponential increase of configurations in projecting the wave function $|\text{Ψ}\rangle =e^{-\text{Ω}}|\text{Φ}\rangle$.

Therefore, an HF projection correction is introduced as a first approximation, to account for the angular momentum projection effect. The projection correction to the ground-state energy is estimated by $\Delta E_{\text{proj}}=\frac{\langle \text{Φ}\left|HP\right|\text{Φ}\rangle }{\langle \text{Φ}\left|P\right|\text{Φ}\rangle }-\frac{\langle \text{Φ}\left|H\right|\text{Φ}\rangle }{\langle \text{Φ}|\text{Φ}\rangle }\mathrm{,}$ (21) where $P_{M{M}^{\prime }}^J=\frac{2J+1}{8{\pi }^2}{\displaystyle \int }\text{d}\omega D_{M{M}^{\prime }}^{J*}\left(\omega \right)R(\omega )$ is the angular momentum projection operator. This provides a D-IMSRG ground-state energy given by $E+\Delta E_{\text{proj}}$ with the projection correction estimated by the HF wave function (here E is obtained by Eq. (19), that is, the ground-state energy without the projection).

A deformed coupled-cluster calculation has estimated that the angular-momentum projection of the HF state reduced the HF energy by approximately 3-5 MeV in the sd shell [82], corresponding to the static correlation, which is not size extensive. Since modern ab initio calculations already include some of the correlations associated with the projection, the energy correction obtained by projecting the ab initio wave function would be slightly smaller than the HF projection correction [82, 83].

In the spherical j-scheme, single-particle levels within the same j shell are degenerate. However, this degeneracy is broken with the onset of deformation although a twofold degeneracy with respect to ±m persists in axially symmetric shapes, significantly increasing the model-space dimension. The dimension of D-IMSRG calculation depends on the nucleon number A and basis-space size N_shell (the number of spherical HO major shells considered in solving the deformed HF). For ⁴⁰Mg, the number of D-IMSRG Hamiltonian matrix elements already exceeds 10⁹ at N_shell=10. Nonetheless, such a large model space may still not be sufficient to make the calculation converged. To estimate the converged ground-state energy, a simple exponential fitting method was applied with respect to N_shell to extrapolate the D-IMSRG result to an infinite basis space, similar to the ones used in NCSM [16, 84-87] and multi-reference IMSRG [19] calculations, $E\left(N_{\text{shell}}\right)=b_0+b_1\text{exp}\left(-b_2{N}_{\text{shell}}\right)\mathrm{,}$ (22) where b₀, b₁, and b₂ are the fitting parameters. The value of $b_0\equiv E\left(N_{\text{shell}}\to \infty \right)$ provides the estimate for the fully converged energy.

The valence-space IMSRG

The spherical j-scheme IMSRG can only treat closed-shell nuclei. The m-scheme D-IMSRG was designed to calculate open-shell nuclei. Unfortunately, the D-IMSRG does not conserve the angular momentum, and the exact angular-momentum projection D-IMSRG has not been available.

The nuclear shell model (SM) has served as one of the most powerful theoretical and computational tools for nuclear structure calculations [38, 88-91]. In the SM, valence nucleons move in the mean field generated by the inert core and interact via residual effective interactions. While the SM has been used predominantly in a phenomenological context [92, 91], there have been efforts dating back to decades ago to derive shell-model parameters based on a realistic interaction between nucleons [93, 94], the ab initio effective shell-model interactions. For open-shell systems, in addition to solving the full A-body problem, such as the D-IMSRG mentioned in Sec. 2.2, it is beneficial to follow the shell-model paradigm to construct and diagonalize the effective Hamiltonian in which the active degrees of freedom are Av valence nucleons confined to a few orbitals near the Fermi level. As for the IMSRG, a valence-space effective interaction can be derived using the spherical symmetry-conserving single-reference IMSRG at a shell closure to perform ab initio shell-model calculations for open-shell nuclei. This method, which combines the IMSRG and SM, is referred to as the VS-IMSRG [39].

The utility of the IMSRG lies in its flexibility to customize the definition of H^od to address specific problems. For the ground state of closed-shell nuclei, all terms that couple the reference state $|\text{Φ}\rangle$ to the rest of the Hilbert space can be eliminated, as in the standard IMSRG. For open-shell nuclei, $|{\text{Φ}}_v\rangle$ is decoupled from states containing non-valence states. This can be achieved by defining the Hod using the following matrix elements, $H^{\text{od}}=｛f_{ph}\mathrm{,}{f}_{pp\text{'}}\mathrm{,}{f}_{hh\text{'}}\mathrm{,}{\text{Γ}}_{pp\text{'}hh\text{'}}\mathrm{,}{\text{Γ}}_{pp\text{'}vh}\mathrm{,}{\text{Γ}}_{pqvv\text{'}}｝+\text{H}\text{.c}\text{.,}$ (23) where p=v,q. These off-diagonal parts of the generators evolve the Hamiltonian to diagonal $H^{\text{d}}$ form, where states outside the valence space are decoupled using the flow equation, as illustrated in Fig. 1, nonperturbatively satisfying the decoupling equation: $P{H}^{\text{d}}\left(\infty \right)Q=Q{H}^{\text{d}}\left(\infty \right)P=\mathrm{0,}$ (24) with $P={\displaystyle \sum {}_v}|{\text{Φ}}_v\rangle \langle {\text{Φ}}_v|$ and $Q=1-P$. After the evolution is complete, the effective shell-model Hamiltonian is obtained. The last step is to use the SM code, such as KSHELL [95], to diagonalize the effective Hamiltonian and express it as a reduced eigenvalue problem in the valence-particle space.

Fig. 1Full-screen View

(Color online) Schematic illustration of the VS-IMSRG decoupling from the initial Hamiltonian H(0) to obtain the final Hamiltonian $H(\infty )$ for the two valence nucleons

The current VS-IMSRG is at two-body approximation without explicitly considering 3NF or three-body correlations. To reduce the residual 3NF effect, a fractional filling of open-shell orbitals in an open-shell nucleus, named ensemble normal ordering (ENO), has been suggested [40]. Using the ENO approximation of the VS-IMSRG, nucleus-dependent valence-space effective Hamiltonian and effective operators of other observables can be obtained.

The Gamow IMSRG with coupling to continuum

Weakly bound and unbound nuclei belong to the category of open quantum systems, where coupling to the particle continuum profoundly affects the system behavior [96, 97]. Many novel phenomena, including halos [98, 99], genuine intrinsic resonances [100, 101], and new collective modes [102-104], have been observed or predicted in exotic nuclei. However, the majority of IMSRG calculations are performed within the HO or real-energy HF basis. Here the real-energy HF means that the HF approach is performed under the HO basis, which is bound and localized and hence isolated from the environment of unbound scattering states because of the Gaussian falloff of the HO functions. Similarly, the real-energy HF basis cannot include the continuum effect in IMSRG calculations.

The complex-energy Berggren basis offers an elegant framework for treating bound, resonant, and scattering continuum states on an equal footing [55, 56]. This basis is a generalization of the standard completeness relation from the real-energy axis to the complex-energy plane. Completeness encompasses a finite set of bound and resonance states together with a complex-energy scattering continuum: $\begin{array}{l}{\displaystyle \sum _n{u}_n}\left(E_n\mathrm{,}r\right)u_n\left(E_n\mathrm{,}{r}^{\text{'}}\right)\\ +{\displaystyle {\int }_L}\text{ d}Eu(E\mathrm{,}r)u\left(E\mathrm{,}{r}^{\text{'}}\right)=\delta \left(r-r^{\text{'}}\right)\text{,}\end{array}$ (25) where $u_n\left(E_n\mathrm{,}r\right)\sim O_l\left(k_{n}r\right)\sim {\text{e}}^{\text{i}{k}_{n}r}\mathrm{,}$ (26) Here, $k_n=\text{i}{\kappa }_n\left({\kappa }_n>0\right)$ for bound states and $k_n={\gamma }_n-\text{i}{\kappa }_n\left({\kappa }_n\mathrm{,}{\gamma }_n>0\right)$ for decaying resonances located in the fourth quadrant of the complex-momentum (complex-k) plane.

In practical applications, it is more convenient to express Eq. (25) in momentum space, ${\displaystyle \sum _{n\in (b\mathrm{,}d)}|u_n\rangle }\langle u_n|+{\displaystyle {\int }_{L^{+}}}|u(k)\rangle \langle u(k)|\text{d}k=1.$ (27)

As depicted in Fig. 2, bound and resonant states appear as poles of the scattering matrix within the complex-k plane, and the scattering continuum is represented by the blue contour. Within the Berggren basis, the GSM [23, 100, 101, 105-115] and complex coupled cluster [116, 18] methods have been well developed and widely applied to the calculations of weakly bound and unbound nuclei.

Fig. 2Full-screen View

(Color online) Location of one-body states in the complex-k plane. The Berggren completeness relation in Eq. (25) involves the bound states (brown-filled circles) lying on the imaginary k-axis, scattering states lying on the contour (solid blue line), and decaying resonant states (blue-filled circles) in the fourth quarter of the complex-k plane lying between the real axis and scattering contour. The capturing states (purple-filled circles) and antibound states (cyan-filled circles) are not included in the present completeness relation

For calculations within the Gamow-Berggren framework, selecting an appropriate one-body potential is essential for generating resonance and the continuum Berggren basis, frequently using the phenomenological Woods-Saxon potential [100, 101, 105, 117, 111]. In our approach, we used the Gamow Hartree-Fock (GHF) method with chiral potentials to produce an ab initio Berggren single-particle basis, which is convenient for computing the Berggren basis using an analytical continuation of Schrödinger’s equation in complex-k space. The complex-k single-particle GHF equation is formulated as follows: $\begin{array}{c}\frac{{\hslash }^2{k}^2}{2\mu }{\psi }_i\left(k\right)+{\displaystyle {\int }_{L^{+}}}\text{d}{k}^{\text{'}}{k}^{\text{'}2}U\left(lj{k}^{\text{'}}k\right){\psi }_i\left(k^{\text{'}}\right)\\ =e_i{\psi }_i\left(k\right)\mathrm{,}\end{array}$ (28) where $\mu =m/\left(1-\frac{1}{A}\right)$, and $k\left(k^{\prime }\right)$ is defined on the scattering contour. $U(lj{k}^{\text{'}}k)$ is the GHF single-particle potential, $U\left(lj{k}^{\text{'}}k\right)=\langle k\left|U\right|k^{\text{'}}\rangle ={\displaystyle \sum _{\alpha \beta }\langle k^{\text{'}}|\alpha \rangle }\langle \alpha \left|U\right|\beta \rangle \langle \beta |k\rangle \mathrm{,}$ (29) where l, j are the orbital and total angular momenta of a single-particle orbital, respectively. Greek letters denote HO states, indicating that 〈β/k〉 is the HO basis wave functions $|\beta \rangle$ expressed in the complex-k plane $\begin{array}{c}\langle \beta |k\rangle ={\left(-i\right)}^{2n+l}{e}^{-l/2b^2{k}^2}{\left(bk\right)}^l\\ \times \sqrt{\frac{2n!b^3}{\text{Γ}\left(n+l+3/2\right)}}{L}_n^{l+1/2}\left(b^2{k}^2\right)\mathrm{.}\end{array}$ (30) Here $L_n^{l+1/2}$ denotes the generalized Laguerre polynomial; $b=\sqrt{\hslash /m\omega }$; and ω is the frequency of the oscillator basis. Since $L_n^{l+1/2}$ is analytic, it can be extended to the complex momentum space, expressing its real part as $\text{Re}\left[L_n\left(x+iy\right)\right]={\displaystyle \sum _{j=0}^{\left|\frac{n}{2}\right|}\frac{{\left(-1\right)}^j{y}^{2j}}{\left(2j\right)!}}{L}_{n-2}^{2j}\left(x\right)\mathrm{,}$ (31) and its imaginary part as $\text{Im}\left[L_n\left(x+iy\right)\right]={\displaystyle \sum _{j=0}^{\left|\frac{n-1}{2}\right|}\frac{{\left(-1\right)}^{j-1}{y}^{2j+1}}{(2j+1)!}}{L}_{n-2j-1}^{2j+1}\left(x\right)\mathrm{,}$ (32) where x and y represent the real and imaginary parts of b²k², respectively; $\langle \alpha \left|U\right|\beta \rangle$ is the HF single-particle potential which can be obtained by solving the real-energy HF equation in the HO basis $\begin{array}{c}\langle \alpha \left|U\right|\beta \rangle ={\displaystyle \sum _{i=1}^A{\displaystyle \sum _{\gamma \delta }\langle \alpha \gamma \left|v^{N\text{}N}\right|\beta \delta \rangle }}{D}_{\gamma i}^{*}{D}_{\delta i}\\ +\frac{1}{2}{\displaystyle \sum _{i\mathrm{,}j=1}^A{\displaystyle \sum _{\gamma \delta ϵ\zeta }\langle \alpha \gamma ϵ\left|v^{3N}\right|\beta \delta \zeta \rangle }}{D}_{\gamma i}^{*}{D}_{\delta i}{D}_{ϵj}^{*}{D}_{\zeta j}\mathrm{,}\end{array}$ (33) where D is the coefficient of the HF single-particle state. In numerical calculations, the GHF equation is solved using the Gauss-Legendre quadrature scheme [118, 107, 112] with discrete points on the contour L⁺, $\frac{{\hslash }^2{k}_{\alpha }^2}{2\mu }{\psi }_i\left(k_{\alpha }\right)+{\displaystyle \sum _{\beta }{\omega }_{\beta }}{k}_{\beta }^2\langle k_{\alpha }\left|U\right|k_{\beta }\rangle {\psi }_i\left(k_{\beta }\right)=e_i{\psi }_i\left(k_{\alpha }\right)\mathrm{.}$ (34) Here kα (kβ) are the discrete momentum points, and ωα (ωβ) are the corresponding Gauss-Legendre quadrature weights. We define ${\psi }_i^{\text{'}}\left(k_{\alpha }\right)={\psi }_i\left(k_{\alpha }\right)k_{\alpha }\sqrt{{\omega }_{\alpha }}\mathrm{.}$ (35) Then, Eq. (34) can be written as ${\displaystyle \sum _{\beta }{h}_{\alpha \beta }}{\psi }_i^{\text{'}}\left(k_{\beta }\right)=e_i{\psi }_i^{\text{'}}\left(k_{\alpha }\right)\mathrm{,}$ (36) with $h_{\alpha \beta }=\frac{{\hslash }^2}{2\mu }{k}_{\alpha }^2{\delta }_{\alpha \beta }+\sqrt{{\omega }_{\alpha }{\omega }_{\beta }}\langle k_{\alpha }\left|U\right|k_{\beta }\rangle \mathrm{.}$ (37) The bound, resonant, and continuum GHF basis can be obtained by diagonalizing the complex-energy in Eq. (36).

Within the GHF basis, the G-IMSRG calculations can be conducted. Notably, the Hamiltonian in real-energy space is Hermitian, $H=H^{†}$. Therefore, in practice, the similarity transformation U(s) is a unitary transformation, satisfying $U\left(s\right)U^{†}\left(s\right)=U\left(s\right)U^{-1}\left(s\right)=1$, and $\eta \left(s\right)=\frac{\text{d}U\left(s\right)}{\text{d}s}{U}^{†}\left(s\right)=-{\eta }^{†}\left(s\right)$ is the anti-Hermitian generator. However, in the G-IMSRG framework, the Hamiltonian is complex symmetric, $H=H^{\text{T}}$ (here T indicates the transpose). Therefore, employing a continuous orthogonal transformation, $U\left(s\right)U^{\text{T}}\left(s\right)=U\left(s\right)U^{-1}\left(s\right)=1$, and the Hamiltonian H(s) can be transformed into a band or block diagonal form, $H\left(s\right)=U\left(s\right)H(0)U^{\text{T}}\left(s\right)\mathrm{.}$ (38) Correspondingly, the generator $\eta (s)$ becomes $\eta \left(s\right)=\frac{\text{d}U(s)}{\text{d}s}{U}^{\text{T}}\left(s\right)=-{\eta }^{\text{T}}\left(s\right)\mathrm{.}$ (39) The G-IMSRG method can directly compute the ground state of a closed-shell nucleus by decoupling the Hamiltonian from the excitations above the closed-shell Fermi surface. To handle open-shell nuclei or excited states, we employed G-IMSRG using the equation-of-motion (EOM) approach [119]. This approach offers a useful alternative to the shell model strategy for calculating excited states, especially when extended valence spaces lead to prohibitively large shell-model basis dimensions. Within the EOM framework, the Schrödinger equation is rewritten using ladder operators, which create excited eigenstates from the exact ground state, $H|{\text{Ψ}}_n\rangle =E_n|\text{Ψ}\rangle \to H{X}_n^{†}|{\text{Ψ}}_0\rangle =E_n{X}_n^{†}|{\text{Ψ}}_0\rangle \mathrm{,}$ (40) where $X_n^{†}$ is given by the dyadic product $|{\text{Ψ}}_n\rangle \langle {\text{Ψ}}_0|$. Further rewriting Eq. (40) as EOM $\left[H\mathrm{,}{X}_n^{†}\right]|{\text{Ψ}}_0\rangle =\left(E_n-E_0\right)X_n^{†}|{\text{Ψ}}_0\rangle \equiv {\omega }_n{X}_n^{†}|{\text{Ψ}}_0\rangle \mathrm{.}$ (41) The amplitudes of $X_n^{†}$ are determined by solving a generalized eigenvalue problem [120].

Coupling EOM methods with G-IMSRG is natural, as the reference state $|{\text{Φ}}_0\rangle$ corresponds to the ground state of $\overline{H}\equiv U\left(\infty \right)H{U}^T\left(\infty \right)$. Multiplying Eq. (41) by $U\left(\infty \right)$ and recalling that $U(\infty )|{\text{Ψ}}_0\rangle =|{\text{Φ}}_0\rangle \mathrm{,}$ (42) we obtain the similarity transformed EOM $\left[\overline{H}\mathrm{,}{\overline{X}}_n^{†}\right]|{\text{Φ}}_0\rangle ={\omega }_n{\overline{X}}_n^{†}|{\text{Φ}}_0\rangle \mathrm{,}$ (43) where ${\overline{X}}_n^{†}\equiv U\left(\infty \right)X_n^{†}{U}^T\left(\infty \right)$. The solutions ${\overline{X}}_n^{†}$ can then be used to obtain the eigenstates of the unevolved Hamiltonian via $|{\text{Ψ}}_n\rangle =U^T\left(\infty \right){\overline{X}}_n^{†}|{\text{Φ}}_0\rangle \mathrm{.}$ (44) Currently, in our applications, we include up to 2p2h excitations in the ladder operator [120] ${\overline{X}}_n^{†}={\displaystyle \sum _{ph}{\overline{X}}_{ph}^{\left(n\right)}}｛a_p^{†}{a}_h｝+\frac{1}{4}{\displaystyle \sum _{pp\text{'}hh\text{'}}{\overline{X}}_{pp\text{'}hh\text{'}}^{\left(n\right)}}｛a_p^{†}{a}_{p\text{'}}^{†}{a}_{h\text{'}}{a}_h｝\mathrm{.}$ (45) In principle, the EOM ladder operator can include any excitation rank up to ApAh, which would constitute an exact diagonalization of $\overline{H}$ and can be computationally expensive. In practical applications, the EOM-G-IMSRG method is commonly employed in an approximative systematically improbable form, referred to as EOM-G-IMSRG(m,n), where m and n denote the truncation level in EOM and G-IMSRG, respectively. The calculations in the present work were performed using the EOM-G-IMSRG(2,2) approximation.

The Gamow coupled-channel method

To provide a thorough description of open quantum systems, the GCC method [67, 68, 70] has been advanced as an alternative approach, focusing on the few-body decay processes influenced by continuum and structural factors [121, 122]. The methodology involves constructing a robust three-body framework, utilizing the Berggren basis [55]. As elucidated in Sec. 2.4, the Berggren basis is specifically designed to incorporate continuum effects, thereby facilitating the analysis of weakly bound and unbound nuclear systems.

2.5.1

Spherical system with an inert core

In the context of the three-body GCC model, the nucleus comprises a core and two valence nucleons or clusters. The GCC Hamiltonian is formulated as follows: $H={\displaystyle \sum _{i=1}^3\frac{{\widehat{\overline{p}}}_i^2}{2m_i}}+{\displaystyle \sum _{i>j=1}^3{V}_{ij}}\left({\overrightarrow{r}}_{ij}\right)-{\widehat{T}}_{\text{c}\mathrm{.}\text{m}\text{.}}\mathrm{,}$ (46) where Vij denotes the interaction between clusters i and j, and ${\widehat{T}}_{\text{c}\text{.m}\text{.}}$ represents the kinetic energy of the center-of-mass. Each i-th cluster ($i=c\mathrm{,}{n}_1\mathrm{,}{n}_2$) is characterized by its position vector ${\overrightarrow{r}}_i$ and linear momentum ${\overrightarrow{k}}_i$. To accurately describe three-body asymptotics and eliminate the spurious center-of-mass motion, it is advantageous to utilize the relative (Jacobi) coordinates: $\begin{array}{l}\overrightarrow{x}=\sqrt{{\mu }_x}\left({\overrightarrow{r}}_{i_1}-{\overrightarrow{r}}_{i_2}\right)\mathrm{,}\\ \overrightarrow{y}=\sqrt{{\mu }_y}\left(\frac{A_{i_1}{\overrightarrow{r}}_{i_1}+A_{i_2}{\overrightarrow{r}}_{i_2}}{A_{i_1}+A_{i_2}}-{\overrightarrow{r}}_{i_3}\right)\mathrm{,}\end{array}$ (47) where $i_1=n_1$, $i_2=n_2$ and $i_3=c$ for T-coordinates, whereas $i_1=n_2$, $i_2=c$ and $i_3=n_1$ for Y-coordinates, as depicted in Fig. 3. Here, Ai represents the mass number of the i-th cluster; ${\mu }_x=\frac{A{i}_1{A}_{i_2}}{A_{i_1}+A_{i_2}}$ and ${\mu }_y=\frac{\left(A_{i_1}+A_{i_2}\right)A_{i_3}}{A_{i_1}+A_{i_2}+A_{i_3}}$ are the reduced masses associated with $\overrightarrow{x}$ and $\overrightarrow{y}$, respectively. For analytical convenience, the hyperradius $\rho =\sqrt{x^2+y^2}$, which remains invariant across different Jacobi coordinate transformations, is introduced.

Fig. 3Full-screen View

(Color online) Illustration of the coordinate and momentum configurations in a core + nucleon + nucleon system: (a) Jacobi T (solid lines) and Y (dashed lines) coordinates, where the former is used to describe the interactions between the nucleons (n₁ and n₂) and the latter for interactions involving the core (c). (b) Corresponding momentum scheme within the c.m. frame. Here, A denotes the mass number; μij represents the reduced mass between clusters i and j; and k₁, k₂, and kc indicate the momenta of nucleons n₁, n₂, and core c, respectively. The figure is taken from [70]

Experimental measurements in the momentum space necessitate the definition of relative momenta as follows: $\begin{array}{l}{\overrightarrow{k}}_x={\mu }_x\left(\frac{{\overrightarrow{k}}_{i_1}}{A_{i_1}}-\frac{{\overrightarrow{k}}_{i_2}}{A_{i_2}}\right)\mathrm{,}\\ {\overrightarrow{k}}_y={\mu }_y\left(\frac{{\overrightarrow{k}}_{i_1}+{\overrightarrow{k}}_{i_2}}{A_{i_1}+A_{i_2}}-\frac{{\overrightarrow{k}}_{i_3}}{A_{i_3}}\right)\mathrm{.}\end{array}$ (48) In the absence of c.m. motion, it is evident that ${\displaystyle {\sum }_i{\overrightarrow{k}}_i}=0$, and notably, ${\overrightarrow{k}}_y$ is oriented in the opposite direction to ${\overrightarrow{\theta }}_{i_3}$. The angles θk and ${{\theta }^{\prime }}_k$ represent the opening angles of the vectors (${\overrightarrow{k}}_x$, ${\overrightarrow{k}}_y$) in T- and Y-Jacobi coordinates, respectively (refer to Fig. 3). For example, in the scenario of two-nucleon decay, the kinetic energy associated with the relative motion of the emitted nucleons is expressed as $E_{\text{pp/nn}}=\frac{{\hslash }^2{k}_x^2}{2{\mu }_x}$, with E_core-p/n pertaining to the kinetic energy of the core-nucleon pair. The distribution types T (θk, E_pp/nn) and Y (${{\theta }^{\prime }}_k$, E_core-p/n) elucidate the nucleon-nucleon correlations and provide insights into the structural dynamics of the progenitor nucleus. The total momentum k is defined as $\sqrt{\frac{k_x^2}{{\mu }_x}+\frac{k_y^2}{{\mu }_y}}$, which asymptotically approaches $\frac{\sqrt{2m{Q}_{\text{2p/2n}}}}{\hslash }$ as time progresses, where Q_2p/2n is the two-nucleon decay energy derived from the binding energy difference between parent and daughter nuclei.

The presence of Pauli forbidden states in three-body models represents a challenge arising from the lack of antisymmetrization between core and valence particles. To address this issue, the orthogonal projection method [123-125], which entails the inclusion of a Pauli operator in the GCC Hamiltonian, was adopted and formulated as: $\widehat{Q}=\text{Λ}{\displaystyle \sum _c|{\phi }^{j_c{m}_c}\rangle }\langle {\phi }^{j_c{m}_c}|\mathrm{,}$ (49) where Λ is a large constant, and $|{\phi }^{j_c{m}_c}\rangle$ represents a two-body state comprising forbidden single-particle (s.p.) states of core nucleons. By setting Λ to high values, Pauli forbidden states are elevated to higher energies, which effectively suppresses their influence within the physical spectrum of the system.

This standard projection technique [67] may introduce minor numerical errors in the asymptotic region because of coordinate transformations. The supersymmetric transformation method [125-127] offers an alternative solution for the exclusion of Pauli-forbidden states. This method employs an auxiliary repulsive “Pauli core” within the original core-proton interaction, thereby effectively eliminating the influence of Pauli-forbidden states from the system.

The total wave function is expressed using hyperspherical harmonics as: ${\text{Ψ}}^{JM\pi }\left(\rho \mathrm{,}{\text{Ω}}_5\right)={\rho }^{-5/2}{\displaystyle \sum _{\gamma K}{\psi }_{\gamma K}^{J\pi }}\left(\rho \right){\mathcal{Y}}_{\gamma K}^{JM}\left({\text{Ω}}_5\right)\mathrm{,}$ (50) where K denotes the hyperspherical quantum number. The set $\gamma =｛s_1\mathrm{,}{s}_2\mathrm{,}{s}_3\mathrm{,}{S}_{12}\mathrm{,}S\mathrm{,}{\mathcal{l}}_x\mathrm{,}{\mathcal{l}}_y\mathrm{,}L｝$ encapsulates the quantum numbers other than K. Here, s represents spin, and $\mathcal{l}$ denotes orbital angular momentum. The function ${\psi }_{\gamma K}^{J\pi }\left(\rho \right)$ specifies the hyperradial wave function, and ${\mathcal{Y}}_{\gamma K}^{JM}\left({\text{Ω}}_5\right)$ represents the corresponding hyperspherical harmonic [125].

The resulting Schrödinger equation for the hyperradial wave functions can be written as a set of coupled-channel equations: $\begin{array}{l}\left[-\frac{{\hslash }^2}{2m}\left(\frac{{\text{d}}^2}{\text{d}{\rho }^2}-\frac{(K+3/2)(K+5/2)}{{\rho }^2}\right)-\tilde{E}\right]{\psi }_{\gamma K}^{J\pi }\left(\rho \right)\\ +{\displaystyle \sum _{K^{\prime }{\gamma }^{\prime }}{V}_{K^{\prime }{\gamma }^{\prime }\mathrm{,}K\gamma }^{J\pi }}\left(\rho \right){\psi }_{{\gamma }^{\prime }{K}^{\prime }}^{J\pi }\left(\rho \right)\\ +{\displaystyle \sum _{K^{\prime }{\gamma }^{\prime }}{\displaystyle {\int }_0^{\infty }}}{W}_{K^{\prime }{\gamma }^{\prime }\mathrm{,}K\gamma }\left(\rho \mathrm{,}{\rho }^{\prime }\right){\psi }_{{\gamma }^{\prime }{K}^{\prime }}^{L\pi }\left({\rho }^{\prime }\right)\text{d}{\rho }^{\prime }=\mathrm{0,}\end{array}$ (51) where $V_{K^{\prime }{\gamma }^{\prime }\mathrm{,}K\gamma }^{L\pi }\left(\rho \right)=\langle {\mathcal{Y}}_{{\gamma }^{\prime }{K}^{\prime }}^{JM}\left|{\displaystyle \sum _{i>j=1}^3{V}_{ij}}\left({\overrightarrow{r}}_{ij}\right)\right|{\mathcal{Y}}_{\gamma K}^{JM}\rangle$ (52) and $W_{K^{\prime }{\gamma }^{\prime }\mathrm{,}K\gamma }\left(\rho \mathrm{,}{\rho }^{\prime }\right)=\langle {\mathcal{Y}}_{{\gamma }^{\prime }{K}^{\prime }}^{JM}\left|\widehat{Q}\right|{\mathcal{Y}}_{\gamma K}^{JM}\rangle$ (53) is the non-local potential generated by the Pauli projection operator, as defined in Eq. 49.

To properly treat the positive-energy continuum space, the Berggren expansion technique is utilized for the hyperradial wave function: ${\psi }_{\gamma K}^{J\pi }\left(\rho \right)={\displaystyle \sum _{\text{n}}{C}_{\gamma \text{n}K}^{J\pi M}}{\mathcal{B}}_{\gamma \text{n}}^{J\pi }\left(\rho \right)\mathrm{,}$ (54) where ${\mathcal{B}}_{\gamma \text{n}}^{J\pi }\left(\rho \right)$ denotes an s.p. state within the Berggren ensemble [55] (detailed in Sec. 2.4). To compute radial matrix elements using the Berggren basis, exterior complex scaling [128] is employed, whereby integrals are evaluated along a complex radial trajectory: $\begin{array}{l}\langle {\mathcal{B}}_{\text{n}}|V(\rho )|{\mathcal{B}}_{\text{m}}\rangle ={\displaystyle {\int }_0^R}{\mathcal{B}}_{\text{n}}\left(\rho \right)V\left(\rho \right){\mathcal{B}}_{\text{m}}\left(\rho \right)\text{d}\rho \\ +{\displaystyle {\int }_0^{+\infty }}{\mathcal{B}}_{\text{n}}\left(R+\rho e^{i\theta }\right)V\left(R+\rho e^{i\theta }\right){\mathcal{B}}_{\text{m}}\left(R+\rho e^{i\theta }\right)\text{d}\rho \mathrm{.}\end{array}$ (55) For potentials that decay as $O\left(1/{\rho }^2\right)$ (such as the centrifugal potential) or more rapidly (such as the nuclear potential), R should be large enough to circumvent all singularities, with the scaling angle θ selected to ensure that the integrals converge (see [129] for further details). Since the Coulomb potential is not square-integrable, its matrix elements diverge when the complex momenta $k_n=k_m$. To address this, the “off-diagonal method” introduced in [130], where a slight offset $\pm \delta k$ is added to the linear momenta of the involved scattering wave-functions, was applied to facilitate the convergence of the resulting diagonal Coulomb matrix elements. The complex-momentum representation has also been adopted in other methods, e.g., in mean-field calculations [131, 132], to describe the continuum effect.

Time-dependent approach

To tackle the decay process, a time-dependent formalism was developed, to allow precise, numerically stable, and transparent investigations of a broad range of phenomena, such as configuration evolution [145], decay rates [146], and fission [147]. For two-nucleon decay, the measured inter-particle correlations can be interpreted by propagating the solutions over long times. Despite previous efforts in this direction [148-150], capturing the asymptotic correlation of emitted particles still requires a precise description of the resonance wave function at large distances. To this end, we utilized the complex-momentum state ${\text{Ψ}}_0^{J\pi }$ obtained using the GCC method. This state can be decomposed into real-momentum scattering states using the Fourier-Bessel series expansion in the real-energy Hilbert space [151]. The resulting wave packet is propagated by the time evolution operator through the Chebyshev expansion [152, 153]: $e^{-i\frac{\widehat{H}}{\hslash }t}={\displaystyle \sum _{n=0}^{\infty }{\left(-i\right)}^n}\left(2-{\delta }_{n0}\right)J_n\left(t\right)T_n\left(\widehat{H}/\hslash \right)\mathrm{,}$ (57) where Jn are the Bessel functions of the first kind and Tn are the Chebyshev polynomials.

The time evolution was limited to the real momentum space, to restore the Hermitian property of the Hamiltonian matrix and ensure conservation of total density. Our implementation of the time-dependent approach is based on the integral equation, which allows maintaining high numerical precision by utilizing the Chebyshev expansion’s good convergence rate [154, 153]. Furthermore, the evolving wave packet has an implicit cutoff at large distances, preventing the divergence of the Coulomb interaction in momentum space. In practice, we considered interactions within a sphere of radius of approximately 500 fm, and the wave function remained defined in momentum space beyond this cutoff, preventing unwanted reflections at the boundary.

The D-IMSRG calculations of deformed light nuclei

In [53], Yuan and his collaborators developed D-IMSRG, as formulated in Sec. 2.2, which better reflects the intrinsic structure of the deformed nucleus and captures more correlations through symmetry restoration. As a test ground, the D-IMSRG was first performed to calculate the ground-state energies of ⁸Be and ¹⁰Be, which are exotic nuclei with 2α cluster structure or elongated shapes, benchmarked against NCSM and VS-IMSRG. Subsequently, D-IMSRG was applied to describe the ground-state properties of even-even nuclei ranging from light beryllium to medium-mass magnesium isotopes. The optimized chiral NN interaction NNLO_opt [156, 157], which gives good descriptions of nuclear binding energies, excitation spectra and neutron matter equation of state without the inclusion of the 3N force, was used during the calculation in [53] with $\hslash \omega =24\text{ MeV}$.

The ground-state energies of ⁸Be and ¹⁰Be were first studied though D-IMSRG, as shown in Fig. 4, with and without the approximate angular momentum projection. It was found out that the trend of calculated energy by D-IMSRG is similar to those of NCSM [16, 84-87] and multi-reference IMSRG [19] calculations, exhibiting an exponential convergence with respect to the basis-space size N_shell. The energies extrapolated to infinite model space using an exponential fit based on Eq. (22) are depicted in Fig. 4. The extrapolated “Extrap” results fitted with different data points are provided in Fig. 4 along with their uncertainties, verifying that the calculation results of D-IMSRG converge exponentially with an increase in N_shell, thereby demonstrating the reliability of the calculations.

Fig. 4Full-screen View

(Color online) Ground-state energies for ⁸Be and ¹⁰Be calculated by D-IMSRG, with and without projection correction, are shown as a function of basis-space size N_shell. Symbols below “Extrap” represent the energies extrapolated to the infinite basis space using an exponential fit, based on different data points: N_shell=3-7, 3-10, and 6-10. The fitting uncertainties are indicated by error bars. Extrapolation uncertainties in NCSM and VS-IMSRG calculations are also represented by error bars. Experimental data were taken from AME2020 [155], and the figure was taken from [53]

In Fig. 4, the angular momentum projection corrections are -7.2 MeV and -5.9 MeV for ⁸Be and ¹⁰Be, respectively. These significant corrections are caused by the large deformations of these two nuclei. The results of NCSM and VS-IMSRG and the experimental data are also shown in Fig. 4.

The results of VS-IMSRG underestimate the ground-state energy in ^8,10Be, which may be attributed to the omission of higher-order collective excitations that are not handled well in VS-IMSRG at the IMSRG(2) level, as discussed in [158, 159]. However, this omission can be compensated by the angular momentum projection correction through D-IMSRG, as illustrated in Fig. 4.

The ground-state energies and two-neutron separation energies calculated by D-IMSRG for ^6-16Be are shown in Fig. 5 (top panel and bottom panel, respectively), along with VS-IMSRG calculations and experimental data. The angular momentum projection lowered the ground-state energies of ^8-16Be by about 5-6 MeV, making the calculated energies closer to the data. Both D-IMSRG and VS-IMSRG calculations indicated that the neutron dripline of beryllium isotopes was at ¹²Be, contrary to the experimental position of ¹⁴Be, which may be due to the absence of a continuum effect [105, 107, 57, 160].

Fig. 5Full-screen View

(Color online) Ground-state energies (upper panel) and two-neutron separation energies S_2n (lower panel) of ^6-16Be calculated by D-IMSRG with and without projection correction. The VS-IMSRG results and experimental data [155] are also presented for comparison. The figure is taken from [53]

In [53], the heavier nuclei of C, O, Ne, and Mg isotopes were also calculated by D-IMSRG, as shown in Fig. 6 along with VS-IMSRG calculations and experimental data. The D-IMSRG calculations with the projection correction agreed well with VS-IMSRG results and experimental data. The angular momentum projection corrections were zero for the closed-shell nuclei ¹⁴C and ^{14,16,22,24,28}O, indicating the spherical characteristics of these nuclei. However, for the expected closed-shell nuclei of ^12,22C, the angular momentum projection corrections were not zero but -5.5 MeV and -2.7 MeV, respectively, indicating their deformation. For Ne and Mg isotopes, the projection results provided energy gains of about 3-6 MeV near the neutron number N=20 island of inversion [161, 162]. There is strong configuration mixing between sd and pf shells in nuclei located in the region of the island of inversion. This cross-shell mixing is missing in the VS-IMSRG calculations although a multi-shell VS-IMSRG has been proposed [163]. Therefore, it can be concluded that deformation effectively brings the deformation orbitals into the wave function of the state in D-IMSRG calculations.

Fig. 6Full-screen View

(Color online) Ground-state energies of C, O, Ne, and Mg isotopes. D-IMSRG results are extrapolated to the infinite basis space using the N_shell=6-10 data points, whereas the VS-IMSRG results are extrapolated based on N_shell=8-13. In VS-IMSRG calculations, the model space includes both protons and neutrons in 0p_3/2,1/2 for ^6-14C; protons in 0p_3/2,1/2 and neutrons in 1s_1/20d_5/2,3/2 for ^14-22C; and both protons and neutrons in 1s_1/20d_5/2,3/2 for O, Ne, and Mg isotopes. The experimental data were taken from AME2020 [155]. The figure was taken from [53]

For nuclei, another important observable is the charge radius. The radii of the studied isotopes were also calculated in [53]. Compared with the calculation of the ground-state energy, the convergence of the radius calculation showed a different trend with increasing basis-space size, as discussed in [32, 168, 87] and as also observed in the D-IMSRG calculations [53], indicating that the exponential fit was not applicable to the radius. Therefore, the authors in [53] did not use extrapolation to fit the radius in the calculations. The angular momentum projection correction was also estimated by the HF wave function in the study. As shown in Fig. 7, the charge radii of Be, C, O, Ne, and Mg isotopes were investigated, along with the VS-IMSRG calculations and experimental data. The projection correction to the charge radius is small, making the D-IMSRG radii with and without the projection correction close to each other and in good agreement with the VS-IMSRG calculations except for ⁸Be, where the D-IMSRG radius is larger. This difference can be attributed to the large deformation of ⁸Be. Therefore, it can be concluded that the calculated charge radii by D-IMSRG and VS-IMSRG are reasonable compared with experiment data, as shown in Fig. 7 although the NNLO_opt interaction tends to underestimate the nuclear radii, as noted in [157].

Fig. 7Full-screen View

(Color online) Charge radii calculated by D-IMSRG in a basis space with N_shell=10 and using VS-IMSRG with N_shell=13 for Be, C, O, Ne, and Mg isotopes. Experimental data were taken from [157, 164-167]. The figure was taken from [53]

δV_pn bifurcation by the VS-IMSRG

The VS-IMSRG was first introduced by Tsukiyama et al. in 2012 [39]. This method combines the SM and IMSRG to nonperturbatively derive effective valence-space Hamiltonians and operators, as detailed in Sec. 3.2. Recently, it has been applied to describe the δV_pn values in the upper fp shell, incorporating a chiral three-nucleon force (3NF), as reported in [169].

Nuclear binding energy B(Z,N) represents the total interaction energy of interacting nucleons in the nucleus with Z protons and N neutrons. Differences in binding energy can isolate specific types of interactions and provide insights into modifications in nuclear structure [170, 171]. The double binding energy difference denoted as δV_pn, has been used as an important mass filter to extract the residual proton-neutron (pn) interaction [172-174], particularly for the N=Z nuclei. The residual proton-neutron interaction δV_pn can be extracted using $\begin{array}{c}\delta V_{\text{pn}}^{\text{ee}}\left(Z\mathrm{,}N\right)=\frac{1}{4}[B(Z\mathrm{,}N)-B(Z\mathrm{,}N-2)-B(Z-\mathrm{2,}N)\\ +B(Z-\mathrm{2,}N-2)],\end{array}$ (58) for the nucleus with N=Z= an even number, and $\begin{array}{c}\delta V_{\text{pn}}^{\text{oo}}\left(Z\mathrm{,}N\right)=[B\left(Z\mathrm{,}N\right)-B(Z\mathrm{,}N-1)-B(Z-\mathrm{1,}N)\\ +B(Z-\mathrm{1,}N-1)]\mathrm{,}\end{array}$ (59) for the nucleus with N=Z= an odd number.

Weakly bound proton-rich nuclei are attracting interest in novel structure [175]. In [169], the masses of ⁶²Ge, ⁶⁴As, ⁶⁶Se, and ⁷⁰Kr were measured for the first time. Additionally, the masses of six N=Z-1 nuclides ⁶¹Ga, ⁶³Ge, ⁶⁵As, ⁶⁷Se, ⁷¹Kr, and ⁷⁵Sr were redetermined with improved accuracy, using a novel method of isochronous mass spectrometry conducted at the Heavy Ion Research Facility in Lanzhou (HIRFL). These newly measured masses provide updated δV_pn values, which offer great test ground for state-of-the-art theoretical calculations. The updated δV_pn values show a clear increasing trend in $\delta V_{\text{pn}}^{\text{oo}}$ beyond Z=28, which is interpreted as an indication of the restoration of pseudo-SU(4) symmetry in the fp shell, as suggested in [176, 177]. In contrast, $\delta V_{\text{pn}}^{\text{ee}}$ shows a decreasing trend that was previously observed in the lower mass region [173, 178]. These δV_pn values were extracted using predicted masses from frequently used mass models; however, none of these models successfully reproduce the bifurcation in δV_pn values [169], except for the ab initio VS-IMSRG calculations. Within the ab initio VS-IMSRG calculation, a chiral 2NF plus 3NF, labeled by EM1.8/2.0 [24], is adopted, which can reproduce well the ground-state energies up to A≈100 region nuclei [179, 7, 28]. The effective Hamiltonian in the full fp-shell above the ⁴⁰Ca core was derived using the VS-IMSRG, and the final diagonalization of the valence-space Hamiltonian was realized using KSHELL [95].

As shown in the lower panel of Fig. 8, VS-IMSRG calculations with chiral 2NF plus 3NF reproduce the experimental δV_pn for the N=Z+2 nuclei exceptionally well. For the odd-odd N=Z nuclei, ${}^{62}\text{Ga}{\mathrm{,}}^{66}\text{As}{\mathrm{,}}^{70}\text{Br}$ and ⁷⁴Rb, the ground states have been identified as $\left(T\mathrm{,}{J}^{\pi }\right)=$ $\left({\mathrm{1,0}}^{+}\right)$ [155]. VS-IMSRG calculations, which inherently incorporate both T=0 and T=1 pn correlations, achieve a commendable match with the experimental $\delta V_{\text{pn}}^{\text{oo}}$ values for nuclei ranging from ⁵⁸Cu to ⁷⁰Br. Particularly noteworthy is that the calculation successfully reproduces the observed increasing trend in $\delta V_{\text{pn}}^{\text{oo}}$ with an increase in nucleon number A. Our calculations consistently attributed an isospin of T=1 to the ground states of these odd-odd nuclei, aligning with experimental results, with the exception of ⁵⁸Cu. Moreover, the decreased trend in the even-even $\delta V_{\text{pn}}^{\text{ee}}$ was also well reproduced by our VS-IMSRG calculations.

Fig. 8Full-screen View

(Color online) Experimental δV_pn for (a) N=Z and (b) N=Z+2 nuclei beyond A=56, compared with the ab initio VS-IMSRG calculations. Data uncertainties are indicated by the size of symbols. δV_pn values from ab initio calculations using 2NF + 3NF and only 2NF plotted as red and blue lines (solid lines for even-even and dashed lines for odd-odd), respectively. The figure is taken from [169]

In mass regions with extremely asymmetric N/Z ratios, 3NF usually provides a repulsive effect on the neutron-neutron (nn) and proton-proton (pp) interactions [180, 28], which is essential for the emergence of new magic numbers [180] and also for reproducing the neutron or proton dripline [28]. To understand the role played by 3NF in the δV_pn of upper fp-shell nuclei, we performed calculations using only a chiral 2NF at N3LO. The results using only 2NF show significant deviations from those calculated with 3NF included, as demonstrated in Fig. 8. Specifically, the agreement with experimental δV_pn values was markedly poor in the calculations without the 3NF included. Additionally, the predicted isospins of ground states for odd-odd nuclei ranging from ⁶²Ga to ⁷⁴Rb were all erroneously identified as T=0 in the calculations with the 3NF included, contradicting experimental data. Furthermore, without the inclusion of 3NF, the calculated $\delta V_{\text{pn}}^{\text{oo}}$ values of N=Z nuclei were lower than $\delta V_{\text{pn}}^{\text{ee}}$ calculated with 3NF included. The 3NF enhances the pn correlations in N=Z nuclei favoring a T=1 isospin coupling, which changes the δV_pn behavior.

The G-IMSRG with the coupling to the continuum

A novel G-IMSRG [57] was developed by Hu et al., using the complex-energy Berggren representation, as introduced in Sec. 2.4. This advanced G-IMSRG is capable of describing the weakly bound and unbound nature of nuclei in the vicinity of nuclear driplines. We applied G-IMSRG to oxygen and carbon isotopes. Recent experiments [181-184] highlight that ²²C is a Borromean halo nucleus, with an experimentally deduced root-mean-squared matter radius of 3.44±0.08 fm [184]. The continuum coupling plays a vital role in generating the extended density distribution. Notably, experimental information about the excited states of ²²C, which can offer additional insights into its halo structure, is lacking. In this study, we performed an ab initio G-IMSRG calculation of the halo ²²C, using both chiral 2NF NNLO_opt and 2NF plus 3NF NNLO_sat interactions. The NNLO_opt interaction matrix elements were expanded within 12 major HO shells at a frequency of $\hslash \omega =20\text{ MeV}$, whereas the NNLO_sat interaction was truncated at 13 major HO shells with $\hslash \omega =22\text{ MeV}$ [185, 158]. The NNLO_opt potential provides a good description of nuclear structure, including binding energies, excitation spectrum, and dripline position without the need for 3NF [156]. The NNLO_sat interaction can provide accurate descriptions of charge radii in light- and mid-mass isotopes [186].

For the sd shell, the neutron 0d_3/2 is a narrow-resonance orbital. With no centrifugal barrier of the l=0 s partial wave, a weakly bound 1s_1/2 orbital can significantly affect the spatial distribution of the wave functions of weakly bound nuclei. Therefore, the 0d_3/2 and 1s_1/2 orbitals should be treated in the Berggren basis, which includes coupling to the continuum, whereas other orbitals can be treated in the real-energy HF basis to reduce the computational cost, as in [57].

Although the Hamiltonian (5) is intrinsic, the IMSRG wave functions are expressed in laboratory coordinates. Therefore, considering center-of-mass (CoM) motion corrections may be necessary. Our previous work has indicated that the CoM effect on an intrinsic Hamiltonian is small for low-lying states [105]. Thus, the approximation method suggested in [187, 37] can be adopted to estimate the CoM effect in IMSRG calculations. Fig. 9 presents real-energy EOM-IMSRG calculations without and with the CoM multiplication term $\beta H_{\text{CoM}}=\beta \left(\frac{{\text{P}}^2}{2mA}+\frac{1}{2}mA{\tilde{\omega }}^2{\text{R}}^2-\frac{3}{2}\hslash \tilde{\omega }\right)$. Note that the value of the CoM vibration frequency $\tilde{\omega }$ can differ from the ω frequency of the HO basis [187]. As illustrated in Fig. 9, the CoM effect remains negligible for these low-lying states.

Fig. 9Full-screen View

(Color online) ²⁴O spectra calculated using NNLO_opt and NNLO_sat interactions. The first two columns display the results from real-energy EOM-IMSRG calculations (denoted as R-IMSRG) without and with CoM treatment $\beta H_{\text{CoM}}$, using the multiplier β = 5. The subsequent three columns present the EOM-G-IMSRG calculation (denoted as G-IMSRG), which are compared with the data from [188, 189]. Resonant states are highlighted by shading, and their widths (in MeV) are annotated nearby. The gray shading indicates the continuum region above the particle emission threshold. The figure is taken from [57]

As described in Sec. 2.4, the equation-of-motion approach can be used to calculate nuclear excited states. Fig. 9 presents the calculated spectrum of ²⁴O, showing resonant excited states. The EOM Gamow IMSRG (indicated by EOM-G-IMSRG) calculations reproduced the experimental excitation energies and resonance widths of the observed states well. A high excitation energy of the first 2⁺ state supports the shell closure at N = 16 in the oxygen chain. Additionally, the calculation predicted three resonant states around the excitation energies of 8 MeV with $J^{\pi }{=2}^{+}-4^{+}$, aligning with the experimentally ambiguous states observed around 7.6 MeV [189]. This finding is consistent with the complex coupled cluster calculation which uses a schematic 3NF [18].

The Borromean halo nucleus ²²C poses significant challenges for many theoretical models [58, 59, 60]. Our GHF calculations suggested that the neutron $\nu 1s_{1/2}$ orbital is weakly bound. Therefore, the two-neutron configuration ${\left(\nu 1s_{1/2}\right)}^2$ is responsible for the formation of the halo structure [181-184]. Figure 10 presents the ground-state densities of ²²C calculated by the real-energy R-IMSRG and complex-energy G-IMSRG using two different chiral interactions. The density was computed by an effective density operator evolving within the Magnus framework (16, 17). The G-IMSRG calculation revealed a long tail in the density distribution, supporting the halo structure of ²²C.

Fig. 10Full-screen View

²²C ground-state densities calculated by real-energy IMSRG (R-IMSRG) and Gamow IMSRG (G-IMSRG), displayed on a logarithmic scale. The inset provides a detailed view of the densities in the central region of the nucleus on a standard scale. The figure is taken from [57]

To assess the continuum effect of the s channel on the properties of ²²C, we performed two types of G-IMSRG calculations using: (i) discrete s states obtained from the real-energy HF calculation, and (ii) Berggren s states obtained from the complex-energy GHF calculation. In both calculations, the neutron d_3/2 channel remained in the GHF basis. Calculations using NNLO_opt with discrete real-energy HF s states yielded a radius of 2.798 fm for the ²²C ground state, which increased to 2.928 fm upon incorporating the continuum s wave. Similarly, the calculations using NNLO_sat yielded a radius of 2.983 fm for the real-energy discrete s states and 3.139 fm for the continuum GHF s wave. The experimentally estimated radius was reported to be 5.4±0.9 fm in an earlier work [181], and later works found it to be 3.44±0.08 fm [184] and 3.38±0.10 fm [190]. These findings highlight the crucial role of the continuum s wave in predicting the radius and understanding the halo structure.

Currently, no experimental data are available for the excited states in ²²C. Figure 11 displays the EOM-G-IMSRG predictions of low-lying states, benchmarked against results from complex CC calculations. Both methods yielded consistent results. The first 2⁺ excited state was bound in both G-IMSRG and coupled cluster calculations. We found that the $2_1^{+}$ state was dominated by the proton 1p1h excitation from the proton 0p3/2 hole to proton 0p1/2 particle orbits. The proton $2_1^{+}$ excited state was lower in energy than the neutron 2⁺ state calculated by the real-energy SM with the ¹⁴C core [191, 192]. The real-energy R-IMSRG results in Fig. 11 show a neutron 2⁺ energy similar to that in [191, 192]. Additionally, there were superposed resonant states with $J^{\pi }{=1}^{+}-4^{+}$ at energies of 3.5-4.0 MeV and widths of 0.15-0.25 MeV. The NNLO_sat calculations yielded slightly higher excitation energies and broader resonant widths for these states, as illustrated in Fig. 11. The resonances were primarily characterized by neutron 1p1h excitations from the v0d_5/2 hole to v0d_3/2 particle orbitals. Structure and decay modes of loosely-bound nuclei are of interest in many respects [193, 194].

Fig. 11Full-screen View

Excited states in ²²C predicted by R-IMSRG and G-IMSRG using two different chiral interactions compared with complex coupled cluster results. The channels listed at the top of the panel indicate that the partial waves are treated in the resonance and continuum Berggren representation. The other labels are the same as those in Fig. 9. The figure is taken from [57]

Few-body decay by GCC

In this section, the discussion centers on the decay properties of open quantum systems, with a particular focus on two-proton (2p) emitters, exploring how deformation and continuum effects influence the decay dynamics of these systems. Previous research highlighted the significant role that these factors play in shaping the decay characteristics of 2p emitters. A critical aspect of our analysis involves extracting structural information from these systems through the measurement of asymptotic nucleon-nucleon correlations, which are experimentally accessible.

In analyzing these correlations, the objective is to deepen our understanding of the universal properties inherent to open quantum systems. This approach not only elucidates the fundamental interactions within these systems but also provides a framework for interpreting experimental results in terms of underlying nuclear structure and dynamics. Such insights are invaluable for advancing our comprehension of the complex behaviors exhibited by open quantum systems under various conditions.

3.4.1

Impact of structure on the decay process

As the heaviest 2p emitter identified to date, ⁶⁷Kr serves as a pivotal case study for examining the influence of structure on decay properties. Notably, shape coexistence is often observed with Kr isotopes. In addition, the deformation effects can significantly impact the lifetime of a decaying system, as evidenced by prior research on one-proton (1p) emitters [135-137, 144, 195-199]. This offers a good opportunity to study how the 2p decay properties change as a function of deformation.

As discussed in [68], Fig. 12a illustrates the proton Nilsson levels (labeled by asymptotic quantum numbers Ω[${\text{Nn}}_z\text{Λ}$]) within the core-p potential. At modest deformations, specifically $\left|{\beta }_2\right|\le 0.1$, the valence protons predominantly occupy the f_5/2 shell. The half-life predicted under vibrational conditions, calculated as T_1/2>218 ms, significantly surpasses the experimentally observed value by more than an order of magnitude, as depicted in Fig. 12b. This discrepancy underscores the need for further theoretical refinement and is consistent with prior theoretical estimates [200, 201], suggesting an intricate interplay between nuclear structure and decay dynamics in ⁶⁷Kr.

Fig. 12Full-screen View

Top: Nilsson levels Ω[${\text{Nn}}_z\text{Λ}$] of the deformed core-p potential as a function of the oblate quadrupole deformation β₂ of the core. The dotted line indicates the valence level primarily occupied by the two valence protons. Bottom: Decay width (half-life) of the 2p ground state radioactivity of ⁶⁷Kr. The solid and dashed lines mark the results for the rotational and vibrational couplings, respectively. The rotational-coupling calculations were carried out by assuming that the 1/2[321] orbital is either occupied by the core (9/2[404]-valence) or valence (1/2[321]-valence) protons. The figure is taken from [68]

As the core deformation increases, a notable oblate gap at Z=36 emerges due to the descending 9/2[404] Nilsson level, which stems from the 0g_9/2 shell. This gap plays a crucial role in shaping the oblate ground state (g.s.) configurations of proton-deficient Kr isotopes [202-204]. As the oblate deformation intensifies, the structure of the valence proton orbital transitions from the 9/2[404] ($\mathcal{l}=4$) state to the 1/2[321] orbital, featuring a significant $\mathcal{l}=1$ component. The precise level crossing between 1/2[321] and 9/2[404] is contingent on the specifics of the core-proton parametrization, yet the overarching pattern remains consistent as depicted in Fig. 12a: a shift from the 2p wave function dominated by $\mathcal{l}=4$ components towards $\mathcal{l}=1$ components as the oblate deformation progresses.

Figure 12b illustrates the predicted 2p decay width within the rotational model for two scenarios: (i) the 1/2[321] level is integrated into the core with the valence protons predominantly residing in the 9/2[404] level, and (ii) the valence protons primarily occupy the 1/2[321] level. Consequently, the reduction in $\mathcal{l}$ content within the 2p wave function markedly enhances the decay width.

At a deformation of ${\beta }_2\approx -0.3$, aligned with estimates from mirror nuclei and various theoretical models [205-208], the calculated 2p ground state half-life of ⁶⁷Kr is ${24}_{-7}^{+10}$ ms. This estimation not only concurs with experimental findings [209] but also underscores the significant impact of nuclear shape and structural dynamics on the decay properties of ⁶⁷Kr.

3.4.2

Dynamics of two-proton decay

For the light 2p emitters, both direct and sequential decays are possible [210, 211], providing a good opportunity to study the decay dynamics as well as the interplay between nucleon-nucleon correlation and single-particle emission.

An illustrative example is the decay of ⁶Be, where the neighboring g.s. of ⁵Li exists as a broad resonance, characterized by a proton decay width of Γ=1.23 MeV [212]. The complex three-body decay dynamics of ⁶Be remain an area of active research and debate, with existing studies indicating unresolved aspects of the decay [141, 211, 213-219]. Theoretical predictions of a diproton structure and experimental observations of broad angular correlations between emitted protons lead to conflicting interpretations.

Based on the time-dependent approach, Fig. 13a,b depicts the temporal evolution of the 2p density and momentum distribution in the ground state of ⁶Be across an extensive time frame. Initially, at t=0, the wave function inside the nucleus shows a localized character with a density distribution exhibiting two maxima. These maxima represent diproton (compact) and cigar-like (extended) configurations based on the relative distances between the valence protons [220].

Fig. 13Full-screen View

The density and momentum distributions of two-nucleon decays from the g.s. of ⁶Be (left) and ⁶He′ (right) for four different time slices. The density distributions are shown in the Jacobi-T coordinates (see Fig. 3). The momentum distribution of the second nucleon is depicted with respect to the momentum of the first nucleon. To clearly show the asymptotic wave function, all the particle densities (in fm^-1) are multiplied by the polar Jacobi coordinate ρ. The dimensionless momentum (angular) distributions are divided by the total momentum k. The figure is taken from [70]

As the decay progresses, Fig. 13b highlights two predominant flux branches. The primary branch shows protons emitted at narrow angles, indicating the presence of a diproton structure during the tunneling phase. This phenomenon is interpreted through nucleonic pairing, which favors low angular momentum states, reducing the centrifugal barrier and enhancing the 2p tunneling likelihood [216, 149, 150, 220]. The secondary branch illustrates protons emitted in nearly opposite directions. Despite their spatial separation, these protons exhibit simultaneous decay, suggesting three-body decay dynamics characterized by correlated decay pathways of the protons with respect to the core. This configuration reveals intricate interplays within the decay process, shedding light on the multifaceted nature of three-body decays in light 2p emitters.

After tunneling through the Coulomb barrier, the two emitted protons tend to gradually separate due to Coulomb repulsion. This is reflected in the bent trajectory of the diproton decay branch and gradual reduction of the momentum alignment observed in Fig. 13a,b. Eventually, the 2p density becomes spatially diffuse, which is consistent with the broad angular distribution measured in [217]. One may notice that even beyond 100 fm (at t≈2 pm/c), the Coulomb repulsion tends to reduce the inter-proton correlation. According to our calculations, the angular correlation starts to stabilize only after very long times greater than 9 pm/c. Therefore, to obtain meaningful estimates of asymptotic observables, very long propagation times are required.

After the protons tunnel through the Coulomb barrier, the Coulomb repulsion between them leads to an increasing spatial separation, impacting their trajectories and momentum alignment. This dynamic is depicted in the bending trajectory of the diproton decay branch, and a reduction in momentum alignment is observed in Fig. 13a,b. Over time, the 2p density distribution becomes increasingly diffuse, aligning with the broad angular distributions reported in experimental studies such as [217].

To provide insights into the nuclear-Coulomb interplay in two-nucleon decay processes, an artificially unbound variant of ⁶He, denoted as ⁶He′, was built to study its two-neutron decay. Initially, the density distributions of ⁶He′ and ⁶Be display similarities, as depicted in Fig. 13, a consequence of the isospin symmetry inherent in the nuclear force.

However, the absence of Coulomb repulsion in the ⁶He′ scenario significantly influences the decay dynamics. In the case of ⁶He′, the dineutron decay branch is more pronounced, with the emitted neutrons maintaining their spatial correlations over time more robustly than their proton counterparts in ⁶Be. This differential behavior leads to distinct nucleon-nucleon correlation patterns in the asymptotic regime. As a result, the nucleon-nucleon correlations observed in ⁶He′, as illustrated in Fig. 14, diverge markedly from those in ⁶Be, underscoring the critical role of Coulomb forces in shaping the decay pathways and final-state interactions in these mirror nuclei.

Fig. 14Full-screen View

Asymptotic energy (left) and angular (right) correlations of emitted nucleons from the g.s. of ⁶Be (top) and ⁶He′ (bottom) calculated at t=15 pm/c with different strengths of the Minnesota interaction [221]: standard (solid line), strong (increased by 50%; dashed line), and weak (decreased by 50%; dash-dotted line). Also shown are the benchmarking results obtained using the Green’s function method (GF; dotted line) with standard interaction strength; θk is the opening angle between ${\overrightarrow{k}}_x$ and ${\overrightarrow{k}}_1$ in the Jacobi-Y coordinate system, and E_pp/nn is the kinetic energy of the relative motion of the emitted nucleons (see Fig. 3 for definitions). The figure is taken from [70]

3.4.3

Nucleon-nucleon correlations

The experimental measurements of nucleon-nucleon correlations among emitted nucleons provide data pertinent to the nuclear structure. This methodology serves as a distinctive avenue for investigating the internal structure and decay mechanisms of 2p emitters. However, the nucleon-nucleon correlations recorded by detectors are influenced by final-state interactions and distort the original nuclear correlations. Consequently, self-consistent theoretical frameworks are urgently required to establish a linkage between the nuclear structure and observed asymptotic correlations. In this section, we elucidate how the GCC method and time-dependent (TD) approach have been effectively employed to address such challenges.

Compared to the neutron dripline, the proton dripline is relatively closer to the line of β-stability, facilitating the acquisition of 2p correlation data in several instances, as evidenced by the findings reported in studies [222, 223]. Remarkably, the energy correlation of protons emitted from the ground state of ¹²O closely resembles those observed in other sd-shell 2p emitters, such as ¹⁶Ne and ¹⁹Mg. This resemblance suggests potential structural similarities in the configurations of valence protons across these isotopes despite their differing proton numbers. This observation underscores the intricate interplay between nuclear structure and decay processes, hinting at underlying uniformities in the spatial and energetic distributions of valence protons within this specific shell.

To elucidate the 2p correlation patterns observed in ¹²O and its isotonic neighbor ¹¹O, a time-dependent approach was employed, as detailed in recent research [71]. The initial 2p density configurations in ^11,12O exhibit a prominent diproton arrangement alongside a secondary, cigar-like structure [224, 220], bearing resemblance to configurations typical of p-shell nuclei. However, in the case of ¹²O, the protons emerging from these configurations coalesce, leading to a broad distribution as reported in [71]. This pattern contrasts starkly with the decay dynamics of ⁶Be, as shown in Fig. 13 and [70]. These observations align with the flux current calculations presented in [220], which indicate a competitive interplay between diproton and cigar-like decay modes, culminating in a so-called democratic decay process. This comparative analysis underscores the complex interdependencies within nuclear decay pathways and highlights the distinctive decay characteristics of ¹²O relative to its nuclear peers.

Consequently, the asymptotic 2p correlations for ¹²O, as depicted in Fig. 15, demonstrate robust alignment with the empirical data [223] and corroborate prior theoretical investigations [225]. Notably, the subtle discrepancies between the experimental results and theoretical predictions, enhanced via Monte Carlo simulations incorporating experimental acceptance and resolution, might be mitigated by refining the Minnesota force initially utilized for modeling the nucleon-nucleon interactions. More importantly, the absence of distinct diproton emissions during the temporal evolution suggests that a low-E_pp correlation does not inherently imply a diproton decay. Moreover, recognizing that the 2p system may manifest as a subthreshold resonance characterized by a broad decay width around 1 MeV [226], is critical. This continuum characteristic potentially influences the energy correlation observed in 2p emitters possessing minimal decay energies, underscoring the complex interplay between nuclear structure and decay dynamics.

Fig. 15Full-screen View

Asymptotic (a) energy and (b) angular correlations of the protons emitted from the two-proton unbound ¹²O isotope. Also shown is (c) the momentum scheme for three-body system. Theoretical distributions were obtained within the time-dependent approach (TD, red line) at t=15 pm/c. MC (green line chart) labels the Monte Carlo simulation of TD results which include experimental resolution and efficiency [223]. The calculated 2p correlations (TD and MC) are compared with experimental data (Exp, blue histogram) of [223]. A is the mass number and k₁, k₂, and kc are the momenta of the nucleons n₁ and n₂, and the core c, respectively in the c.m. coordinate frame. The figure is taken from [71]

To further illustrate, the lightest oxygen isotope, ¹¹O, has been characterized as a broad structure encompassing multiple resonances [224, 227-230]. Utilizing GCC calculations, four low-lying states with quantum numbers $J^{\pi }{=3/2}_1^{-}$, ${5/2}_1^{+}$, ${3/2}_2^{-}$, and ${5/2}_2^{+}$ were predicted within the experimental energy range [220]. Each state exhibits a substantial decay width ranging from 1 MeV to 2 MeV. In the time-dependent calculations [71], these states were propagated individually, effectively disregarding potential interference effects. This approach underscores their significant decay widths, which foster robust continuum coupling, resulting in a more homogeneous density distribution throughout the decay process compared to ¹²O, as noted in [71].

The emergent Y-type correlations exhibit a pronounced dependency on angular momentum, which can be instrumental in experimentally determining spin assignments. To simulate the asymptotic correlations of the emitted valence protons, the correlations of these four states were amalgamated, utilizing the weights derived from resonance-shape fitting [224]. This composite correlation, depicted in Fig. 16(a), closely aligns with the experimental observations shown in Fig. 16(b), lending credence to the hypothesis that the observed broad structure integrates states with $J^{\pi }{=3/2}^{-}$ and 5/2⁺ [224, 227]. Moreover, Figs. 16(g) and (h) illustrate the energy and angular correlations, respectively, adjusted for experimental resolution and efficiency. The qualitative agreement between these corrected simulations and experimental data provides additional support for the multi-resonance composition within the observed peak, further validating the complex resonance structure of ¹¹O.

Fig. 16Full-screen View

(a) Theoretical and (b) experimental Jacobi-Y correlations of the two protons emitted from the broad low-energy structure in ¹¹O; (c)-(f) the corresponding contributions predicted from each low-lying state. The experimental resolution and efficiency have been considered in (a) Monte Carlo simulations. Also shown are the corresponding (g) energy and (h) angular correlations obtained in the time-dependent (TD) calculations (red line), Monte Carlo (MC) simulations (green step chart), and experiments (Exp, blue histogram). The figure is taken from [71]

3.4.4

Exotic decays in open quantum systems

In addition to 2p decay, there exist several other exotic decay modes in nuclear physics, such as two-neutron decay, two-alpha decay, and multi-particle emission. These decay processes exhibit intriguing behaviors that are unique to the quantum realm. The classical understanding of radioactive decay is based on the principle that the rate of decay is directly proportional to the quantity of the radioactive substance present. This model fundamentally assumes that the decay is a stochastic process at the level of individual particles, meaning that the probability of decay is independent of the system’s previous history.

However, this classical perspective is challenged within the quantum framework due to phenomena such as memory effects and quantum entanglement. These effects introduce deviations from the exponential decay law, traditionally used to describe radioactive decay. Notably, quantum mechanics reveals that decay probabilities can exhibit non-exponential behavior at very short and very long timescales. Theoretical and experimental studies have demonstrated that quantum systems do not always follow the expected exponential decay pattern [151, 231-241].

These findings underscore the complex nature of decay processes in quantum systems, where the inherent properties of the particles and their interactions can lead to observable departures from classical predictions. The implications of these quantum behaviors are profound, impacting our understanding of fundamental decay processes and the predictive models used in nuclear physics.

In analyzing the non-exponential decay in quantum systems, one crucial concept is the survival amplitude $\mathcal{A}(t)$. This amplitude is defined as the overlap between the initial quantum state Ψ(0) and the state at a later time Ψ(t). Mathematically, the survival amplitude can be expressed and computed using the Fourier transform of the spectral function ρ(E), as illustrated in [242]: $\mathcal{A}(t)=\langle \text{Ψ}(0)|\text{Ψ}(t)\rangle ={\displaystyle {\int }_0^{+\infty }}\rho (E)e^{-i\frac{E}{\hslash }t}\text{d}E\mathrm{.}$ (60) In this formulation, ρ(E), which is the probability density of finding a system with energy E, is derived from the squared modulus of the projection of the state vector Ψ onto the real-energy eigenstates $\langle E|\text{Ψ}\rangle$. Thus, $\rho (E)={\left|\langle E|\text{Ψ}\rangle \right|}^2$ represents the distribution of the initial state over the various energy eigenstates.

The survival probability $\mathcal{S}(t)$, which is a measure of the likelihood that the system remains in its initial state at time t, is then calculated in a straightforward manner from the survival amplitude: $\mathcal{S}(t)=|\mathcal{A}(t{)|}^2\mathrm{.}$ (61) This probabilistic measure reflects how the state evolves over time, deviating from its initial configuration. This deviation is a key indicator of quantum mechanical effects in decay processes and provides insight into the complex nature of quantum dynamics.

Non-exponential decay of a threshold resonance. The survival probability of a quantum state is intrinsically linked to the energy distribution described by the spectral function of the system. Typically, for a system exhibiting exponential decay, the spectral function is expected to follow a Breit-Wigner distribution, characterized by a Lorentzian shape centered around the resonance energy with a width corresponding to the decay rate. However, this does not hold for near-threshold states, particularly those with large decay widths [233, 151, 240, 241, 243].

These near-threshold states often display significant deviations from the exponential decay law. The temporal evolution of resonance states has been shown to involve both exponential and non-exponential components [244, 245]. The exponential components, characterized by a rapid decrease in probability, dominate the early time behavior of the decay process. However, as these components decay, the non-exponential elements become more prominent.

Over time, as the influence of the exponential decay wanes, a transition to a power-law regime becomes inevitable. This regime is indicative of the long-time tail behavior common to quantum systems with broad spectral distributions. Such transitions are crucial for understanding the full dynamics of decay processes, particularly in scenarios where classical exponential decay laws fail to capture the complexities introduced by quantum mechanical principles.

Figure 17a,c illustrates this transition, highlighting how the decay initially follows an exponential decrease before transitioning to a power-law decay. This behavior underscores the complex nature of quantum decay processes and the need for a deeper exploration into the underlying physics, particularly for states near the energy threshold or those with significant quantum mechanical interactions.

Fig. 17Full-screen View

Survival probability $\mathcal{S}(t)$ as a function of time (relative to T_1/2) for (a) the 1/2^- state of ⁹He for different depths V₀ of the Woods-Saxon (WS) potential, and (c) the low-lying states of ⁹N. The near-threshold behavior of the spectral function ρ (relative to the Breit-Wigner distribution) is shown for (b) neutron and (d) proton s,p, d partial waves. The polar angle φ indicates the location of the resonant state in the complex-k plane. Also shown is the survival probability for the virtual 1/2⁺ state in ⁹He. For this state, T_1/2 was assumed to be 20 fm/c. The figure is taken from [246]

The universality of the transition from exponential to non-exponential decay is a key characteristic of quantum decay processes, the specific dynamics of which are deeply influenced by several factors. These include the structure of the initial state, the chosen decay channel, and notably, the nature of the scattering continuum which drives the post-exponential decay behavior [246].

To provide a concrete example of these dynamics, the survival probability of the 1/2^- resonant state in ⁹He has been analyzed by adjusting the depth of core-nucleon potential used in the calculations [246]. This adjustment affects the resonant states, which can be characterized in the complex-k plane. The positioning of these states is determined using the polar angle $\phi =-{{\displaystyle \cot }}^{-1}(2E/\text{Γ})/2$, which offers insights into the relative contributions of the exponential and non-exponential decay components in state evolution.

In this analysis, as depicted in Fig. 17a, the deviations from exponential decay become increasingly pronounced as φ shifts toward -45°. This angular movement suggests a strengthening of the non-exponential decay component, particularly in threshold resonances where the resonant energy Er is approximately equal to the decay width Γ. This result is consistent with [240, 247, 248], which indicate that post-exponential decay features tend to dominate more rapidly in systems where the resonance lies near the decay threshold. Such resonances provide a clearer and more readily observable transition to non-exponential decay, making them ideal subjects for experimental and theoretical studies aiming to explore quantum decay dynamics beyond the conventional exponential model.

Interference between near-lying states. Besides the threshold effect, the decay dynamics of quantum systems can also be significantly influenced by the interaction between closely lying resonances, particularly when these states share the same spin-parity configuration. This scenario leads to an intricate interplay due to the interference between overlapping resonances, which in turn modify the decay characteristics [249-251].

The mechanism underlying this behavior is related to the quantum interference effects, whereby the wavefunctions of the resonant states overlap and coherently interact with the continuum states. This interplay can lead to a redistribution of decay widths among the resonances, with one or more states experiencing an enhancement in decay widths due to the increased coupling [252-255]. This enhanced coupling is a critical factor in the non-exponential decay characteristics observed in such systems, as it directly impacts the decay pathways and probabilities.

To provide a detailed examination of how continuum coupling affects the spectral functions of overlapping resonances, a hypothetical study was conducted on a two-level 0⁺ system in an artificial nucleus labeled as ⁶He′ [246]. This recent study [246] explored how two 0⁺ states, lying close in energy, interact with each other and the continuum, highlighting that the interference between these states not only affects their individual decay rates but also alters the overall spectral shape of the system. This interaction leads to one of the resonances showing a collective enhancement in decay width, which is a direct manifestation of the increased coupling with the continuum.

In this scenario, the excited state $|1\rangle$ predominantly features a d² configuration, whereas the ground state $|2\rangle$ mainly consists of a p² configuration. Figure 18 illustrates the evolution of the spectral functions and the corresponding survival probabilities for different energy splittings, $\Delta E=\left|E_r\left(1\right)-E_r\left(2\right)\right|$, of the doublet states.

Fig. 18Full-screen View

Interference between two close-lying 0⁺ resonances in ⁶He′ for the three values of the doublet Δ E (in MeV) energy splitting. Left: Spectral functions versus decay energy. The arrow indicates the suppression of the spectral function of $|2\rangle$. Right: Time dependence of the corresponding survival probabilities. The decay widths (in keV) of the doublet (Γ₁, Γ₂) are (34, 60), (30, 52), and (6, 68) for large, moderate, and small values of Δ E, respectively. The figure is taken from [246]

When the energy splitting Δ E is large, only a minor suppression occurs at the tail of the spectral function for state $|2\rangle$, and both states exhibit comparable decay widths. However, as the states begin to overlap, significant interference effects arise, which dramatically impact the spectral functions of the doublet (Fig. 18e, f). This interference leads to pronounced deviations from the exponential decay regime in survival probabilities.

Specifically, state $|1\rangle$ decays much more rapidly than its intrinsic decay width would suggest, whereas state $|2\rangle$ exhibits a remarkably slow decay. These observations are consistent with the findings of previous studies [253, 255], which discuss how such exponential deviations during the decay process can occur between any near-lying resonances of the same spin-parity. This phenomenon is driven by virtual transitions governed by the scattering continuum and differing orbital angular momentum structures of the doublet states.

This behavior highlights the complex dynamics that can arise from interference effects in quantum decay processes, particularly when closely spaced resonances are involved. The interplay between the initial state configurations, decay channels, and the nature of the continuum coupling leads to non-trivial modifications in decay rates and survival probabilities, underscoring the intricate nature of quantum mechanical decay processes.