Neutron-rich nuclei are of great importance for understanding the astrophysical r-process [1]. For instance, those around the closed neutron shell N = 126 can provide information on the solar r-abundance distribution [2]. The neutron-rich radioactive isotopes can also be used as projectiles for synthesizing super-heavy nuclei which is one of the most interesting challenges in nuclear physics. However, new heavy neutron-rich nuclei out of limits of the present nuclear chart are hardly to be produced via traditional ways, such as fusion reactions with stable beams at low energies or fragmentation of heavy projectile at relativistic energies.

In recent years, theoretical predictions indicate that multi-nucleon transfer (MNT) reactions would be a possible route to produce heavy neutron-rich nuclei far away from the stability line [3, 4]. The experimental results of ¹³⁶Xe + ¹⁹⁸Pt at the incident energy of 8 MeV/nucleon [2] show that the production cross sections of neutron-rich nuclei with N = 126 are orders of magnitude larger than those obtained in fragmentation reaction of ²⁰⁸Pb (1 AGeV) + Be [5]. But in some other experiments of MNT reactions with stable beams such as ⁶⁴Ni + ²⁰⁷Pb [6] and ¹³⁶Xe + ²⁰⁸Pb [7, 8], no new neutron-rich nuclei is detected. Based on the N/Z equilibrium concept, theoretical predictions show that using neutron-rich radioactive beams can remarkably improve the production cross sections of neutron-rich isotopes along N = 126 [9, 10]. However, due to the fact that the experimental intensities of radioactive beams are orders of magnitude lower than stable beams, the advantage of using radioactive beams may be canceled. Opportunities will arise in the near future on the second generation radioactive beam facilities like the European EURISOL [11]. At the present time, to find an optimum stable projectile-target combination are highly appealed.

To describe MNT reactions, many theoretical models are developed. For example, the multidimensional Langevin model shows great success on predicting the isotopic production cross sections, even for those nuclei far away from projectile and target [3, 4, 12]. Semi-classical GRAZING model [13, 14] and the complex Wentzel-Kramers-Brillouin (CWKB) model [15] can well describe transfer process in peripheral collisions. In GRAZING model, the reactants move on classical trajectories in the combined field of Coulomb repulsion and nuclear surface-surface attraction. Surface modes of the colliding nuclei are taken into account. Independent single nucleon transfer between the projectile-like and target-like nuclei during the collision is governed by the quantum coupled equations [13, 14]. Neutron evaporation is considered for the excited primary fragments. The mass, charge, energy, and angular momentum distributions of the products produced in grazing collisions can be obtained. For details of GRAZING model and its applications (see, e.g., [16–19], and references therein). The dinuclear system model can reproduce the experimental data related with quasi-fission or deep-inelastic collisions [10, 20–23]. The improved quantum molecular dynamics (ImQMD) model [24–26] is capable of describing collisions from central to very peripheral regions on a microscopic basis, the widths of isotopic distributions can be reproduced in ImQMD since stochastic two-body collisions are taken into account [27–29]. The time-dependent Hartree-Fock (TDHF) theory shows success in describing few-nucleon transfer process [30–34]. Recently, the stochastic mean-field (SMF) approach beyond TDHF [35] has been proposed to investigate the damped MNT reaction ¹³⁶Xe + ²⁰⁸Pb, the experimental broad mass distribution can be reproduced.

In this work, the GRAZING model and TDHF theory incorporating with GEMINI++¹ [36] are adopted to investigate the production of neutron-rich nuclei in MNT reactions. The TDHF theory [37] is based on the independent particle picture and is a good approximation to the nuclear many-body problem. It is capable of describing low-energy heavy-ion reactions and provides insight on the average behavior of the dynamics. The state-of-art TDHF calculations are performed in a three-dimensional (3D) framework without any symmetry constraints due to the advances in computational power. It has been applied for investigations on various subjects, for instance, collective vibration [38, 39], fusion reaction [40, 41], fission dynamics [42–44], dissipation mechanism [45–48]. Recently, 3D TDHF is applied to MNT reactions [30–33, 49]. Fluctuation and dissipation can not be properly described in TDHF since two-body collisions and internucleon correlations are not included. Thus, widths of the mass or isotopic distributions in MNT reactions at incident energies above the Coulomb barrier are underestimated in TDHF. To consider these effects in TDHF is beyond the scope of this paper and further studied will be carried out in the future.

This paper is outlined as follows. In section 2, the TDHF approach and particle-number projection (PNP) method as well as the numerical details are introduced. In section 3, numerical results of the isotopic production cross sections in both GRAZING and TDHF approach are shown. Dynamical properties of the reactions are also discussed for TDHF. Finally, a summary is drawn in section 4.

