# Data-Augmentation for Graph Neural Network Learning of the Relaxed Energies
of Unrelaxed Structures

Jason B. Gibson  Department of Materials Science and Engineering, University
of Florida, Gainesville, Florida 32611, USA  Quantum Theory Project,
University of Florida, Gainesville, Florida 32611, USA  Ajinkya C. Hire
Department of Materials Science and Engineering, University of Florida,
Gainesville, Florida 32611, USA  Quantum Theory Project, University of
Florida, Gainesville, Florida 32611, USA  Richard G. Hennig<EMAIL_ADDRESS>Department of Materials Science and Engineering, University of Florida,
Gainesville, Florida 32611, USA  Quantum Theory Project, University of
Florida, Gainesville, Florida 32611, USA

###### Abstract

Computational materials discovery has continually grown in utility over the
past decade due to advances in computing power and crystal structure
prediction algorithms (CSPA). However, the computational cost of the ab initio
calculations required by CSPA limits its utility to small unit cells, reducing
the compositional and structural space the algorithms can explore. Past
studies have bypassed many unneeded ab initio calculations by utilizing
machine learning methods to predict formation energy and determine the
stability of a material. Specifically, graph neural networks display high
fidelity in predicting formation energy. Traditionally graph neural networks
are trained on large data sets of relaxed structures. Unfortunately, the
geometries of unrelaxed candidate structures produced by CSPA often deviate
from the relaxed state, which leads to poor predictions hindering the model’s
ability to filter energetically unfavorable prior to ab initio evaluation.
This work shows that the prediction error on relaxed structures reduces as
training progresses, while the prediction error on unrelaxed structures
increases, suggesting an inverse correlation between relaxed and unrelaxed
structure prediction accuracy. To remedy this behavior, we propose a simple,
physically motivated, computationally cheap perturbation technique that
augments training data to improve predictions on unrelaxed structures
dramatically. On our test set consisting of 623 Nb-Sr-H hydride structures, we
found that training a crystal graph convolutional neural networks, utilizing
our augmentation method, reduced the MAE of formation energy prediction by 66%
compared to training with only relaxed structures. We then show how this error
reduction can accelerates CSPA by improving the model’s ability to filter out
energetically unfavorable structures accurately.

##  I Introduction

The discovery of novel functional materials drives innovation. The process of
discovery has dramatically accelerated over the past decade, partially as a
product of growing crystal structure databases Jain  _et al._ (2013); Kirklin
_et al._ (2015); Draxl and Scheffler (2018); Curtarolo  _et al._ (2012) and
improved computationally based CSPA Oganov  _et al._ (2019), such as genetic
algorithms (GA) Revard  _et al._ (2018), basin hoping Wales and Doye (1997),
elemental substitution Noh  _et al._ (2019), and particle swarm techniques
Chen  _et al._ (2021). CSPA have played a prominent role in successfully
predicting the structure and establishing the stability of high-pressure,
high-temperature superconducting binary hydrides Duan  _et al._ (2014); Liu
_et al._ (2017); Peng  _et al._ (2017). Many recent studies have started
looking for stable ternary hydride superconductors. The addition of a third
element to the binary hydrides can potentially stabilize these materials at
much lower pressureSun  _et al._ (2019); Cataldo  _et al._ (2021); Hilleke and
Zurek (2022). Complex ternary and quaternary materials systems are also
promising candidates for hydrogen storage applicationsHuang  _et al._ (2021).

In particular, GAs have proven their utility to identify thermodynamically
stable phases efficiently; successfully identifying novel materials for
applications such as Li-Ge batteries Tipton  _et al._ (2014) and solar cells
Nguyen  _et al._ (2013). Unfortunately, finding thermodynamically stable
phases in ternary and quaternary systems is notoriously difficult due in part
to the computationally expensive ab initio calculations required to relax and
calculate the energies of GA produced structures, accounting for 99% of the
algorithm’s computational cost. Heiles and Johnston (2013) This places
restrictions on the size and the composition of the unit cells, inhibiting the
exploration of complex material systems.

This computational cost can be reduced by bypassing many of the costly ab
initio calculations via a machine-learned filter or by implementing a more
computationally efficient machine-learned surrogate potential to pre-relax
structures Xie  _et al._ (2021). Wu et al. Wu  _et al._ (2013) fitted a
classical potential to the structures evaluated by density functional theory
(DFT). They used the potential to pre-relax the structures in the GA and only
evaluated the best structures with DFT. Jennings et al. Jennings  _et al._
(2019) used a machine learning (ML) model to predict a structure’s fitness
directly and then only used DFT to evaluate structures that improved the
current population. These methods are still somewhat hindered because many DFT
evaluated structures are required to train a ML model to an adequate fidelity.
Further, the models are specific to the materials’ space the GA is searching,
restricting their application to the given GA search.

Alternatively, there has been work to create universal ML models that
determine a material’s stability by predicting the formation energy of
structures containing elements across the periodic table. Most notably, Xie et
al. Xie and Grossman (2018) predicted formation energy using a crystal graph
convolutional neural network (CGCNN) trained on the materials project (MP)
database Jain  _et al._ (2013). The CGCNN represents a crystal structure as a
multi-graph and builds a graph convolutional neural network on top of the
multi-graph. This enables the model to learn the best features to represent
the structure as opposed to the typical ”handcrafted” feature approach Ahmad
_et al._ (2018) and achieve a formation energy validation MAE of 39 meV/atom
Xie and Grossman (2018). More recently, the MAE of formation energy prediction
of graph-based models continued to decrease to 21-39 meV/atom Choudhary and
DeCost (2021); Cheng  _et al._ (2021); Chen  _et al._ (2019); Park and
Wolverton (2020); Noh  _et al._ (2020).

However, Park et al. Park and Wolverton (2020) found that the model’s
dependence on a structure’s atomic coordinates hinders the model’s predictive
fidelity on structures that strongly deviate from their relaxed states. Given
that to obtain a structure in a relaxed state, a DFT relaxation and hence
energy calculation is needed, the reported MAEs do not represent the model’s
ability to accurately identify unrelaxed structures that would relax to stable
structures. On an unseen test set of 311 stable ThCr2Si2-type compounds, the
CGCNN obtained a reasonable formation energy MAE on relaxed structures of 56
meV/atom. However, the prediction MAE was 370 meV/atom for the same test set
on unrelaxed structures. This high error led to a true positive rate (TPR) of
0.48 when filtering the unrelaxed compounds in the data set Park and Wolverton
(2020).

Noh et al. Noh  _et al._ (2020) directly addressed the high formation energy
prediction MAE of unrelaxed structures by adding two forms of regularization
to the CGCNN in their CGCNN-HD approach. Replacing the softplus activation
function in the convolution function with a hyperbolic tangent and adding
dropout layers Srivastava  _et al._ (2014) between each fully connected layer,
reduced the MAE of the formation energy for a test set of unrelaxed Mg-Mn-O
compounds from 518 meV/atom to 296 meV/atom.

Figure 1: Data augmentation for learning the potential energy surface (PES).
The red line denotes a 2D representation of the continuous PES of materials.
The blue line illustrates the effective PES, which describes the energy of a
relaxed structure for a given unrelaxed input structure. Data augmentation
aims to improve the machine learning of this effective PES by better sampling
the configuration space. The black circle indicates the relaxed structures
contained in the data set, and the blue circles symbolize artificially
generated structures for the data augmentation.

The significant errors for unrelaxed structures are due to the limited
sampling of the complex multi-dimensional configuration space of the potential
energy surface (PES), with the relaxed structures only describing the minima
of the surface. Since unrelaxed structures are not located at these minima,
predicting a structure’s formation energy at an unrelaxed configuration is a
qualitatively different task Goodall  _et al._ (2021). To sample
configurations near the minima of the PES, Smith et al. Smith  _et al._ (2017)
applied a data augmentation technique known as normal mode sampling to a large
dataset of molecular structures, resulting in a high fidelity neural network
potential. However, determining the normal modes requires millions of phonon
calculations which, while attainable for molecular structures, is infeasible
for crystal structures.

Recently, Honrao et al. Honrao  _et al._ (2020) showed that GA data could be
used to predict relaxed formation energies of unrelaxed structures to high
accuracy. This high accuracy was achieved by augmenting the training data set,
setting the formation energy of every structure within a basin of attraction
to the minima of the respective basin of attraction, essentially modeling the
continuous PES as a step function.

We propose leveraging these findings to improve formation energy prediction of
unrelaxed structures by augmenting our data using a simple, physically
motivated perturbation technique. Figure 1 illustrates our augmentation
approach, which perturbs the atomic coordinates of a relaxed structure to
generate additional training points that describe the regions surrounding the
minima of the PES. We then map these perturbed structures to the energy of the
relaxed structure, which requires no additional ab initio calculations. We
utilize the CGCNN Xie and Grossman (2018) and CGCNN-HD Noh  _et al._ (2020) to
analyze how the augmentation affects formation energy predictions. We train
these models on the MP database Jain  _et al._ (2013) augmented by perturbed
structures. The resulting CGCNN models have similar prediction errors as the
original ones for relaxed structures. To show the improvement in formation
energy prediction, we apply the models to a test set consisting of 623
unrelaxed Nb-Sr-H hydride structures produced from a GA structure search. We
find that compared to training on only relaxed structures, training with the
augmented data set, consisting of relaxed and perturbed structures, reduced
the formation energy prediction MAE from 251 meV/atom to 86 meV/atom for CGCNN
and from 172 meV/atom to 82 meV/atom for CGCNN-HD, as compared to the models
trained only on relaxed structure.

##  II Methods

For training, we use two data sets derived from the MP database Jain  _et al._
(2013) accessed on December 10, 2021. The 1st data set, referred to as the
relaxed data set, consists of 126k relaxed structures from the MP database,
20% of this data is held out for validation. The 2nd data set, referred to as
the augmented data set, consists of the relaxed set and one perturbed
structure for every relaxed structure.

We augment the data by perturbing the coordinates, $R_{i}$ of all atoms, $i$,
in each relaxed structures using a displacement vector

$\Delta\vec{R}_{i}=(M_{i}^{x},M_{i}^{y},M_{i}^{z})$ (1)

where each Cartesian component is obtained by sampling a displacement
distribution. The displacement distribution was determined by analyzing the
displacements of atoms during relaxation in three separate GA structure
searches. The details of the GA structure searches can be found in the
supplemental materials. The distance was determined by first taking the
difference between the initial and final structure’s fractional coordinates
using the minimum image convention Deiters (2013). These differences were then
multiplied by the lattice vector matrix of the relaxed structure to obtain the
cartesian displacement vector and the euclidean norm for this vector. Figure 2
shows the resulting distribution of displacements, which was then fitted to a
Gaussian mixture model (GMM) as implemented in the open python library scikit-
learn Pedregosa  _et al._ (2011). The value of $M_{pi}$ was then selected by
randomly sampling the GMM. Information about the change in lattice vectors and
volume can be found in Figure S3.

Figure 2: Distribution of the displacement of atoms during structural
relaxation in three separate GA runs. The displacement is measured by the
change in fractional coordinates multiplied by the lattice vector matrix of
the relaxed structure.

The CGCNN and CGCNN-HD are trained on both the relaxed and augmented data. The
model’s architecture was determined by performing a grid search on the CGCNN
trained on the augmented data. Table S2 provides the range of parameters
considered in the grid search. The architecture that minimized the validation
error consists of 3-graph convolutional layers followed by 6-hidden layers
with 64 neurons each. This architecture was then used for all models.
Interestingly we found that models with more than eight hidden layers suffered
from the vanishing gradient problem Goh  _et al._ (2017) when training on the
relaxed data while training on the augmented data allowed for deeper models.
The remaining model hyperparameters are set to the values reported in Ref. Xie
and Grossman (2018) for CGCNN and Ref. Noh  _et al._ (2020) for the CGCNN-HD.

Figure 3: The convex hull of the energies of the predicted structures for the
NbH2-SrH2-H2 materials system. The distance from the convex hull measures the
thermodynamics stability of the various candidate compounds. The green shading
indicates the formation energy of the thermodynamically stable compounds and
mixtures relative to the three compounds NbH2, SrH2, and H2. We identify a
previously unknown ternary hydride, Sr2NbH9.

To provide test data for our models, we performed a GA search over the ternary
system formed by H2-Sr6NbH16-Nb6SrH16. We used the Genetic Algorithm for
Structure and Phase Prediction (GASP) python package Tipton and Hennig (2013);
Revard  _et al._ (2016) for performing the GA search. Our aim with the search
was to produce high hydrogen-containing structures that might show
superconductivity. We remove the elemental hydrogen structures and partitioned
this data into two test sets consisting of the relaxed and unrelaxed hydrides,
referred to as “H-Relaxed” and “H-Unrelaxed,” respectively. Additionally.
since the MP database contains few hydride structures, this data provides a
challenging test case.

Figure 4: Learning curves for the CGCNNN (1st row) and CGCNN-HD (2nd row)
trained on the relaxed (1st column) and augmented (2nd column) data. The faded
curves show the exact loss values, while the solid curves show the smoothed
values. The red and green curves denote the loss on the training and
validation data, respectively. The orange and blue curves display the loss for
the H-Relaxed and H-Unrelaxed test sets, respectively. Note that the H-Relaxed
and H-Unrelaxed data sets were not used in the training or validation of the
model. The $r$ values are the Pearson correlation coefficients between the
stated trends.

To relax and evaluate the energies of the candidate structures generated by
GASP, we use VASP Kresse and Hafner (1993, 1994); Kresse and Furthmüller
(1996, 1996) with the projector augmented wave method Blöchl (1994) and the
Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation for the
exchange-correlation functional Perdew  _et al._ (1996). We use a $k$-point
density of 40 per inverse Å with the Methfessel-Paxton scheme and a smearing
of 100 meV for the Brillouin zone integration, and a cutoff energy of 250 eV
for the plane-wave basis set. The GA search was terminated after 771 DFT
relaxations. We recomputed all energies using the VASP inputs generated by the
MPRelaxset class of pymatgen to ensure consistency between the training and
test sets and computed the formation energies Wang  _et al._ (2021). Figure 3
shows the ternary convex hull of the Nb-Sr-H system produced using the GA-
generated structure and the known competing phases from the MP database.
Noteworthy, our structure search found a new, previously unreported,
thermodynamically stable ternary hydride, Sr2NbH9. Preliminary analysis of
Sr2NbH9 suggest that the band gap closes at around 100 GPa.

##  III Results

Figure 4 shows the training, validation, and test (H-Relaxed/H-Unrelaxed) RMSE
for each training epoch of the respective models. For interpretability, the
trends are smoothed using an exponential moving average with a smoothing
weight of 0.95. The Pearson correlation coefficients Freedman  _et al._
(2007), with a value of 1 for perfect correlation and -1 for perfect anti-
correlation, are computed between smoothed trends of H-Relaxed/H-Unrelaxed,
H-Relaxed/validation, and H-Unrelaxed/validation.