In this section, we briefly introduce TDHF formalism, PNP method and the coupling to GEMINI++. Details of the TDHF theory can be found in, e.g., [50, 51], and references therein. In TDHF approach, time evolution for the single-particle states are described by a set of coupled non-linear equations

where ψ_{α, ν}(r, t) is the single-particle state and N is the total number of states. ĥ is the self-consistent mean-field Hamiltonian of single-particle motion and it is always related to Skyrme energy density functional (EDF) which depends on local densities [52]. We underline that there are no adjustable parameters in the TDHF approach. The uncertainty in TDHF calculations may arise from the uncertainty of fundamental nuclear properties, such as the distribution of shape deformation of the reactants in their ground states. This can be obtained by preparing reactants with deformation constraints. The above TDHF equations can be derived from the variation of the action [53]

where

is the correlated many-body wave function of the system and is given by a single Slater determinate.

In the present work, the 3D unrestricted TDHF code Sky3D [51] with Skyrme SLy5 parametrization [54] is adopted for both the static and dynamic calculations. The static HF calculations are performed on 32 × 32 × 32 Cartesian grids with 1.0 fm grid spacing in all three directions. In dynamical calculations, the meshes are extended to 70 × 32 × 70 with the same grid spacing in static HF. The projectile and target are initially placed at a separation distance of 24 fm and then boosted with the associated center-of-mass energy E_c.m. and the impact parameter b. The time step △t is set to be 0.2 fm/c and six-order Taylor series expansion is employed. The TDHF simulations are stopped when the separation distance of the primary fragments after collision reaches 30 fm. Since single-particle wave functions are partially exchanged between the projectile-like and target-like nuclei in the collision process, the outgoing states are not the eigenstates of the particle-number operators (Ẑ for protons and N^ for neutrons) but superpositions of them. One can only get the expectation (mean) values of the charge and mass numbers for each fragment after collision. If one wants to get the distributions of proton or neutron numbers in one of the primary fragments, one should project the many-body states on good particle numbers by introducing the PNP operator [50, 55]

where q = n, p labels the nucleon species, and Nⁿ means N neutrons while N^p means Z protons. The subscript V denotes the region of coordinate space encompassing one of the primary fragments, N^Vq=∑α=1NqΘV(r) and Θ_V(r) = 1 if r ∈ subspace V and 0 elsewhere. The integral is performed with an M−point uniform discretion. For convergence, M is set to be 300.

The probability to find N^q particles in subspace V is then obtained accordingly. The cross section of primary fragment with neutron number N and charge number Z at a certain incident energy is

where P(N, Z; b) = P(N; b)P(Z; b) represents the probability to find N neutrons and Z protons at the impact parameter b; b_max is a cutoff impact parameter which depends on the incident energy and should be large enough to guarantee that most of the transfer cross sections are included. But it is not necessary to set b_max to be too large, because the transfer probability is extremely small and elastic collision dominates in very peripheral collisions. In this work, we set it to be 10 fm for ¹³⁶Xe+²⁰⁸Pb at E_c.m. = 450 MeV and 13 fm for ¹³⁶Xe+¹⁹⁴Ir at E_c.m. = 720 MeV. b ranges from 0 to b_max with the interval Δb = 1 fm.

Deexcitation of the primary fragments are considered by using the statistical code GEMINI++ with default parameters [36, 56]. All possible sequential binary decay modes, from emission of nucleons and light particles through asymmetric to symmetric fission as well as the γ-emission are included in GEMINI++. The code needs information of a primary fragment including charge and mass number as well as the excitation energies and the angular momentum as the inputs. Detailed calculations of these quantities can be found in Jiang and Wang [34]. The deexcitation of a certain projectile-like fragment (PLF) or target-like fragment (TLF) should be repeated M_trial times due to the statistical nature of GEMINI++. Here M_trial = 1, 000 is used. The number of events in which final fragment with (N_final, Z_final) is counted and denoted as M(N_final, Z_final). Then the final production cross section is given as

Owing to the intrinsic stochastic nature of the Monte Carlo method employed in GEMINI++, Type A standard uncertainties for the isotopic production cross sections are calculated in the simplest case. The deexcitation process of a certain fragment is performed ten times repeatedly with M_trial = 1, 000 for each time. Average values of σ(N_final, Z_final) and the uncertainties can be obtained straightforwardly. We find the uncertainties are very small. For the sake of simplicity, the deexcitation process is performed only once with M_trial = 1, 000 in this work.