Figure 5: The Parity plots for the CGCNNN (1st row) and CGCNN-HD (2nd row)
trained on the relaxed (1st column) and augmented (2nd column) data. The
$x$-axis denotes the DFT computed formation energies, while the $y$-axis
denotes the predicted formation energies. The values reported in the lower
right are the coefficient of determination, $R^{2}$, the MAE, and the RMSE.

The Pearson correlation coefficients show that the CGCNN and CGCNN-HD trained
on only relaxed structures (Figures 4(a) and (b)) result in anti-correlated
trends between predictions on the relaxed and unrelaxed structures of the test
set. These findings demonstrate that accurate predictions on relaxed
structures do not lead to accurate predictions on unrelaxed structures and
provide insight into the high prediction error for unrelaxed structure inputs
reported in the literature. Furthermore, the CGCNN and CGCNN-HD trained on
only relaxed structures display anti-correlated trends between the H-Unrelaxed
test set and the validation set. This anti-correlation is detrimental to the
model’s predictive performance on unrelaxed structures because the validation
error shows that predictions are improving and training should continue, while
in actuality, unrelaxed structure predictions are getting worse. Additionally,
the anti-correlation leads to an inability to correctly optimize a model’s
hyperparameters.

Figures 4(c) and (d) illustrate the effectiveness of training with the
augmented dataset. Simply perturbing the atomic coordinates of each relaxed
structure structure and then training on both the relaxed and perturbed
structures dramatically improved the models’ predictive ability on unrelaxed
structures. Both the CGCNN and CGCNN-HD trained on the augmented data set show
a high Pearson correlation between H-Unrelaxed/validation, enabling the
effective optimization of a model’s hyperparameters and the implementation of
early stopping.

Figure 5 compares the models’ formation energy predictions on the unrelaxed
hydride test data to the DFT-computed formation energies. The CGCNN in Figure
5 (a) trained on only relaxed data tends to over predict $\mathrm{E_{f}}$ for
the higher energy hydrides, which leads to a significant prediction MAE of 251
meV/atom. The added regularization applied to the CGCNN-HD in Figure 5 (b)
improves the predictions on the higher energy hydrides. Still, the model tends
to over predict $\mathrm{E_{f}}$ leading to an MAE of 172 meV/atom when
training on relaxed.

Training with the augmented dataset substantially improves the prediction MAEs
for both the CGCNN and CGCNN-HD, reducing the prediction MAE to 86 meV/atom
and 82 meV/atom, respectively. While the CGCNN-HD trained on augmented data
has the lowest testing MAE, the dense region of underpredicted formation
energies seen in Figure 5 (d) leads to a substantial number of misclassified
unstable structures, hindering the model’s ability to filter unstable
structures. This will be discussed further in the proceeding section.

Due to the large perturbations, many augmented structures have likely moved to
neighbor basins of attraction. Intuitively this would seem to yield
substantial errors. However, as shown previously, the prediction of unrelaxed
structures improved substantially. We suspect the error associated with
perturbing a structure to a neighboring basin of attraction is mitigated due
to the tendency of neighboring basins to cluster around similar minima on the
PES Revard  _et al._ (2014). Still, the large perturbations likely introduce
some prediction error. A more sophisticated augmentation method could likely
reduce the number of structures perturbed into neighboring basins and further
improve predictions.

Interestingly, the models trained on the augmented data also display improved
predictions on H-Relaxed (Figure S1). However, this improvement seems to be
specific to our test data as predictions on the relaxed validation data were
better when the model was trained only on relaxed structures(Figure S2). The
CGCNN-HD trained on only relaxed structures under predicted the structures in
H-Relaxed likely because the training data contains relatively few transition
metal hydrides and the bounded hyperbolic activation function, utilized in the
CGCNN-HD’s convolutional layers, has poor predictive power on unseen domains
Kim  _et al._ (2021). This poor predictive power on unseen domains is also the
reason the CGCNN-HD models make poor predictions on the the high and low
formation energy structures of the MP data.

###  III.1 Filtering unstable hydrides

To evaluate the models’ ability to filter energetically unfavorable
structures, we removed the MP database’s correction applied to hydrogen-
containing compounds and constructed a convex hull using the five known
competing phases of the Nb-Sr-H system. Then, based on this constructed convex
hull, we computed the hull distance $\mathrm{E_{Hull}^{DFT}}$of all the
structures in the test set, utilizing their DFT-computed formation energy, and
defined all structures with $\mathrm{E_{Hull}^{DFT}<0}$ as stable. A total of
10 structures in the test set met the stability criteria.

To construct the receiver operating characteristic (ROC) curve, shown in
Figure 6 the predicted formation energies was used to compute hull distance
($\mathrm{E_{Hull}^{ML}}$). To compute a range of true positives, false
positives, true negatives, and false negatives, we varied the stability
criteria of $\mathrm{E_{Hull}^{ML}}$ over a range of hull distances that
ensure a completed ROC curve for. We defined a true positive as a stable
structure predicted as stable, a false positive as an unstable structure
predicted as stable, a true negative as an unstable structure predicted as
unstable, and a false negative as a stable structure predicted as unstable.

Figure 6: Receiver operating characteristic curve for the CGCNN and CGCNN-HD.
the $aug$ subscript represents the model was trained on the augmented data.
the $rel$ subscript represents the model was trained on the relaxed data. The
dashed-black line represents a random classifier. The AUC is reported for all
models.

To further evaluate the models, we assume our test data is randomly generated
and consider two hypothetical cases in which the models may be utilized. Case
1 emulates a study where identifying all stable structures is desirable (TPR =
1.0). Case 2 emulates a study where missing some stable structures is
acceptable (TPR = 0.7).

For case 1, the CGCNN trained on the augmented data performed best,
successfully identifying all stable structures at a filtration criteria of
$\mathrm{E_{Hull}^{ML}<39~{}meV/atom}$ while misclassifying 130 structures,
yielding a 5-fold reduction in the number of energy calculations needed to
identify all stable structures. The CGCNN-HD trained on the augmented data
misclassified 320 structures with a filtration criteria of
$\mathrm{E_{Hull}^{ML}<85~{}meV/atom}$. The CGCNN-HD trained on relaxed data
misclassified 442 structures with optimal filtration criteria of
$\mathrm{E_{Hull}^{ML}<328~{}meV/atom}$. The CGCNN trained on relaxed data
misclassified 491 structures with optimal filtration criteria of
$\mathrm{E_{Hull}^{ML}<576~{}meV/atom}$.

For case 2, again, the CGCNN trained on the augmented data performed best,
obtaining a TPR of 0.7 and only misclassifying 54 structures with a filtration
criteria of $\mathrm{E_{Hull}^{ML}<6~{}meV/atom}$. The CGCNN-HD trained on
augmented data misclassified 89 structures with a filtration criteria of
$\mathrm{E_{Hull}^{ML}<21~{}meV/atom}$. The CGCNN trained on relaxed data
misclassified 110 structures with a filtration criteria of
$\mathrm{E_{Hull}^{ML}<159~{}meV/atom}$. The CGCNN-HD trained on relaxed data
misclassified 130 structures a filtration criteria of
$\mathrm{E_{Hull}^{ML}<91~{}meV/atom}$.

The naive approach of setting the stability criteria to be the same for
$\mathrm{E_{Hull}^{ML}}$ and $\mathrm{E_{Hull}^{DFT}}$, restricts the ability
to select a balance of accuracy and computational cost. for example at
stability criteria of $\mathrm{E_{Hull}^{ML}<0}$, the CGCNN trained with the
augmented dataset obtained a TPR of 0.6, misclassifying 47 structures. While
this performance is acceptable, as shown previously, at stability criteria of
$\mathrm{E_{Hull}^{ML}<39~{}meV/atom}$ the model can correctly identify all
stable structures which may be more desirable for a given application.

##  IV Conclusion

We proposed a simple, physically motivated, computationally cheap perturbation
technique that augmented our data to better represent the PES, dramatically
improving unrelaxed structure formation energy predictions. When compared to
training on only relaxed structures, training with an augmented data set
consisting of one relaxed and one perturbed structure for every relaxed
structure, prediction MAEs of the CGCNN and CGCNN-HD was reduced from 251
meV/atom and 172 meV/atom to 86 meV/atom and 82 meV/atom, respectively.
Further, we showed that predictions on relaxed structures showed an anti-
correlation to predictions of unrelaxed structures. We then showed that
augmenting the data enabled the predictions of relaxed structures to follow a
similar trend to that of unrelaxed structures. Finally, we utilize a ROC curve
to show two cases where our method may be useful in accelerating CSPA. While
there likely exist more advanced augmentation techniques, this work showed the
surprising effectiveness of a relatively simple method of augmentations, that
outperformed the current state of the art in formation energy prediction of
unrelaxed structures.

##  V Code Availability

Code for implementing the model on both cpus and gpus, training the models,
augmenting training data are available at
https://github.com/JasonGibsonUfl/Augmented_CGCNN. Data will be made available
upon reasonable requests.

##  VI Acknowledgments

This work was supported by the National Science Foundation under grants Nos.
PHY-1549132, the Center for Bright Beams, and the software fellowship awarded
to J.B.G. by the Molecular Sciences Software Institute funded by the National
Science Foundation (Grant No. ACI-1547580). Computational resources were
provided by the University of Florida Research Computing Center.

##  VII Author Contributions

JBG and RGH conceived the Augmentation technique and training strategy. JBG
implemented the Augmentation technique and performed the model training and
analysis. ACH performed the genetic algorithm structure search to provide
testing data. JBG, ACH, and RGH contributed to the writing of the manuscript.

## References

  * Jain  _et al._ (2013) A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. a. Persson, APL Materials 1, 011002 (2013).
  * Kirklin  _et al._ (2015) S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, npj Computational Materials 1, 15010 (2015).
  * Draxl and Scheffler (2018) C. Draxl and M. Scheffler, MRS Bulletin 43, 676–682 (2018).
  * Curtarolo  _et al._ (2012) S. Curtarolo, W. Setyawan, G. L. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Mehl, H. T. Stokes, D. O. Demchenko, and D. Morgan, Computational Materials Science 58, 218 (2012).
  * Oganov  _et al._ (2019) A. R. Oganov, C. J. Pickard, Q. Zhu, and R. J. Needs, Nature Reviews Materials 4, 331–348 (2019).
  * Revard  _et al._ (2018) B. Revard, W. Tipton, and R. Hennig, “Genetic algorithm for structure and phase prediction,” https://github.com/henniggroup/GASP-python (2018).
  * Wales and Doye (1997) D. J. Wales and J. P. K. Doye, The Journal of Physical Chemistry A 101, 5111 (1997), https://doi.org/10.1021/jp970984n .
  * Noh  _et al._ (2019) J. Noh, S. Kim, G. h. Gu, A. Shinde, L. Zhou, J. M. Gregoire, and Y. Jung, Chem. Commun. 55, 13418 (2019).
  * Chen  _et al._ (2021) B. Chen, L. J. Conway, W. Sun, X. Kuang, C. Lu, and A. Hermann, Phys. Rev. B 103, 035131 (2021).
  * Duan  _et al._ (2014) D. Duan, Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao, H. Yu, B. Liu, W. Tian, and T. Cui, Scientific Reports 4 (2014), 10.1038/srep06968.
  * Liu  _et al._ (2017) H. Liu, I. I. Naumov, R. Hoffmann, N. W. Ashcroft, and R. J. Hemley, Proceedings of the National Academy of Sciences 114, 6990 (2017).
  * Peng  _et al._ (2017) F. Peng, Y. Sun, C. J. Pickard, R. J. Needs, Q. Wu, and Y. Ma, Physical Review Letters 119 (2017), 10.1103/physrevlett.119.107001.
  * Sun  _et al._ (2019) Y. Sun, J. Lv, Y. Xie, H. Liu, and Y. Ma, Physical Review Letters 123 (2019), 10.1103/physrevlett.123.097001.
  * Cataldo  _et al._ (2021) S. D. Cataldo, C. Heil, W. von der Linden, and L. Boeri, Physical Review B 104 (2021), 10.1103/physrevb.104.l020511.
  * Hilleke and Zurek (2022) K. P. Hilleke and E. Zurek, Journal of Applied Physics 131, 070901 (2022).
  * Huang  _et al._ (2021) Y. Huang, Y. Cheng, and J. Zhang, Industrial & Engineering Chemistry Research 60, 2737 (2021).
  * Tipton  _et al._ (2014) W. W. Tipton, C. A. Matulis, and R. G. Hennig, Computational Materials Science 93, 133 (2014).
  * Nguyen  _et al._ (2013) M. C. Nguyen, J.-H. Choi, X. Zhao, C.-Z. Wang, Z. Zhang, and K.-M. Ho, Phys. Rev. Lett. 111, 165502 (2013).
  * Heiles and Johnston (2013) S. Heiles and R. L. Johnston, International Journal of Quantum Chemistry 113, 2091 (2013), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24462 .
  * Xie  _et al._ (2021) S. R. Xie, M. Rupp, and R. G. Hennig, “Ultra-fast interpretable machine-learning potentials,” (2021), arXiv:2110.00624 [cond-mat.mtrl-sci] .
  * Wu  _et al._ (2013) S. Q. Wu, M. Ji, C. Z. Wang, M. C. Nguyen, X. Zhao, K. Umemoto, R. M. Wentzcovitch, and K. M. Ho, Journal of Physics Condensed Matter 26 (2013), 10.1088/0953-8984/26/3/035402.
  * Jennings  _et al._ (2019) P. C. Jennings, S. Lysgaard, J. S. Hummelshøj, T. Vegge, and T. Bligaard, npj Computational Materials 5, 1 (2019).
  * Xie and Grossman (2018) T. Xie and J. C. Grossman, Physical Review Letters 120 (2018), 10.1103/PhysRevLett.120.145301.
  * Ahmad  _et al._ (2018) Z. Ahmad, T. Xie, C. Maheshwari, J. C. Grossman, and V. Viswanathan, ACS Central Science 4, 996 (2018), pMID: 30159396, https://doi.org/10.1021/acscentsci.8b00229 .
  * Choudhary and DeCost (2021) K. Choudhary and B. DeCost, npj Computational Materials 7 (2021), 10.1038/s41524-021-00650-1.
  * Cheng  _et al._ (2021) J. Cheng, C. Zhang, and L. Dong, Communications Materials 2 (2021), 10.1038/s43246-021-00194-3.
  * Chen  _et al._ (2019) C. Chen, W. Ye, Y. Zuo, C. Zheng, and S. P. Ong, Chemistry of Materials 31, 3564 (2019), https://doi.org/10.1021/acs.chemmater.9b01294 .
  * Park and Wolverton (2020) C. W. Park and C. Wolverton, Phys. Rev. Materials 4, 063801 (2020).
  * Noh  _et al._ (2020) J. Noh, G. H. Gu, S. Kim, and Y. Jung, Journal of Chemical Information and Modeling 60, 1996 (2020), pMID: 32208718, https://doi.org/10.1021/acs.jcim.0c00003 .
  * Srivastava  _et al._ (2014) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, Journal of Machine Learning Research 15, 1929 (2014).
  * Goodall  _et al._ (2021) R. E. A. Goodall, A. S. Parackal, F. A. Faber, R. Armiento, and A. A. Lee, “Rapid discovery of novel materials by coordinate-free coarse graining,” (2021), arXiv:2106.11132 [cond-mat.mtrl-sci] .
  * Smith  _et al._ (2017) J. S. Smith, O. Isayev, and A. E. Roitberg, Chemical Science 8, 3192 (2017).
  * Honrao  _et al._ (2020) S. J. Honrao, S. R. Xie, and R. G. Hennig, Journal of Applied Physics 128, 085101 (2020).
  * Deiters (2013) U. K. Deiters, Zeitschrift für Physikalische Chemie 227, 345 (2013).
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, J. Mach. Learn. Res. 12, 2825–2830 (2011).
  * Goh  _et al._ (2017) G. B. Goh, N. O. Hodas, and A. Vishnu, Journal of Computational Chemistry 38, 1291 (2017), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.24764 .
  * Tipton and Hennig (2013) W. W. Tipton and R. G. Hennig, Journal of Physics: Condensed Matter 25, 495401 (2013).
  * Revard  _et al._ (2016) B. C. Revard, W. W. Tipton, A. Yesypenko, and R. G. Hennig, Physical Review B 93, 054117 (2016).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner (1994) G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Computational Materials Science 6, 15 (1996).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Wang  _et al._ (2021) A. Wang, R. Kingsbury, M. McDermott, M. Horton, A. Jain, S. P. Ong, S. Dwaraknath, and K. A. Persson, Scientific Reports 11, 15496 (2021).
  * Freedman  _et al._ (2007) D. Freedman, R. Pisani, and R. Purves, Pisani, R. Purves, 4th edn. WW Norton & Company, New York (2007).
  * Revard  _et al._ (2014) B. C. Revard, W. W. Tipton, and R. G. Hennig, “Structure and stability prediction of compounds with evolutionary algorithms,” in  _Prediction and Calculation of Crystal Structures: Methods and Applications_, edited by S. Atahan-Evrenk and A. Aspuru-Guzik (Springer International Publishing, Cham, 2014) pp. 181–222.
  * Kim  _et al._ (2021) Y. Kim, Y. Kim, C. Yang, K. Park, G. X. Gu, and S. Ryu, npj Computational Materials 7, 140 (2021).