¹³⁶Xe + ²⁰⁸Pb is a candidate reaction for the production of neutron-rich nuclei at N = 126. The experiment at E_c.m. = 450 MeV was performed at Argonne in 2015 [8]. In Figure 1 we plot the calculated isotopic production cross sections of the PLFs in this reaction. The results of GRAZING are shown as blue solid lines while those of TDHF+GEMINI are presented as black dashed lines. The experimental data are shown for comparison. The production cross sections of the TLFs are also calculated but have already been published in another paper [34]. One can find that for Z = 51 − 53 and Z = 55 − 58, the magnitude of the peak values can be reproduced by TDHF+GEMINI. Whilst those for Z = 52, 53, 55, 56 can be reproduced by GRAZING. For Z = 54, the peak values are overestimated in both models. This is because the results of (quasi)elastic channels are not excluded in our calculations. One can also find that better predictions are obtained for proton pickup channels than stripping channels in both models. As the number of transferred nucleons increases, discrepancies between model predictions and the experimental data get larger. These isotopes far away from the entrance channel may be produced in strongly damped collisions. Such processes can not be well estimated by the two models: two-body collisions are not considered in TDHF [57, 58] while the GRAZING model only takes grazing collisions into account.

However, no new neutron-rich nuclei were detected in ¹³⁶Xe + ²⁰⁸Pb. The experiment results of ¹³⁶Xe+¹⁹⁸Pt at E_c.m. = 645 MeV [2] indicate that this reaction is a better candidate to produce neutron-rich nuclei with N = 126. Calculations of this reaction are performed by using both GRAZING and TDHF+GEMINI. Unfortunately, no neutron-rich nuclei with N = 126 is obtained in TDHF+GEMINI. The production cross sections of those nuclei predicted by GRAZING are given in Figure 2. Detailed discussions about the predictions of TDHF on this reaction will be reported elsewhere. At the end of this section, we will show some results of ¹³⁶Xe + ¹⁹⁴Ir given by TDHF+GEMINI and discuss the limits of TDHF approach on investigating MNT reactions when the incident energy is much above the Coulomb barrier. In the following, we concentrate on the predictions of GRAZING on neutron-rich nuclei with N = 126.

In order to find another optimum projectile-target combination for producing new exotic neutron-rich nuclei with stable reactants, we carry out a systematic study on 15 MNT reactions with projectile around ¹³⁶Xe and target around ¹⁹⁸Pt at various incident energies by using GRAZING. The results indicate that ¹³⁶Xe + ¹⁹⁴Ir is a surrogate reaction to produce neutron-rich nuclei around N = 126. The production cross sections of isotones with N = 126 in ¹³⁶Xe + ¹⁹⁴Ir are compared with those of ¹³⁶Xe + ²⁰⁸Pb and ¹³⁶Xe + ¹⁹⁸Pt. The center-of-mass energy for all the three systems is set to be E_c.m. = 645 MeV (this energy is the same as the experiment of ¹³⁶Xe + ¹⁹⁸Pt [2], and we will show later that it is also an optimum energy for ¹³⁶Xe + ¹⁹⁴Ir). The simulation results are plotted in Figure 2 as open symbols. The experimental data of ¹³⁶Xe + ¹⁹⁸Pt taken from [2] are shown as black solid squares for comparison. One can first see that more N = 126 isotones with charge number Z ⩽ 78 can be produced in both ¹³⁶Xe + ¹⁹⁴Ir and ¹³⁶Xe + ¹⁹⁸Pt rather than ¹³⁶Xe + ²⁰⁸Pb. Particularly, the system ¹³⁶Xe + ¹⁹⁴Ir has huge advantages for producing N = 126 isotones with Z ⩽ 74. Those nuclei are out of the limits of the present nuclear landscape and are of great interest for nuclear and astro-nuclear physics. One can also find that the experimental data are underestimated by GRAZING. This can be understood since only peripheral collisions are treated in GRAZING and those nuclei with large number of nucleon transferred are produced in damped collisions at small impact parameters.

To find an optimum incident energy to produce more neutron-rich nuclei with N = 126 in ¹³⁶Xe + ¹⁹⁴Ir, the production cross sections of N = 126 isotones with Z = 72 − 77 at various incident energies from E_c.m. = 1.07V_B to 1.85V_B (V_B is the Bass barrier [59] and it is around 410 MeV for this reaction) are calculated by using GRAZING. The results are presented in Figure 3. It can be seen that very neutron-rich isotones with Z ⩽ 74 can not be produced if E_lab/A ⩽ 5.5 MeV. The cross sections of all these isotones increase with the increasing incident energy when E_lab/A < 6.5 MeV. However, when E_lab/A is in the range of 7–9 MeV, plateau-like structures are observed for all the curves. This phenomenon indicates that the production cross sections of those nuclei are insensitive to the incident energy if 7 ⩽ E_lab/A ⩽ 9 MeV.