# Production of heavy neutron-rich nuclei in transfer reactions within the
dinuclear system model

Long Zhu1,2,111Email<EMAIL_ADDRESS>Zhao-Qing Feng3,222Email:
<EMAIL_ADDRESS>and Feng-Shou Zhang1,2,4,333Corresponding author:
<EMAIL_ADDRESS>1The Key Laboratory of Beam Technology and Material
Modification of Ministry of Education, College of Nuclear Science and
Technology, Beijing Normal University, Beijing 100875, China
2Beijing Radiation Center, Beijing 100875, China
3Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000,
China
4Center of Theoretical Nuclear Physics, National Laboratory of Heavy Ion
Accelerator of Lanzhou, Lanzhou 730000, China

###### Abstract

The dynamics of nucleon transfer processes in heavy-ion collisions is
investigated within the dinuclear system model. The production cross sections
of nuclei in the reactions 136Xe+208Pb and 238U+248Cm are calculated and the
calculations are in good agreement with the experimental data. The transfer
cross sections for the 58Ni+208Pb reaction are calculated and compared with
the experimental data. We predict the production cross sections of neutron-
rich nuclei 165-168Eu, 169-173Tb, 173-178Ho, and 181-185Yb based on the
reaction 176Yb+238U. It can be seen that the production cross sections of the
neutron-rich nuclei 165Eu, 169Tb, 173Ho, and 181Yb are 2.84 $\mu$b, 6.90
$\mu$b, 46.24 $\mu$b, and 53.61 $\mu$b, respectively, which could be
synthesized in experiment.

##  1 Introduction

The production of superheavy nuclei (SHN) have been studied extensively in the
past three decades. Much progress has been made experimentally [1, 2, 3, 4]
and theoretically [5, 6, 7, 8, 9, 10, 11]. There are mainly two sorts of
reaction mechanism to produce SHN, which are fusion reactions and
multinulcleon transfer reactions. The synthesis of superheavy nuclei is
motivated with respect to searching the “island of stability”. The heavy
nuclei produced in fusion reactions are not so neutron-rich. Therefore, it is
quite hard to reach the center of the “island of stability” through the fusion
process. The damped collisions between heavy nuclei could be one possible
approach.

The light and medium mass nuclei far from the stability line (neutron- and
proton-rich) can be produced in fragmentation process or fission of heavy
nuclei. The heavy neutron-rich nuclei including those located at the “island
of stability” are not synthesized and studied yet, which could be synthesized
by using multinucleon transfer reactions.

The fundamental mechanisms of the deep inelastic collisions have been
investigated many years ago [12, 13]. In recent years, the production of
neutron-rich heavy and superheavy nuclei have been discussed by Zagrebaev and
Greiner [14, 15]. The multinucleon transfer processes are also investigated by
the time-dependent Hartree-Fock (TDHF) approach [16, 17] and improved quantum
molecular dynamics (ImQMD) model [18]. The dinuclear system (DNS) model has
been successfully used in investigating mechanisms of the synthesis of SHN in
fusion reactions [19, 20, 22, 21, 23] and the transfer processes of the damped
collisions [24, 25, 26].

In this work, the transfer cross sections are calculated within the framework
of DNS model and compared with the available experimental data. The cross
sections of transferring nucleons from light to heavy fragments and from heavy
to light fragments are both investigated. We predict the production cross
sections of the unknown neutron-rich nuclei through the transfer reaction
176Yb+238U.

The article is organized as follows. In Sec. 2, we describe the theoretical
model. The results and discussion are presented in Sec. 3. We summarize the
main results in Sec. 4.

##  2 Model description

The diffusion process is treated along proton and neutron degrees of freedom.
The distribution probability is obtained by solving a set of master equations
in the potential energy surface (PES) of the DNS. The time evolution of the
mass asymmetry is described by the following master equation [23]:

$\displaystyle\frac{dP(Z_{1},N_{1},E_{1},t)}{dt}=\sum_{Z_{1}^{{}^{\prime}}}W_{Z_{1},N_{1};Z_{1}^{{}^{\prime}},N_{1}}(t)[d_{Z_{1},N_{1}}P(Z_{1}^{{}^{\prime}},N_{1},E_{1}^{{}^{\prime}},t)-d_{Z_{1}^{{}^{\prime}},N_{1}}P(Z_{1},N_{1},E_{1},t)]$
$\displaystyle+\sum_{N_{1}^{{}^{\prime}}}W_{Z_{1},N_{1};Z_{1},N_{1}^{{}^{\prime}}}(t)[d_{Z_{1},N_{1}}P(Z_{1},N_{1}^{{}^{\prime}},E_{1}^{{}^{\prime}},t)-d_{Z_{1},N_{1}^{{}^{\prime}}}P(Z_{1},N_{1},E_{1},t)]$
$\displaystyle-\Lambda_{qf}(\Theta(t))P(Z_{1},N_{1},E_{1},t).$ (1)

Here $P(Z_{1},N_{1},E_{1},t)$ is the probability distribution function of the
fragment with proton number $Z_{1}$ and neutron number $N_{1}$ with
corresponding excitation energy $E_{1}$ at time $t$.
$W_{Z_{1},N_{1};Z_{1}^{{}^{\prime}},N_{1}}=W_{Z_{1},N_{1};Z_{1},N_{1}^{{}^{\prime}}}$
is the mean transition probability from the channel ($Z_{1}$, $N_{1}$,
$E_{1}$) to ($Z_{1}^{{}^{\prime}}$, $N_{1}$, $E_{1}^{{}^{\prime}}$) [or
($Z_{1}$, $N_{1}$, $E_{1}$) to ($Z_{1}$, $N_{1}^{{}^{\prime}}$,
$E_{1}^{{}^{\prime}}$)]. $d_{Z_{1},N_{1}}$ denotes the microscopic dimension
corresponding to the macroscopic state ($Z_{1}$, $N_{1}$, $E_{1}$). The
details of $W$ and $d_{Z_{1},N_{1}}$ are given in Ref. [27]. The sum is taken
over all possible proton and neutron numbers that fragment
$Z_{1}^{{}^{\prime}}$ and $N_{1}^{{}^{\prime}}$ may take, but only one nucleon
transfer is considered in the model. $\Lambda_{qf}$ describes the quasifission
rate.

The evolution of the DNS along the relative distance $R$ leads to QF of the
DNS. The QF rate $\Lambda_{qf}(\Theta(t))$ can be treated with the one-
dimensional Kramers rate [24],

$\displaystyle\Lambda_{qf}(\Theta(t))=$
$\displaystyle\frac{\omega}{2\pi\omega^{B_{qf}}}[\sqrt{(\frac{\Gamma}{2\hbar})^{2}+(\omega^{B_{qf}})^{2}}-\frac{\Gamma}{2\hbar}]\times
exp[-\frac{B_{qf}(Z_{1},N_{1})}{\Theta(t)}].$ (2)

The QF rate exponentially depends on the QF barrier $B_{qf}(Z_{1},N_{1})$. The
local temperature $\Theta(t)$ is calculated by using Fermi-gas expression
$\Theta=\sqrt{E_{1}/a}$ with the local excitation energy $E_{1}$ and the
level-density parameter $a=A/12$ $\textrm{MeV}^{-1}$. The frequency
$\omega^{B_{qf}}$ of the inverted harmonic oscillator approximates the
potential $V$ in $R$ at the top of the quasifission barrier, and $\omega$ is
the frequency of the harmonic oscillator approximating the potential in $R$
around the bottom of the pocket. The $\Gamma$ determines the friction
coefficients. Here, $\Gamma=2.8$ MeV, $\hbar\omega^{B_{qf}}=2.0$ MeV, and
$\hbar\omega=3.0$ MeV.

The local excitation energy is defined as

$\displaystyle E_{1}=E_{diss}-$
$\displaystyle(U(Z_{1},N_{1},Z_{2},N_{2},J)-U(Z_{P},N_{P},Z_{T},N_{T},J))$ (3)
$\displaystyle-\frac{(J-M)^{2}}{2\zeta_{rel}}-\frac{M^{2}}{2\zeta_{int}}.$

Here, $U(Z_{1},N_{1},Z_{2},N_{2})$ and $U(Z_{P},N_{P},Z_{T},N_{T})$ are the
driving potentials of fragment $A_{1}$ and the entrance point of the DNS.
$E_{diss}$ is the excitation energy of the composite system, which is
converted from the relative kinetic energy loss. $M$ denotes the intrinsic
angular momentum derived from the dissipation of the relative angular
momentum, and $\zeta_{int}$ is the corresponding moment of inertia. $J$
denotes the initial angular momentum. $\zeta_{rel}$ is the relative moment of
inertia of the DNS, which is given by $\zeta_{rel}=\mu R^{2}_{m}$. $R_{m}$ is
the distance of two colliding nuclei located at bottom of potential pocket.
For heavy systems with no potential pocket the transfer processes take place
at the touching configuration.

During the diffusion process, the relative kinetic energy will dissipate into
the DNS system. The local excitation energy is determined by the excitation
energy of DNS and the PES [23]. The PES of DNS is given as

$\displaystyle
U(Z_{1},N_{1},Z_{2},N_{2},R,\beta_{1},\beta_{2},\theta_{1},\theta_{2},J)=U_{LD}(Z_{1},N_{1})+U_{LD}(Z_{2},N_{2})-$
$\displaystyle
U_{LD}(Z,N)-V^{CN}_{rot}(J)+V_{CN}(Z_{1},N_{1},Z_{2},N_{2},R,\beta_{1},\beta_{2},\theta_{1},\theta_{2},J).$
(4)

$U_{LD}(Z_{1},N_{1})$, $U_{LD}(Z_{2},N_{2})$ and $U_{LD}(Z,N)$ are the binding
energies of the fragments $A_{i}$ and compound nucleus $A$, respectively. The
value of
$U(Z_{1},N_{1},Z_{2},N_{2},R,\beta_{1},\beta_{2},\theta_{1},\theta_{2},J)$ is
normalized to the energy of the rotating compound nucleus by
$U_{LD}(Z,N)+V^{CN}_{rot}(J)$. $V_{CN}$ is the interaction potential of two
fragments, which depends on the deformation parameter and the orientation of
the deformed fragments ($\theta_{1}$ and $\theta_{2}$). The details of
$V_{CN}$ are given in Ref. [24]. The distance R between the center of two
fragments is chosen to be the value at $R_{m}$ or the touching configuration,
in which the DNS is assumed to be formed.

The sharing of the excitation energy between primary fragments is assumed to
be proportional to their masses. The code GEMINI is used to treat the
sequential statistical evaporation of excited fragments. Subsequent de-
excitation cascades of the excited fragments via emission of light particles
(neutron, proton, and $\alpha$) and gamma-rays competing with fission process
are taken into account, which lead to the final mass distribution of the
reaction products.

Figure 1:  Driving potentials of the tip-tip (solid line) and side-side
(dashed line) collisions in the reaction 136Xe+208Pb.

The production cross section of a primary fragment with charge $Z_{1}$ and
mass number $A_{1}$ can be calculated as follows:

$\displaystyle\sigma^{pr}_{Z_{1},A_{1}}(E_{c.m.})=\sum^{J_{max}}_{J=0}\sigma^{pr}_{Z_{1},A_{1}}(E_{c.m.},J)=\sum^{J_{max}}_{J=0}\sigma_{cap}(E_{c.m.},J)P(Z_{1},A_{1},J,E_{1},\tau_{int}),$
(5)

where $\sigma_{cap}(E_{c.m.},J)$ is the partial capture cross section, which
defines the probability of overcoming the Coulomb barrier. $E_{c.m.}$ is the
incident energy in center of mass system. The details of calculation of
$\sigma_{cap}(E_{c.m.},J)$ is given in Ref. [23]. $\tau_{int}$ is the
interaction time of the diffusion process at the bottom of the potential
pocket, which is dependent on the incident energy $E_{c.m.}$ and the incident
angular momentum $J$. $\tau_{int}$ cannot be measured directly. In this work,
$\tau_{int}$ is determined by using the deflection function method [28].

##  3 Results and discussion

Figure 2:  The time evolution of the primary mass distributions of the
reaction 136Xe+208Pb at $E_{c.m.}=514$ MeV with $t_{0}=10^{-22}$ s. $t_{df}$
is the interaction time calculated by using deflection function method [28].