In Figure 4 we show the isotopic production cross sections of the TLFs with Z = 72 − 77 in ¹³⁶Xe + ¹⁹⁸Pt and ¹³⁶Xe + ¹⁹⁴Ir at E_c.m. = 645 MeV (E_lab/A ≈ 8 MeV) as black dashed and red solid lines, respectively. It can be found that for nuclei on the proton-rich side, larger production cross sections are obtained in ¹³⁶Xe + ¹⁹⁴Ir. For isotopes with Z = 76 and 77, more neutron-rich ones are produced in ¹³⁶Xe + ¹⁹⁸Pt. Whilst for Z = 75 the two systems give comparable results for isotopes on the neutron-rich side. As Z decreases, the advantage to produce more neutron-rich nuclei arises in ¹³⁶Xe + ¹⁹⁴Ir. Note that the results of Z = 72 in ¹³⁶Xe + ¹⁹⁸Pt are not given by GRAZING. The above results of the two reactions are difficult to be understood through the N/Z equilibrium process because N/Z is around 1.52 for ¹³⁶Xe and ¹⁹⁴Ir while it is about 1.54 for ¹⁹⁸Pt and ²⁰⁸Pb, respectively. It might be related with the Q_gg value effect, where Q_gg is the ground state-to-ground state Q value. Other quantum effect like shell effect might also play a role. Further studies should be carried out to check these discussions.

Finally, in Figure 5 we show the density contour plots of ¹³⁶Xe + ¹⁹⁴Ir at some special reaction stages in TDHF for two different initial configurations at E_c.m. = 720 MeV (E_lab/A ≈ 9) and b = 6 fm. Isodensities at half the saturation density (ρ₀/2 = 0.08 fm⁻³) is plotted as black solid lines. We should mention that in TDHF the projectile and target are both deformed in their ground states with β₂ = 0.064 and 0.154, respectively. Since the deformation of ¹³⁶Xe is very small, we only take into account the orientation effect of the deformed ¹⁹⁴Ir at the beginning of the dynamical calculations. The two initial configurations are named “tip collision” (the symmetry axis of ¹⁹⁴Ir is set parallel to the bombarding direction: z-axis) and “side collision” (the symmetry axis of ¹⁹⁴Ir is set parallel to x-axis). It can be found that in both configurations, the composite system is strongly elongated before rupture of the neck. The primary fragments produced at the end of the dynamical calculations also have large deformations. Similar results are obtained for TDHF+GEMINI in ¹³⁶Xe + ¹⁹⁸Pt at E_c.m. = 645 MeV. Such large deformation in the exit channel makes the primary fragments have very large excitation energies. This leads to strong evaporation of nucleons in the deexcitation process. So no neutron-rich isotope with N = 126 is observed in TDHF+GEMINI for ¹³⁶Xe + ¹⁹⁴Ir and ¹³⁶Xe + ¹⁹⁸Pt reactions when the incident energy is much larger than the Coulomb barrier. It is well-known that mean-field model such as TDHF is suitable for low-energy reactions, however, the lack of two-body collisions limits the predictive power of TDHF on estimating the yields of MNT reactions when the incident energy is much above the barrier.

Employing the semi-classic GRAZING model and the microscopic TDHF approach, we have investigated the MNT reactions of ¹³⁶Xe + ¹⁹⁴Ir, ²⁰⁸Pb. Neutron evaporation is considered in GRAZING while GEMINI++ is coupled with TDHF to deal with the deexcitation process. The calculated production cross sections of the PLFs for ¹³⁶Xe + ²⁰⁸Pb at E_c.m. = 450 MeV are compared with the experimental data. The model predictions can well describe the yields of nuclei near the projectile. The predictions of GRAZING on ¹³⁶Xe + ¹⁹⁴Ir show that this reaction has the advantage to produce more neutron-rich nuclei with N = 126 compared with ¹³⁶Xe + ¹⁹⁸Pt. The latter one has been carried out at GANIL [2] and the results are inspiring. No new neutron-rich isotope is observed for these two reactions at energies much above the barrier in TDHF+GMINI which might be interpreted as the lack of two-body collisions in the mean-field theory. Further investigations on ¹³⁶Xe + ¹⁹⁴Ir by using other theoretical models with two-body collisions included, such as ImQMD, are in progress.

All datasets generated for this study are included in the article/supplementary material.

XJ performed numerical calculations. XJ and NW wrote the manuscript.