Figure 1 shows the driving potentials as a function of the mass number of
fragments for two cases of tip-tip and side-side orientations in the reaction
136Xe+208Pb. The arrows show the positions of the masses of projectile and
target in entrance channel. A valley can be seen around the nucleus 208Pb for
both orientations, which is because 208Pb is a double magic nucleus. It is
also can be seen that the valley is deeper for side-side orientation. In this
work, we investigate the dissipation process under the potential energy
surface with tip-tip orientation.

Figure 3:  The comparison of calculated production cross sections for primary
(thin lines) and survived (thick lines) isotopes of Po (solid lines), Rn
(dashed lines), and Ra (dash-dotted lines) with the experimental data [29] in
the reaction 136Xe+208Pb at $E_{c.m.}=514$ MeV. The solid points correspond to
the experimental cross sections for 210Po (square), 222Rn (triangle), and
224Ra (star).

The time evolution of the primary mass distributions, for the reaction
136Xe+208Pb at $E_{c.m.}=514$ MeV, is shown in Fig. 2. It can be seen that the
distribution is around the entrance channel with $\tau_{int}=5\times 10^{-22}$
s. With the increasing interaction time the distributions become wider.
Equilibrium is reached at about $\tau_{int}=40\times 10^{-22}$ s. The
distribution with $\tau_{int}=t_{df}$ is also shown. $t_{df}$ is calculated by
using deflection function method. It is indicated that the system also can
reach equilibrium with $\tau_{int}=t_{df}$, which can be seen from mass
asymmetry region. The mass distribution around entrance channel for
$\tau_{int}=t_{df}$ is higher than that for $\tau_{int}=40\times 10^{-22}$ s.
This is because $t_{df}$ depends on the angular momentum $J$. The interaction
time decreases with the increasing angular momentum.

Figure 4:  Transfer cross sections for the 58Ni+208Pb reaction at
$E_{c.m.}=$256.8 MeV. The number of transferred protons from 58Ni to 208Pb
(negative number) is indicated. The measured cross sections have been reported
in Ref. [16].

In Fig. 3, the calculated production cross sections of isotopes of Po, Rn, and
Ra in transfer reaction 136Xe+208Pb are shown. The bombarding energy
$E_{c.m.}=$514 MeV. The calculated results are the case of tip-tip collisions,
which have the height of the interaction potential at the touching
configuration with the value 421 MeV. The larger qusifission rate of such
systems with no potential pocket results in the DNS quickly decaying into two
fragments. Therefore, the neutron transfer process is governed by the driving
potential with short interaction time. Within the error bars, the calculated
results are in good agreement with the experimental data [29] for production
of 210Po and 224Ra, while the cross section of 222Rn are underestimated.

Figure 5:  Transfer cross sections for the 58Ni+208Pb reaction at
$E_{c.m.}=$256.8 MeV. The number of transferred neutrons from 208Pb to 58Ni
(positive number) is indicated. The measured cross sections have been reported
in Ref. [16].

Figure 4 shows transfer cross sections of the reaction 58Ni+208Pb at
$E_{c.m.}=$256.8 MeV. Each panel of Fig. 4 shows cross sections according to
the change of proton number of projectilelike fragment (PLF) from 58Ni, as
functions of neutron number of PLF. The calculated results show a good
agreement with the experimental data at neutron rich region, while the
calculated results underestimate the experimental data at neutron lack region.
As the transferred proton number decreases, the peak position of curve shifts
towards larger neutron number, while the peaks of the experimental data
locates at a neutron number of 30. The reason probably is that we only
consider the driving potential with tip-tip orientation.

Figure 6:  Cross sections for the formation of isotopes of Einsteinium
($Z=99$) (solid lines), Fermium ($Z=100$) (dashed lines), and Mendelevium
($Z=101$) (dash-dotted lines) in the reaction 238U+248Cm at $E_{c.m.}=800$
MeV. The thin and thick lines are distribution of primary and final fragments,
respectively. The experimental data are taken from Ref. [30]. The circles
denote the unknown nuclei.

Figure 5 shows transfer cross sections of the reaction 58Ni+208Pb at
$E_{c.m.}=$256.8 MeV according to the change of neutron number of PLF from
58Ni, as functions of proton number of PLF. From each panel of Fig. 5 one can
see that the calculated results are in good agreement with the experimental
data except that the proton number is 28 in first panel. Therefore, the
process of transferring some neutrons from 208Pb to 58Ni and transferring some
protons from 58Ni to 208Pb can be better described.

Figure 7:  (a) Production cross sections of isotopes of Yb in the transfer
reaction 176Yb+238U at $E_{c.m.}=570$ (dash-dotted line), 600 (solid line),
650 (dashed line), and 690 MeV (dotted line). (b) Cross sections for the
formation of isotopes of the elements Ytterbium (solid line), Holmium (dashed
line), Terbium (dotted line), and Europium (dash-dotted line) in the reaction
176Yb+238U at $E_{c.m.}=600$ MeV. The circles denote the unknown neutron-rich
nuclei.

In order to further test the model, the cross sections for the production of
heavy actinides in damped collisions 238U+248Cm at $E_{c.m.}=800$ MeV are
shown in Fig. 6. The interaction time of this system is very short due to the
strong Coulomb repulsion. One can see a good agreement between theoretical and
experimental cross sections [30] for production of isotopes of the elements
Es, Fm, and Md. The cross sections decrease drastically with increasing atomic
numbers of the fragments. The nuclei 261Md and 262Md could be synthesized in
this reaction with cross sections of about 1.52 nb and 0.17 nb, respectively.
Reaction fragments formed in transfer reactions of heavy ions are strongly
excited. It can be seen that the survival probabilities of most of primary
fragments are quite low due to a dominant fission channel and the cross
sections of primary fragments are a few hundreds times larger than those of
the final fragments. The incident energy should also satisfy the condition
that two colliding heavy nuclei can come into contact and have enough
interaction time for nucleon transfer. Therefore, the optimal incident energy
should be found for the largest yield of heavy neutron-rich nuclei.

In Fig. 7 the calculated yields of heavy neutron-rich nuclei in the transfer
reaction 176Yb+238U are shown. Figure 7 (a) shows the distributions of
isotopes of the element Yb formed in the reaction 176Yb+238U at different
incident energies. The height of the interaction potential at the touching
configuration is about 565 MeV. The height of curve increases with the
increasing incident energy from $E_{c.m.}=570$ to 650 MeV, while it decreases
from $E_{c.m.}=650$ to 690 MeV. This is because with the increasing incident
energy, the excitation energy of these fragments increases and thus decreases
their survival probabilities. It can be seen that cross sections are
relatively larger in neutron-rich region when $E_{c.m.}=600$ MeV.

The calculated cross sections for production of Europium ($Z=63$), Terbium
($Z=65$), Holmium ($Z=67$), and Ytterbium ($Z=70$) isotopes in the damped
collision 176Yb+238U at $E_{c.m.}=600$ MeV are shown in Fig. 7 (b). The
circles denote the unknown heavy neutron-rich nuclei. For the production of
neutron-rich nuclei, the cross sections decrease dramatically with the
increasing neutron number. We show in Fig. 7 (b) the neutron-rich nuclei with
production cross section larger than 1 nb. The calculated production cross
sections for the neutron-rich nuclei 165,166,167,168Eu are 2.84 $\mu$b, 0.78
$\mu$b, 19.64 nb, and 1.36 nb, 169,170,171,172,173Tb are 6.90 $\mu$b, 3.70
$\mu$b, 0.12 $\mu$b, 44.44 nb, and 1.00 nb, 173,174,175,176,177,178Ho are
46.24 $\mu$b, 16.83 $\mu$b, 2.00 $\mu$b, 0.85 $\mu$b, 73.75 nb, and 18.08 nb,
and 181,182,183,184,185Yb are 53.61 $\mu$b, 4.67 $\mu$b, 2.85 $\mu$b, 84.70
nb, and 4.29 nb, respectively. One can see that the cross sections for
synthesis of the heavy neutron-rich nuclei 165Eu, 169Tb, 173Ho, and 181Yb are
2.84 $\mu$b, 6.90 $\mu$b, 46.24 $\mu$b, and 53.61 $\mu$b, respectively, which
are quite large for experimental detection.

##  4 Conclusions

The calculated production cross sections of heavy nuclei in transfer reaction
136Xe+208Pb are in good agreement with the available experimental data. The
transfer cross sections of the reaction 58Ni+208Pb, according to the change of
neutron (proton) number of PLF from 58Ni, as functions of proton (neutron)
number of PLF are investigated. For the production of heavy neutron-rich
isotopes of actinides, the reaction 238U+248Cm is also studied. The calculated
transfer cross sections show good agreement with the experimental data. It is
reasonable to use the DNS model to study the dynamics of transfer reactions.
By using the transfer reaction 176Yb+238U, the neutron-rich nuclei 165Eu,
169Tb, 173Ho, and 181Yb could be synthesized with the cross sections of 2.84
$\mu$b, 6.90 $\mu$b, 46.24 $\mu$b, and 53.61 $\mu$b, respectively.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China
under Grand Nos. 11025524 and 11161130520, National Basic Research Program of
China under Grant No. 2010CB832903, the European Commission’s 7th Framework
Programme (FP7-PEOPLE-2010-IRSES) under Grant Agreement Project No. 269131,
and the National Natural Science Foundation of China Projects (Nos 11175218
and U1332207).

## References

## References

  * [1] Hofmann S and Münzenberg G 2000 Rev. Mod. Phys. 72 733\.
  * [2] Oganessian Yu Ts 2007 J. Phys. G 34 R165.
  * [3] Morita K, Morimoto K, Kaji D, Haba H, Ozeki K, Kudou Y, Sumita T, Wakabayashi Y, Yoneda A, Tanaka K, Yamaki S, Sakai R, Akiyama T, Goto S, Hasebe H, Huang M, Huang T, Ideguchi E, Kasamatsu Y, Katori K, Kariya Y, Kikunaga H, Koura H, Kudo H, Mashiko A, Mayama K, Mitsuoka S, Moriya T, Murakami M, Murayama H, Namai S, Ozawa A, Sato N, Sueki K, Takeyama M, Tokanai F, Yamaguchi T, and Yoshida A 2012 J. Phys. Soc. Jpn. 81 103201\.
  * [4] Khuyagbaatar J, Yakushev A, Düllmann Ch E, Ackermann D, Andersson L -L, Asai M, Block M, Boll R A, Brand H, Cox D M, Dasgupta M, Derkx X, Nitto A Di, Eberhardt K, Even J, Evers M, Fahlander C, Forsberg U, Gates J M, Gharibyan N, Golubev P, Gregorich K E,Hamilton J H, Hartmann W, Herzberg R -D, Heßberger F P, Hinde D J, Hoffman-n J, Hollinger R, Hübner A, Jäger E, Kindler B, Kratz J V, Krier J, Kurz N, Laatiaoui M, Lahiri S, Lang R, Lommel B, Maiti M, Miernik K, Mina-mi S, Mistry A, Mokry C, Nitsche H, Omtvedt J P, Pang G K, Papadakis P, Renisch D, Roberto J, Rudolph D, Runke J, Rykaczewski K P, Sarmien-to L G, Schädel M, Schausten B, Semchenkov A, Shaughnessy D A, Steinegger P, Steiner J, Tereshatov E E, Thörle-Pospiech P, Tinschert K, De Heiden-reich T Torres, Trautmann N, Türler A, Uusitalo J, Ward D E, Wegrzecki M, Wiehl N, Van Cleve S M, and Yaku-sheva V 2014 Phys. Rev. Lett. 112 172501\.
  * [5] Siwek-Wilczyńska K, Cap T, Kowal M, Sobiczewski A, and Wilczyński J 2012 Phys. Rev. C 86 014611\.
  * [6] Wang N, Zhao E G, Scheid W, and Zhou S G 2012 Phys. Rev. C 85 041601(R).
  * [7] Shen C W, Abe Y, Boilley D, Kosenko G, and Zhao E G 2008 Int. J. Mod. Phys. E 17 66\.
  * [8] Abe Y, Marchix A, Shen C, Yilmaz B, Kosenko G, Boilley D, and Giraud B G 2007 Int. J. Mod. Phys. E 16 491\.
  * [9] Mandaglio G, Giardina G, Nasirov A K, and Sobiczewski A 2012 Phys. Rev. C 86 064607\.
  * [10] Zhu L, Xie W J, and Zhang F S 2014 Phys. Rev. C 89 024615\.
  * [11] Wong C Y 1976 Phys. Rev. Lett. 37 664\.
  * [12] Cassing W and Nörenberg W 1983 Nucl. Phys. A 401 467\.
  * [13] Nörenberg W 1974 Phys. Letts. B 53 289\.
  * [14] Zagrebaev V and Greiner W 2005 J. Phys. G 31 825\.
  * [15] Zagrebaev V and Greiner W 2007 J. Phys. G 34 2265\.
  * [16] Sekizawa K and Yabana K 2013 Phys. Rev. C 88 014614\.
  * [17] Golabek C and Simenel C 2009 Phys. Rev. Lett. 103 042701\.
  * [18] Tian J, Wu X, Zhao K, Zhang Y, and Li Z 2008 Phys. Rev. C 77 064603\.
  * [19] Adamian G G, Antonenko N V, and Scheid W 2004 Phys. Rev. C 69 011601(R).
  * [20] Adamian G G, Antonenko N V, Scheid W, and Volkov V V 1998 Nucl. Phys. A 633 409\.
  * [21] Adamian G G, Antonenko N V, and Scheid W 2000 Nucl. Phys. A 678 24\.
  * [22] Adamian G G, Antonenko N V, and Scheid W 2004 Phys. Rev. C 69 014607\.
  * [23] Zhu L, Feng Z Q, Li C, and Zhang F S 2014 Phys. Rev. C 90 014612\.
  * [24] Adamian G G, Antonenko N V, and Scheid W 2003 Phys. Rev. C 68 034601\.
  * [25] Feng Z Q, Jin G M, and Li J Q 2009 Phys. Rev. C 80 067601\.
  * [26] Kalandarov Sh A, Adamian G G, Antonenko N V, and Scheid W 2010 Phys. Rev. C 82 044603\.
  * [27] Feng Z Q, Jin G M, Fu F and Li J Q 2006 Nucl. Phys. A 771 50\.
  * [28] Li J Q and Wolschin G 1983 Phys. Rev. C 27 590\.
  * [29] Kozulin E M, Vardaci E, Knyazheva G N, Bogachev A A, Dmitriev S N, Itkis I M, Itkis M G, Knyazev A G, Loktev T A, Novikov K V, Razinkov E A, Rudakov O V, Smirnov S V, Trzaska W, and Zagrebaev V I 2012 Phys. Rev. C 86 044611\.
  * [30] Schädel M, Brüchle W, Gäggeler H, Kratz J V, Sümmerer K, Wirth G, Herrmann G, Stakemann R, Tittel G, Trautmann N, Nitschke J M, Hulet E K, Lougheed R W, Hahn R L, and Ferguson R L 1982 Phys. Rev. Lett. 48 852\.