text:
# Fusion process of Lennard-Jones clusters: global minima and magic numbers formation.

## I Introduction

It is well known that the sequence of cluster magic numbers carries essential information about the electronic and ionic structure of the cluster LesHouches . Understanding of the cluster magic numbers is often equivalent or nearly equivalent to the understanding of cluster electronic and ionic structure. A good example of this kind is the observation of the magic numbers in the mass spectrum of sodium clusters Knight84 . In this case, the magic numbers were explained by the delocalized electron shell closings (see MetCl99  and references therein). Another example is the discovery of fullerenes, and, in particular, the $`C_{60}`$ molecule Kroto85 , which was made by means of the carbon clusters mass spectroscopy.

The formation of a sequence of cluster magic numbers should be closely connected to the mechanisms of cluster fusion. It is natural to expect that one can explain the magic number sequence and find the most stable cluster isomers by modeling mechanisms of cluster assembly and growth. On the other hand, these mechanisms are of interest on their own, and the correct sequence of the magic numbers found in such a simulation can be considered as a proof of the validity of the cluster formation model.

The problem of magic numbers is closely connected to the problem of searching for global minima on the cluster multidimensional potential energy surface. The number of local minima on the potential energy surface increases exponentially with the growth cluster size and is estimated to be of the order of $`10^{43}`$ for $`N=100`$ LesHouches . Thus, searching for global minima becomes increasingly difficult problem for large clusters. There are different algorithms and methods of the global minimization, which have been employed for the global minimization of atomic cluster systems.

One of the most widely used global optimization methods calls the simulated annealing Kirkpatrick ; Wille ; Ma93 ; Ma94 ; Tsoo ; Laaarhoven . This method is an extension of metropolis Monte Carlo techniques Kalos . In the simulated annealing one starts calculation from a high energy state of the system and then cools the system down by decreasing its kinetic energy. In standard application, the system is annealed using molecular dynamics or Monte-Carlo based methods, but also more sophisticated variants of this algorithm are used, such as gaussian density annealing and analogues Ma93 ; Ma94 ; Tsoo .

A related method is quantum tunnelling Finnila ; Amara ; Sylvain . This method attempts to reduce the effects of barriers on the potential energy surface. This is done by allowing the system to behave quantum mechanically, leading to the possibility of tunnelling. However, the computational expense of this method grows exponentially with increasing the number of atoms in the system, what makes this method applicable for relatively small clusters Maranas92 ; Maranas94 .

Another class of global optimization methods is based on smoothing the potential energy surface Stillinger ; Piela ; Pillardy . Within the framework of this method, a transformation is applied to the potential energy surface in order to decrease the number of high lying local minima. The global minimum of the smoothed potential energy surface, which is found by steepest descent energy minimization method, or by Monte-Carlo search, is then mapped back to the original surface.

Another successful method of global optimization is the basin-hopping algorithm LesHouches ; Wales ; Leary . This method involves a potential energy transformation, which does not change neither the global minimum, nor the relative energies of local minima as it is usually done in the potential energy surface smoothing methods. In this method, the potential energy surface in any point of the configuration space is assigned to that of the local minimum obtained by the given geometry optimization method and the potential energy surface is mapped onto a collection of interpenetrating staircases with basins corresponding to the set of configurations which lead to a given minimum after optimization.

Successful global optimization results for LJ clusters reported to date were obtained with the use of genetic algorithms Goldberg ; Niesse ; Gregurick ; Deaven ; Romero . These methods are based on the idea that the population of clusters evolves to low energy by mutation and mating of structures with low potential energy. There are different versions of genetic algorithm (see papers cited above). One of them is seeding algorithm, or seed growth method Niesse ; Gregurick . In this method it is assumed that the most stable cluster of $`N`$ particles can be obtained from the most stable cluster of $`N1`$ particle by a random addition of one particle near the boundary of the cluster and then by optimizing the structure with a given optimization method. Genetic algorithms and other methods related to them are based on stochastic minimization technique, and require no quantum calculation or much information about the potential energy surface, as the previously mentioned algorithms.

The new algorithm GrowProcPRL  which we describe in detail in this work for the first time is based on the dynamic searching for the most stable cluster isomers in the cluster fusion process. We call this algorithm as the cluster fusion algorithm (CFA). Our calculations demonstrate that the CFA can be considered as an efficient alternative to the known techniques of the cluster global minimization. The big advantage of the CFA consists in the fact that it allows to study not just the optimized cluster geometries, but also their formation mechanisms.

In the present work we approach the formulated problem in a simple, but general, form. In our simplest scenario, we assume that atoms in a cluster are bound by LJ potential and the cluster fusion takes place atom by atom. At each step of the fusion process all atoms in the system are allowed to move, while the energy of the system is decreased. The motion of the atoms is stopped when the energy minimum is reached. The geometries and energies of all cluster isomers found in this way are stored and analyzed. On the first glance, this method is somewhat similar to the genetic algorithms and to the seed growth method, in particular. However, the principle difference between the CFA and the genetic algorithms consists in the fact that in CFA the new atoms are fused to the system not randomly, but to certain specific points in the system, like the cluster centre of mass, the centres of mass of the cluster faces or some specific points in the vicinity of the cluster surface (see section II for detail). Such an approach, allows one to model various physical situations and scenarios of the cluster fusion process, which are beyond the scope of the stochastic seed-growth method. Thus, with the use of the CFA, we have determined the fusion paths for all known global energy minimum LJ cluster configurations in the size range up to $`N150`$. We have also determined the fusion paths for the chains of clusters based on the icosahedral, octahedral, decahedral and tetrahedral symmetries. We report the striking correspondence of the peaks in the dependence on cluster size of the second derivative of the binding energy per atom calculated for the chain of LJ clusters based on the icosahedral symmetry with the peaks in the abundance mass spectra experimentally measured for the clusters of noble gas atoms Echt81 ; Haberland94 ; Harris84 .

Our paper is organized as follows. In section II, we describe theoretical model for the cluster fusion process and the CFA. In section III, we present and discuss the results of computer simulations based on the CFA. In section IV, we draw a conclusion to this paper. In Appendix A, we provide additional details of the numerical algorithms developed in our work.

## II Theoretical model for cluster fusion process

A group of atoms bound together by interatomic forces is called an atomic cluster (AC). There is no qualitative distinction between small clusters and molecules. However, as the number of atoms in the system increases, ACs acquire more and more specific properties making them unique physical objects different from both single molecules and from the solid state.

Each stable cluster configuration corresponds to a minimum on the multidimensional potential energy surface of the system. The number of local minima on the potential energy surface increases very rapidly with the growth cluster size. Thus, searching for the most stable cluster configurations possessing the absolute energy minimum, the so-called global energy minimum structures, becomes increasingly difficult problem for large clusters.

Searching for such global energy minimum cluster structures is one of the focuses of our paper. This problem by itself is not new as it is clear from introduction. However, in our work we approach this problem at a new perspective and investigate not just the global energy minimum cluster structure by its own, but rather the formation of such structures in the cluster fusion process.

From the physical point of view, it is natural to expect that one can find the most stable cluster isomers by modeling the mechanisms of cluster assembly and growth, i.e. the cluster fusion process GrowProcPRL . This idea is based on the fact that in nature, or in laboratory, the most stable clusters are often obtained namely in the fusion process. This idea is general and is applicable to different types of clusters.

In this work we consider ACs formed by atoms interacting with each other via the pairing force. The interaction potential between two atoms in the cluster can, in principle, be arbitrary. For concrete computations, we use the Lennard-Jones (LJ) potential,

$$U(r)=4\epsilon \left(\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right),$$
(1)
where $`r`$ is the interatomic distance, $`\epsilon `$ is the depth of the potential well ($`\epsilon >0`$), $`2^{1/6}\sigma `$ is the pair bonding length. The constants in the potential allow one to model various types of clusters for which LJ paring force approximation is reasonable. The most natural systems of this kind are the clusters consisting of noble gas atoms Ne, Ar, Kr, Xe formed by van der Waals forces. The constants of the LJ potential appropriate for the noble gas atoms one can find in RS . Thus, for Ne, Ar, Kr, Xe, $`\epsilon =3.6`$, $`12.3`$, $`17.2`$, $`24.3`$ meV and $`\sigma =3.1`$, $`3.8`$, $`4.0`$, $`4.4`$ $`\AA `$ respectively. Note that for the LJ clusters it is always possible to choose the coordinate and energy scales so that $`\sigma =1`$ and $`\epsilon =1`$. It makes all LJ cluster systems scalable. They only differ by the choice of $`\sigma `$, $`\epsilon `$ and the mass of a single constituent (atom). In our paper we use LJ potential with $`\sigma =1`$ and $`\epsilon =1/12`$.

The AC treatment with the use of the LJ forces often implies the applicability of the classical Newton equations for the description of the AC dynamics. Following this line, we describe the atomic motion in the cluster by the Newton equations with the LJ pairing forces. In this case, the information about quantum properties of the system is hidden in the LJ-potential constants $`\sigma `$ and $`\epsilon `$. In computations, the system of coupled Newton equations for all atoms in the cluster is solved numerically using the 4th order Runge-Kutta method.

Experimentally, mass spectra of the Ne, Ar, Kr, Xe clusters have been investigated in Echt81 ; Haberland94 ; Harris84 . In figure 1, we compile the results of these measurements. This figure demonstrates that the mass spectra for the clusters of noble gas atoms have many common features. Some differences in the spectra can be attributed to the differences in the experimental conditions in which the spectra have been measured.

The peaks in the spectra indicate the enhanced stability of certain clusters. These clusters are called the magic clusters and their mass numbers $`N`$ are the magic numbers. The enhanced stability of the magic clusters can have different origin for various types of clusters. For LJ clusters, the origin of magic numbers is connected with the formation and filling the icosahedral shells of atoms. The completed icosahedral shells of atoms correspond to the following sequence of magic numbers:

$$N=\frac{10}{3}z^3+5z^2+\frac{11}{3}z+1$$
(2)
as it is clear from a simple geometry analysis. Here, the integer $`z=1,2,3,4\mathrm{}`$ is the order of the icosahedral shell. The first four icosahedral magic numbers, N, as they follow from (2) are equal to 13, 55, 147, 309. In figure 1, one can see, however, many more peaks corresponding to the formation of magic clusters having the partially completed icosahedral shells. In this paper, we elucidate the origin of all the peaks in the recorded mass spectra of the noble gas ACs in the size range up to $`N150`$ by modeling the mechanisms of cluster fusion and by finding in this way the most stable cluster isomers, i.e. the global energy minimum cluster structures.

To solve this problem, for each $`N`$, we need to find solutions of the Newton equations leading to the stable cluster configurations and then to choose the one, which is energetically the most favourable. The choice of the initial conditions for the simulation and the algorithm for the solution of this problem are described below.

Each stable cluster configuration corresponds to a minimum on the cluster potential energy surface, i.e. to the point in which all the atoms in the system are located in their equilibrium positions. A minimum can be found by allowing atoms to move, starting from a certain initial cluster configuration, and by absorbing all their kinetic energy in the most efficient way. If the starting cluster configuration for $`N+1`$ atoms has been chosen on the basis of the global minimum structure for $`N`$ atoms, then it is natural to assume, and we prove this in the present work, that very often the global minimum structure for $`N+1`$ atoms can be easily found. The success of this procedure reflects the fact that in nature clusters in their global minima often emerge, namely, in the cluster fusion process, which we simulate in such calculation.

We have employed the following algorithm for the kinetic energy absorption. At each step of the calculation we consider the motion of one atom only, which undergoes the action of the maximum force. At the point, in which the kinetic energy of the selected atom is maximum, we set the absolute value of its velocity to zero. This point corresponds to the minimum of the potential well at which the selected atom moves. When the selected atom is brought to the equilibrium position, the next atom is selected to move and the procedure of the kinetic energy absorption repeats. The calculation stops when all the atoms are in equilibrium.

We have considered a number of scenarios of the cluster fusion on the basis of the developed algorithm for finding the stable cluster configurations.

In the simplest scenario clusters of $`N+1`$ atoms are generated from the $`N`$-atomic clusters by adding one atom to the system. In this case the initial conditions for the simulation of $`(N+1)`$-atomic clusters are obtained on the basis of the chosen $`N`$-atomic cluster configuration by calculating the coordinates of an extra atom fussed to the system on a certain rule. We have probed the following paths: the new atom can be fussed either

(A1) to the center of mass of the cluster, or

(A2) randomly outside the cluster, but near its surface, or

(A3) to the centres of mass of all faces of the cluster, or

(A4) to the points that are close to the centres of all faces of the cluster, located from both sides of the face on the perpendicular to it, or

(A5) to the centres of mass of the faces laying on the cluster surface.

Here, the cluster surface is considered as a polyhedron, so that the vertices of the polyhedron are the atoms and two vertices are connected with an edge, if the distance between them is less than the value $`d=1.2d_0`$, where $`d_0=2^{1/6}\sigma `$ is the bonding length of a free pair of atoms. The cluster surface is considered as a polyhedron that covers all the atoms in the system and has the minimum volume. The whole cluster volume is divided on a sum of attached triangle, square and pentagonal pyramids. In the A3 method, all the faces in the cluster are counted as a sum of the faces laying on the cluster surface and the faces of the pyramids filling the cluster volume.

The choice of the method how to fuse atoms to the system depends on the problem to be solved. The A1 method is appropriate in situations when the new atoms are fused into the cluster volume. The A2 simulates the process when atoms are fused to the cluster surface. By both A1 and A2 methods, large clusters consisting of many particles can be generated rather quickly. The A2 method is especially fast, because fusion of one atom to the boundary of the cluster usually does not lead to the recalculation of its central part. The A3 and A4 methods can be used for searching the most stable, i.e. energetically favourable, cluster configurations or for finding cluster isomers with some other specific properties. The A4 method leads to finding more cluster isomers than the A3 one, but it takes more CPU time. The A5 method is especially convenient for modeling the cluster fusion process which we focus on in this paper. Using this method one can generate the cluster growth paths for the most stable cluster isomers.

In the cluster fusion process, new atoms should be added to the system starting from the initially chosen cluster configuration step by step until the desired cluster size is reached. Each new step in the cluster fusion process should be made with the use of the methods A1-A5. By methods A3-A5, one generates many different cluster isomers at each new step of the fusion process, those number grows rapidly as N increases. It is not feasible and often it is not necessary to fuse new atoms to the all found isomers. Instead, one can fuse atoms only to the selected clusters, which are of interest. Below, we outline a number of selection criteria, which can be used for the cluster selection in the fusion process depending on the problem to be solved. At each new step of the fusion process:

(SE1) one of the clusters with the minimum number of atoms is selected, or

(SE2) the cluster with the minimum energy among the already found stable clusters of the maximum size is selected, or

(SE3) a few low energy cluster isomers among the already found stable clusters of the maximum size are selected, or

(SE4) the cluster with the maximum energy among the already found stable clusters of the maximum size is selected, or

(SE5) the cluster possessing a significant structural change among the already found stable clusters is selected.

The SE1 criterion is relevant in the situation, when the full search of cluster isomers is needed. It is applicable to the systems with relatively small number of particles. Note that each cluster can be selected only once. The SE2 criterion is the most relevant for modelling the cluster fusion process and the purposes of the present work. It turns out to be very efficient and often leads to finding the most stable cluster configurations for a given number of particles. However, this is not always the case, particularly for large clusters. In these situations, the better results on the cluster global optimization can be achieved with the use of SE3 method. However, this method takes more CPU time. The SE4 and SE5 criteria turn out to be useful for the redirection of the cluster fusion process towards the lower energy cluster isomer branches.

Calculations performed with the use of the methods described above show that often clusters of higher symmetry group possess relatively low energy. Thus, the symmetric cluster configurations are often of particular interest. The process of searching for the symmetric cluster configurations can be speed up significantly, by the impose of the symmetry constraints in the cluster fusion process. This means that for obtaining a symmetric $`N`$ atomic cluster isomer from the initially chosen symmetric $`(NM)`$-atomic configuration one should fuse $`M`$ atoms to the surface of this isomer symmetrically.

This goal can be achieved by various methods.

(SY1) The planes of symmetry of the parent cluster are to be found. Their intersections determine the axes of symmetry. The atoms of the parent cluster do not move and held their position, then an atom is added outside the cluster surface. The initial coordinates of it are not important, because its stable position in the vicinity of the frozen cluster surface is calculated. Then, this atom is reflected many times with respect to the found symmetry planes, and the atoms are added on the places of its images. The obtained configuration is symmetric also. After that, all the atoms in the system are allowed to move while their kinetic energy is absorbed. If the obtained stable cluster configuration is symmetric, it can be used as the initial cluster for further computations of this type. Using this procedure several times the chains of clusters of certain symmetry can be found.

(SY2) This method is similar to the one described above, but now the atoms are added to the centres of the faces on the cluster surface. If the cluster has two (or more) different types of faces, then, at first, the atoms are added to the centres of the faces of the first type and their positions are optimized with frozen cluster core. Then, the whole system is optimized and checked for symmetry. If it is symmetric, the process continues, atoms are added to the faces of the second type and the same actions are performed. This method makes possible to find more stable configurations than the first one.

(SY3) The atoms are added to the axes of symmetry outside the cluster surface. Then, the atoms of the parent cluster are assumed to be frozen and the optimized coordinates of the added atoms are calculated. After that the optimization of the whole system is performed. If the obtained cluster configuration is symmetric, it can be used as the initial cluster for further computations.

Note that the SY1-SY3 methods do not model any concrete physical scenario for the cluster fusion process. They rather represent efficient mathematical algorithms for the generation symmetrical cluster configurations, which often turn out to be very useful for the cluster structure analysis.

Note also that in addition to the cluster fusion techniques described above any stable cluster configuration can be manually modified by adding or removing a number of atoms. These modifications can also be performed using advantages of the cluster symmetry if it does exist. The manual modification of the cluster configuration is useful when the cluster geometry is already established or nearly established and only small modifications of the system are required. For example, by this method one can perform the surface rearrangement of atoms leading to the growth the number of bonds in the cluster and thus to the enhancement of its stability. This goal can be achieved by the replacement of a few loosely bound atoms on the cluster surface to the positions, in which they have larger number of bonds. As it is demonstrated in the next section, such rearrangements result in finding the cluster isomers, which alternatively can be obtained when performing the fusion simulation with the use of the SE3 criterion. The manual modifications technique requires usually much less CPU time. Cluster configurations found in this way can be used as starting configurations for the subsequent cluster fusion process.

## III Results and discussions

Using our algorithms, we have examined various paths of the cluster fusion process, and determined the most stable cluster isomers for the cluster sizes of up to 150 atoms.

We have generated the chains of clusters based on the icosahedral, octahedral, tetrahedral and decahedral symmetries with the use of the A3-A5 and SE2-SE5 methods. We show that clusters possessing icosahedral and decahedral type of lattice are energetically more favourable. For the calculated cluster chains, we analyze the average interatomic distances in clusters, average number of bonds and binding energies. Using these results we explain the sequence of magic numbers, experimentally measured for the noble gas clusters, and elucidate the level of applicability of the liquid drop model for the description of the LJ cluster systems. Using the A1-A2 methods we calculate chains of energetically unfavourable cluster isomers as an example of the spontaneous cluster growth of such systems. Also, we consider the formation of symmetrical clusters using methods SY1-SY3 and their correspondence to the global energy minimum cluster structures.

### III.1 Fusion of global energy minimum clusters

#### III.1.1 LJ cluster geometries

The growth of the most stable, i.e. possessing the lowest energy, LJ clusters based on the icosahedral symmetry, the so-called global energy minimum cluster structures, is illustrated in figure 2, parts a-d. We analyse the cluster geometries within the size range $`4N150`$ and determine the transition path from smaller clusters to larger ones. The cluster geometries have been determined using the methods described in section II. Namely, we have used A3-A5 methods combined with SE2-SE5 cluster selection criteria.

For each cluster, we show the cluster axis, i.e. a group of atoms located on a line, that goes through the cluster centre of mass. The cluster axis can be linear (see the clusters with higher symmetry, for example, $`LJ_{13}`$, $`LJ_{55}`$, $`LJ_{71}`$, $`LJ_{147}`$) or slightly distorted (see, for example, non-completed cluster configurations $`LJ_{25}`$, $`LJ_{27}`$, $`LJ_{76}`$, $`LJ_{116}`$). Other atoms in the cluster are located around the cluster axis and form the completed or open polygonal rings (see figure 2). The bonds between the atoms laying in the axis and in the rings are shown by thick and thin lines respectively. In some particular cases, we present additional connections in order to stress specific structural properties of some clusters (see, for example, $`LJ_{11}`$, $`LJ_{31}`$, $`LJ_{51}`$, $`LJ_{82}`$).

Figure 2 demonstrates how the cluster of a certain size can be obtained from the smaller one. For each cluster, we draw the added atoms by grey circles. For some clusters, it is necessary additionally to replace one atom, or even a few atoms from one position to another in order to reach the global energy minimum. In these cases, we mark the positions, from which the atoms are removed, by grey rings (see, for example, $`LJ_{96}`$, $`LJ_{102}`$, $`LJ_{126}`$).

Figure 2 demonstrates that the majority of energetically favourable cluster structures can be obtained with the use of the A5 method combined with the SE2 cluster selection criterion, according to which the global energy minimum cluster geometry is obtained from the preceding cluster configuration by fusion a single atom to the cluster surface. Such situation takes place for $`N=`$ 4-17, 19-26, 28-29, 32-34, 36-37, 39-50, 52-64, 67, 70-72, 74-75, 77, 79-81, 83, 89-92, 94-95, 97, 99-101, 103-104, 106-110, 112, 114, 116-120, 122, 124-125, 127, 129, 131-132, 134, 136-150, i.e. in about $`75\%`$ cases.

The global energy minimum cluster structures presented in figure 2 coincide with those found by other cluster global optimization techniques (see LesHouches ; Wales ; Leary  and references therein). Note that the different choice of the constants $`\epsilon `$ and $`\sigma `$ in various calculation schemes does not change the cluster geometry and influences only the scale of the system’s size and its total energy. However, the existing global optimization techniques do not allow to perform the cluster growth analysis, which we carry out in our work for the first time. Some preliminary results of our work have been published in GrowProcPRL ; NanoMeeting .

It is interesting that all the calculated cluster geometries have the structure, in which a number of completed and open polygons round the cluster axis. The maximum number of atoms in polygons depends on the cluster size. Within the size range of $`N150`$, as it is seen from figure 2, the pentagonal, decagonal and pentadecagonal polygons are possible. The fact that this type of rings within the cluster structure dominates is closely connected with the prevalence of the icosahedral type of packing of atoms in the LJ clusters.

Note that only the clusters possessing relatively high symmetry have all the polygons rounding the axis completed. For the most of the clusters, the polygons laying close to the cluster surface are open (see e.g. $`LJ_{61}`$, $`LJ_{89}`$, $`LJ_{124}`$). For such clusters, the atoms added to the system occupy the positions in the open polygons located near the cluster surface.

#### III.1.2 Surface atoms rearrangements

Our simulations demonstrate that the fusion of a single atom to the global energy minimum cluster structure of $`N`$ atoms, in some cases does not immediately lead to a global energy minimum of $`(N+1)`$-atomic cluster. This happens for $`N=`$ 18, 27, 30, 35, 38, 51, 65, 66, 68, 69, 73, 76, 78, 84, 86, 87, 88, 93, 96, 98, 102, 105, 111, 113, 115, 121, 123, 126, 128, 130, 133, 135. The global energy minimum cluster structure for these cluster sizes can not be found directly with the use of A5 method combined with the SE2 criterion from the global energy minimum cluster configuration of the smaller size. In these cases the SE2 criterion leads to finding the higher energy cluster isomers, which we also call the excited state cluster isomers. In order to determine the global energy minimum cluster configurations for the above mentioned $`N`$, it is necessary to rearrange additionally one or a few atoms in the excited state cluster isomer (see example shown in figure 3). In figure 2, we illustrate this procedure by marking the positions of the rearranged atoms by grey circles. The initial positions of the atoms are marked by grey rings. This figure shows that the rearrangement takes place always in the vicinity of the cluster surface. This fact has a simple physical explanation. The surface atoms are bound weaker than those inside the cluster volume, and thus they are movable easier and allow the favourable surface atomic rearrangement. This consideration demonstrates that the surface rearrangement of atoms is an essential component of the cluster growth process.

We illustrate the surface rearrangement of atoms in figure 3. The first row in figure 3 shows that the fusion of a single atom to the global energy minimum cluster structure of 26 atoms does not lead to the global energy minimum of the $`LJ_{27}`$ cluster. Thus, the rearrangement of surface atoms in the $`LJ_{27}`$ cluster leading to the formation of the global energy minimum cluster structure is needed. The necessary rearrangement of atoms is shown in the second row of figure 3. The result of such rearrangement can be obtained if one starts the cluster growth from the excited state cluster isomer of the $`LJ_{25}`$ cluster, as it is illustrated in the third row of figure 3.

In this particular example, only the two smaller cluster sizes are involved in the fusion process of the $`LJ_{27}`$ cluster. However, the cluster fusion via excited states is not always that simple and evident. In some cases, it involves more than 10 intermediate steps. Such situation takes place, for example, for the $`LJ_{69}`$ cluster, which can be obtained from the excited state isomer of the $`LJ_{55}`$ cluster.

#### III.1.3 Average number of bonds in LJ clusters

The rearrangement of the surface atoms can be performed algorithmically by freezing almost all the atoms and allowing only a few atoms on the cluster surface to move. However, the number of possible cluster configurations increases rapidly with the growth cluster size, and it becomes not feasible to optimize and investigate all of them. Thus, one has to use a simple criterion for the selection of some clusters, for which the surface rearrangement of atoms is of interest. For such an analysis, we have introduced the criterion based on counting the number of bonds in the cluster. We consider the two atoms as bound if the distance between them is smaller than the given cut-off distance value $`d`$. Defining bonding between the atoms in such a way, we then calculate the total number of bonds in the cluster. This characteristic can be used in searching for the required cluster rearrangement, because the cluster with the maximum number of bonds, possesses the highest binding energy.

In figure 4, we present the dependence of average number of bonds in the cluster as a function of the cluster size. For this calculation we have used the cut-off value $`d=1.2d_0`$, where $`d_0=2^{1/6}\sigma `$ is the LJ potential bonding length. In figure 4, circles and squares show the average number of bonds for the global energy minimum structures presented in figure 2 and for the excited state cluster isomers with the energy close to the global minimum correspondingly. This figure demonstrates that circles always show the upper limit for the average number of bonds in the cluster.

For some clusters, like $`LJ_{27}`$ or $`LJ_{30}`$, the number of bonds in the global energy minimum cluster structure and the closest excited state cluster isomer appears to be very close. Such situations occur, when the energies of the two clusters turn out to be very close. In such cases, bonding between atoms at distances larger than the chosen cut-off value becomes important and needs to be counted in order to see that the number of bonds in the global energy minimum cluster structure is maximal. For example, the total number of bonds for the two different cluster isomers of $`LJ_{30}`$ calculated with the cut-off value $`d`$ is equal to $`121`$. Doubling the value of $`d`$ results in the total number of bonds $`343`$ and $`339`$ for the global energy minimum and the excited state cluster isomer respectively.

Such an analysis turns out to be very useful for the algorithmic search of the surface rearrangement of atoms leading to the growth of the number of bonds in the cluster and its global optimization. In simple cases, the favourable rearrangement is often obvious and it can be performed manually, by the replacement of a group of atoms from one position to another. This is often the case for larger clusters, see, for example, the rearrangement in the $`LJ_{38}`$ or $`LJ_{69}`$ clusters.

Note that the average number of bonds in the clusters with icosahedral symmetry converges to the bulk limit relatively slowly. As it is clear from a simple geometrical analysis, for an infinitively large icosahedron the average number of bonds should be equal to $`6`$, while for $`LJ_{55}`$, $`LJ_{147}`$ and $`LJ_{309}`$ it is 4.25, 4.73 and 5.00 respectively.

#### III.1.4 LJ cluster lattices

The methods and criteria described above appear to be insufficient for the complete description of the LJ cluster fusion process. At certain cluster sizes, a radical rearrangement of the cluster structure is necessary. Below, we call such rearrangements as the core or lattice rearrangements. Indeed, clusters in the fusion chain form the lattice of a certain type. Different lattices are based on different principles of the cluster packing and symmetry. For LJ clusters, we analyzed four different lattices: icosahedral, decahedral, octahedral and tetrahedral. The cluster configurations found in our simulations have been distinguished between the four groups according to their core symmetry. The cluster core is formed by a certain number of completed shells of atoms. Clusters consisting of the completed shells only are called the principal magic clusters. The mass numbers for the principal magic clusters (principal magic numbers) with the icosahedral, decahedral, octahedral and tetrahedral type of packing can be easily determined from simple geometric consideration and read as:

$`N_i`$ $`=`$ $`{\displaystyle \frac{10}{3}}z^3+5z^2+{\displaystyle \frac{11}{3}}z+1`$ (3)
$`N_d`$ $`=`$ $`{\displaystyle \frac{5}{6}}z^3+{\displaystyle \frac{5}{2}}z^2+{\displaystyle \frac{8}{3}}z+1`$
$`N_o`$ $`=`$ $`{\displaystyle \frac{2}{3}}z^3+2z^2+{\displaystyle \frac{7}{3}}z+1`$
$`N_t`$ $`=`$ $`{\displaystyle \frac{1}{6}}z^3+1z^2+{\displaystyle \frac{11}{6}}z+1`$
where $`N_i`$, $`N_d`$, $`N_o`$ and $`N_t`$ are the principal magic numbers for the icosahedral, decahedral, octahedral and tetrahedral cluster packing, $`z`$ is the number of shells in the cluster or its order.

In figure 5, we show the geometrical structure of a few principal magic clusters of each type. All the clusters are formed by polygons rounding the cluster axis. Their number and type varies for different clusters.

The decahedral principal magic clusters are formed by pentagonal rings, that are located one above another along the cluster axis (see row 1 in figure 5). The number of rings and their size depends on the number of shells, $`z`$. Note that in decahedral clusters, the lines connecting corresponding vertices of pentagonal rings of the same size are parallel to the cluster axis.

The icosahedral principal magic clusters differ from the decahedral ones (see row 2 in figure 5). Depending on the order of cluster, the icosahedral cluster can additionally to pentagonal rings contain decagonal, pentadecagonal and etc ones. These additional rings glue together two equal decahedrons of the same order having the common axis, along which the decahedrons are turned one with respect to another on the angle $`36^o`$. As a result, some of neighbouring rings in the icosahedral clusters are also turned one with respect to another (see figure 5).

In row 3 and 4 of figure 5, we show the geometries of the principal octahedral and tetrahedral magic clusters respectively. In clusters of this type, squares and triangles correspondingly round the cluster axis.

#### III.1.5 LJ cluster lattice rearrangements

Let us now consider the role played by the lattice rearrangement in the formation of the global energy minimum cluster structures. Figure 2 demonstrates that the global energy minimum for $`LJ_{31}`$, $`LJ_{82}`$ and $`LJ_{85}`$ can not be found by the surface rearrangement of atoms.

Figure 2 clearly demonstrates that the $`LJ_{31}`$ cluster has obvious elements of the decahedral packing. In order to stress this fact, we connect the corresponding vertices of the pentagons by thin lines. In the clusters of the smaller size, the neighbouring pentagonal rings (completed or uncompleted) are turned one with respect to another. This feature is characteristic for the icosahedra type of cluster packing. The $`LJ_{31}`$ cluster structure arises in the cluster growing process via a long chain of excited state cluster isomers with $`N13`$. This chain is shown in figure 6. This figure demonstrates that the essential rearrangement occurs already in the $`LJ_{18}`$ cluster leading to the formation of the two co-ordinated pentagonal rings located one above another. After that point, the formation of the $`LJ_{31}`$ cluster takes place in a regular way and can be generated atom by atom.

The next prominent rearrangement of the cluster core takes place for the $`LJ_{82}`$ cluster. The formation of this cluster involves a long chain of excited state cluster isomers with $`N55`$. In figure 7, we present this chain of clusters. The core structure for the $`LJ_{82}`$ cluster emerges in the $`LJ_{59}`$ cluster, in which the three atoms from one pentagonal ring are located just above the three atoms from another pentagonal ring (see figure 7). The further formation of the $`LJ_{82}`$ cluster happens in a regular manner atom by atom.

The next cluster that can not be found from the previous one by the surface rearrangement of atoms is the $`LJ_{85}`$ cluster. However, this cluster configuration can be easily generated starting from $`LJ_{81}`$ (see figure 2, part b) and using a simple rearrangement of atoms located at the cluster surface.

To illustrate that the lattice transformations affect significantly various cluster characteristics, in figures 4 and 8, we present dependences of the average number of bonds $`n`$ and the average bonding length $`R`$ on cluster size. These figures demonstrate that for the cluster sizes $`N=31`$, $`82`$, $`85`$ the dependences have the step-like irregularities. These irregularities are caused by the cluster lattice rearrangements.

In figure 2, we present the geometries of the global energy minimum cluster structures with the icosahedral type of atomic packing. This type of packing is the most favourable for the LJ clusters. However, other type of packing is also possible and in a few cases it becomes even more favourable than the icosahedral one. Within the cluster size range considered, there are only 8 clusters for which packing other than icosahedral results in the cluster energy lower than obtained for the icosahedral clusters. These cluster structures having the octahedral, decahedral and tetrahedral type of atomic packing are presented in figure 9. They can not be obtained from their icosahedral neighbours directly. In order to find these cluster structures, one needs either to consider a chain of excited state cluster isomers analogous to the described above for the formation of the $`LJ_{31}`$ and $`LJ_{82}`$ clusters or to treat separately the growth of the octahedral decahedral and tetrahedral cluster families. Each of the cluster families should have its own magic numbers. Only a few of all these cluster configurations can compete in energy with the clusters based on the icosahedral packing as we demonstrate it in the next section. Within the size range considered such situation arises only for $`N=`$38, 75-77, 98, 102-104.

Note that in figures 4 and 8, we also show the average number of bonds $`n`$ and the average bonding length $`R`$ for non-icosahedral global energy minimum cluster isomers and compare these characteristics with $`n`$ and $`R`$ for clusters with icosahedral symmetry of the same size.

### III.2 Cluster binding energies

The binding energy per atom for LJ clusters is defined as follows:

$$E_b/N=E_N/N$$
(4)
where $`E_N`$ is the total binding energy of $`N`$-atomic cluster. In tables 1 and 2, we compile total energies and point symmetry groups for the global energy minimum cluster structures within the size range considered.

Figure 10 shows the dependence of the binding energy per atom for LJ clusters as a function of cluster size. We have generated the chains of clusters based on the icosahedral, decahedral, octahedral and tetrahedral types of lattices with the use of the A3-A5 methods combined with SE2-SE5 cluster selection criteria.

Figure 10 shows that the most stable clusters are mostly based on the icosahedral type of packing with exceptions for $`N=38`$, $`75N77`$, $`98`$ and $`102N104`$. In these cases the octahedral, decahedral, tetrahedral and again decahedral cluster symmetries respectively become more favourable. To illustrate this, in the insets to figure 10, we show the behaviour of the curves in the regions $`30N45`$ and $`70N110`$ in greater detail. In these regions, energies of the clusters with different type of packing become especially close.

In previous section, we have demonstrated that the structural cluster properties experience dramatic change at $`N=31`$, $`82`$, $`85`$. In the contrary to the structural characteristics, the dependence of binding energy per atom on $`N`$ behaves smoothly in the vicinity of these points. In order to stress the role of the cluster lattice rearrangement in the formation of the global energy minimum cluster structures, in the inset to figure 10, we plot by open squares the energies for the chains of clusters growing up from $`LJ_{30}`$ and $`LJ_{81}`$ without rearrangement of their lattice. Clusters in these chains become more and more energetically unfavourable with the growth cluster size as compared to the global energy minimum structures.

On the top of figure 10 we present the sequence of magic numbers experimentally measured for noble gases Echt81 ; Haberland94 . It is seen that the magic numbers correspond to the peculiarities in the curve for the binding energy per atom as a function of cluster size calculated for the cluster chain based on the icosahedral type of packing. The peculiarities indicate the enhanced stability of the corresponding clusters. Note that in experiment the magic numbers are seen much better in the regions of $`N`$, in which the icosahedrally packed clusters are energetically the most favourable. The most prominent peculiarities arise for the completed icosahedral shells with 13, 55 and 147 atoms. The peculiarities in the dependence of binding energy per atom on $`N`$ diminish with the growth of the cluster size due to approaching to the bulk limit, in which the binding energy per atom is constant.

### III.3 Liquid drop model

The main trend of the energy curves plotted in figure 10 can be understood on the basis of the liquid drop model, according to which the cluster energy is the sum of the volume and the surface energy contributions:

$$E_N=\lambda _VN+\lambda _SN^{2/3}\lambda _RN^{1/3}$$
(5)
Here the first and the second terms describe the volume, and the surface cluster energy correspondingly. The third term is the cluster energy arising due to the curvature of the cluster surface. Choosing constants in (5) as $`\lambda _V=0.71554`$, $`\lambda _S=1.28107`$ and $`\lambda _R=0.5823`$, one can fit the global energy minimum curve plotted in figure 10 with the accuracy less than $`1\%`$. The deviations of the energy curves calculated for various chains of cluster isomers from the liquid drop model (5) are plotted in figure 11. We have fited the energies of the icosahedral clusters in the cluster size range $`4N150`$ using the nonlinear least squares fitting method.

The peaks on these dependences indicate the increased stability of the corresponding magic clusters. The ratio between the volume and surface energies in (5) can be characterised by the dimensionless parameter $`\delta =\lambda _V/\lambda _S`$, being equal in our case to $`\delta =0.559`$. Note, that this value slightly differs from the previously reported in GrowProcPRL , $`\delta _{PRL}=0.555`$, which was obtained by fitting the cluster global energy minima within the same cluster size range.

This result can be compared with those published in Doye . In this paper the energy of the first four completed icosahedral shells was fitted using equation (5). This fitting gives the value $`\delta =0.576`$ differing from our result due to the difference in the fitting method. In Doye  the energies of the 4 selected clusters have been used, while we fit the energies of all the clusters within the size range $`4N150`$

Note that a similar model is used in nuclear physics for the description of the nuclei binding energy. For the nuclear matter, the constant $`\delta `$ is equal to $`\delta =0.903`$ BM .

### III.4 Cluster magic numbers

The dependence of the binding energies per atom for the most stable cluster configurations on $`N`$ allows one to generate the sequence of the cluster magic numbers. In figure 12, we plot the second derivatives $`\mathrm{\Delta }^2E_N`$ for the cluster chains with icosahedral, decahedral, octahedral and tetrahedral type of packing.

We compare the obtained dependences with the experimentally measured abundance mass spectra for the noble gas clusters Echt81 ; Haberland94  (see figure 1) and establish the striking correspondence of the peaks in the experimentally measured mass spectra with those in the $`\mathrm{\Delta }^2E_N`$ dependence calculated for the icosahedral type of clusters. The sequence of the magic numbers for the Ar, Xe and Kr clusters reads: 13, 19, 23, 26, 29, 32, 34, 43, 46, 49, 55, 61, 64, 71, 74, 81, 87, 91, 101, 109, 116, 119, 124, 131, 136, 147 Echt81 ; Haberland94 . The most prominent peaks in this sequence 13, 55 and 147 correspond to the closed icosahedral shells, while other numbers correspond to the filling of various parts of the icosahedral shell.

The connection between the second derivatives $`\mathrm{\Delta }^2E_N`$ and the peaks in the abundance mass spectrum of clusters one can understand using the following simple model. Let us assume that the mass spectrum of clusters is formed in the evaporation process. This means that changing the number of clusters $`n_N`$ of the size $`N`$ in the cluster ensemble takes place due to the evaporation of an atom by the clusters of the size $`N`$ and $`N+1`$, i.e.

$$\mathrm{\Delta }n_Nn_{N+1}W_{N+1N}n_NW_{NN1}$$
(6)
where the evaporation probabilities are proportional to

$`W_{N+1N}`$ $``$ $`e^{\frac{E_N+E_1E_{N+1}}{kT}}`$ (7)
$`W_{NN1}`$ $``$ $`e^{\frac{E_{N1}+E_1E_N}{kT}}`$
Here $`T`$ is the cluster temperature, $`k`$ is the Bolzmann constant. In our model $`E_1=0`$, so after simple transformations the equation for $`\mathrm{\Delta }n_N`$ reads as:

$`\mathrm{\Delta }n_N`$ $``$ $`n_{N+1}e^{\frac{E_NE_{N+1}}{kT}}(1{\displaystyle \frac{n_N}{n_{N+1}}}e^{\frac{E_{N1}2E_N+E_{N+1}}{kT}})`$
$``$ $`n_Ne^{\frac{E_NE_{N+1}}{kT}}(1e^{\frac{E_{N1}2E_N+E_{N+1}}{kT}})`$
Here we assumed that $`n_{N+1}n_Nn_{N1}1`$. Let us now estimate the relative abundances in the mass spectrum for Ar clusters for temperatures about 100 K. The exponent $`e^{\frac{E_NE_{N+1}}{kT}}`$ influences the absolute value of $`\mathrm{\Delta }n_N`$. This factor becomes exponentially small at $`kTE_NE_{N+1}`$, which for the Ar clusters means $`T800K`$, because $`\mathrm{\Delta }E_N^{Ar}=0.071eV800K`$. The small value of this factor results in the growth of the characteristic period of the evolution of $`n_N`$ with time. The factor in the brackets determines the relative cluster abundances. Indeed, its positive value for certain $`N`$ leads to the growth of the corresponding clusters in the system, while the negative value of the factor to the opposite behaviour. The factor in the brackets is characterised by $`\mathrm{\Delta }^2E_N=E_{N1}2E_N+E_{N+1}`$, which is for the Ar clusters $`|\mathrm{\Delta }^2E_N|=0.008eV90K`$. Thus, for temperatures $`T90K`$ the exponent in the brackets can be expanded. In this case one derives

$$\mathrm{\Delta }n_Nn_Ne^{\frac{E_NE_{N+1}}{kT}}\frac{\mathrm{\Delta }^2E_N}{kT}$$
(9)
This equation demonstrates, that positive values of $`\mathrm{\Delta }^2E_N`$ leads to the enhanced abundance of the corresponding clusters.

For Xe and Kr clusters this approximation becomes applicable at somewhat larger $`T`$, because the depth of the LJ potential, $`\epsilon `$, for this type of clusters, is larger than for Ar clusters. For Ar clusters $`\epsilon =12.3`$ meV while for Kr and Xe it is $`17.2`$ and $`24.3`$ meV respectively.

Within the size range considered there are only five magic numbers, which can not be directly explained on the basis of $`\mathrm{\Delta }^2E_N`$ calculated for the chain of icosahedral clusters. This situation takes place for $`N=`$34, 81, 87, 91 and 136. However, the origin of these magic numbers can also be understood.

The magic number $`N=`$34 is not seen in the chain of the energetically most favourable icosahedral clusters, because of the lattice rearrangement that occurs for $`N=`$31. Beginning from $`LJ_{31}`$ a new cluster growing chain becomes energetically more favourable, but the chain of clusters growing from $`LJ_{30}`$ without its lattice rearrangement turns out to be rather close to the global energy minimum chain (see figure 10). As a result, in the vicinity of the transition point both chains influence the magic number formation. Thus, for $`N=`$34, the positive peak in the $`\mathrm{\Delta }^2E_N`$ for the chain of the global minimum clusters is absent, but it for present in the chain of clusters growing from the $`LJ_{30}`$ cluster without its lattice rearrangement (see open squares in figure 12). This magic number can also be seen for the octahedral and tetrahedral cluster chains, which come energetically rather close to the icosahedral chain in this region of $`N`$.

For the $`LJ_{81}`$ cluster the situation is similar. Here again a structural rearrangement of the cluster lattice takes place influencing the magic number formation. Thus, the magic numbers 81, 87, 91 arise in the chain of clusters growing from the $`LJ_{81}`$ cluster without its lattice rearrangement. The number 81 arises also in the chain of the decahedral clusters, which in this region of $`N`$ is energetically close to the global energy minimum (see figure 12). The chain of tetrahedral clusters affects the formation of the 91 magic number due to the same reason. Though this magic number is absent in the sequence of magic numbers generated for global energy minimum cluster structures, it is well pronounced in the magic number sequence for the tetrahedral clusters and the chain of clusters growing from $`LJ_{81}`$ without its lattice rearrangement.

Another magic number, which is masked in the magic number sequence for the icosahedral clusters is $`N=`$136. This number is masked because of the radical surface rearrangement of atoms needed for obtaining its configuration from the $`LJ_{135}`$ one. The $`LJ_{135}`$ cluster is the truncated 147 icosahedron. It possesses the $`I_h`$ point symmetry group. It can be obtained from the $`LJ_{134}`$ cluster by a rearrangement of 7 atoms located on the cluster surface. This cluster configuration can be obtained via a long chain of excited state cluster configurations which starts from the $`LJ_{55}`$ icosahedral cluster. In experiment, the formation of the global energy minimum $`LJ_{135}`$ cluster structure should occur with relatively low probability because of the reasons outlined above. Instead, it is feasible to assume that for this particular cluster size different cluster isomers with energies close to the global minimum can be readily created. This influences the 135 magic number formation resulting in its shift to $`N=136`$.

Note that $`\mathrm{\Delta }^2E_N`$ values calculated for the chain of icosahedral clusters in the size range $`136N146`$ are very close to zero. This fact can be explained as follows. There are 12 vacancies in the outer icosahedral shell of the $`LJ_{135}`$ cluster, which are located at the vertices of the $`LJ_{147}`$ icosahedron. These vacancies are sequentially filled within the size range $`136N146`$. Within this size range the energy difference between the two neighbouring clusters practically does not vary and is determined by the energy of an atom placed in one of the icosahedron vertice vacancies. Such a situation takes place, because the distances between any of two vertices are much larger than the pairing bond length in the LJ potential and the interaction between the atoms placed in the icosahedron vertices is neglible. As a result, of this the total energy grows nearly linearly and $`\mathrm{\Delta }^2E_N0`$

In figure 13, we plot images of the magic global energy minimum cluster structures. In this figure it is clearly seen how one icosahedral shell emerges from another. This figure demonstrates that the process of the $`LJ_{55}`$ icosahedron cluster formation has similarities with the formation of the $`LJ_{147}`$ cluster. In both cases, at the first stage a cap is formed on the top of the completed icosahedron shell, see the $`LJ_{19}`$ and $`LJ_{71}`$ magic clusters. Then, in both cases a structural rearrangement of the cluster takes place. When the rearrangement of the cluster lattice is completed ($`LJ_{31}`$, $`LJ_{82}`$), the growth of the cluster occurs by filling the atomic vacancies in the vicinity of the cluster surface up to the point when a new icosahedral shell is formed.

We have demonstrated that the magic numbers for the icosahedral clusters map well enough the experimentally observed magic numbers. In figure 9, we present the images of the non icosahedral global energy minimum cluster structures. Experimentally these clusters are not found to be the magic clusters, although, they are the global minimum clusters being magic for the chains of clusters based on the decahedral, octahedral or tetrahedral type of packing (see figures 10 and 12). This fact can be understood if one takes into account that the chains of clusters with different type of atomic packing are formed independently and the transition of clusters from one chain to another at certain $`N`$ occurs with the small probability. From the binding energy analysis, it is clear that the chain of icosahedral clusters prevails. Thus, in experiment, the icosahedral clusters mask the clusters belonging to the chains of non-icosahedraly packed clusters, even when the non-icosahedral cluster structures happen to become energetically more favourable.

### III.5 Spontaneous cluster fusion

Using the A2 method and the SE2 cluster selection criterion we have generated a chain of energetically unfavourable cluster isomers with $`N84`$ as an example of the spontaneous cluster growth.

In figure 14, we compare the binding energies per atom for the chain of the spontaneously produced clusters with the binding energies per atom calculated for the global energy minimum cluster structures. This figure shows that the binding energy per atom for the spontaneously generated cluster structures are systematically lower than those for the global energy minimum structures. However, the two dependences behave quite differently. Thus, at certain cluster sizes, the binding energy per atom for the spontaneously generated chain of clusters increases in a step-like manner. In the size range considered, this occurs for the cluster sizes $`N=1920,`$ $`3738`$, $`5051`$, $`6869`$, $`8182`$. These irregularities originate from the significant structural rearrangement taking place in the spontaneously generated cluster structures. In figure 15, we illustrate such a rearrangement for the clusters with $`N=1920`$ and $`N=6869`$.

The structural rearrangement in the $`LJ_{19}`$ cluster leads to the formation of the compact 13-atomic icosahedron shell. Indeed, in the initial $`LJ_{19}`$ cluster isomer, only a part of the icosahedron can be identified (see figure 15). Fusion of just one atom to the cluster boundary rearranges its structure and forms the compact icosahedral shell as a cluster core. This transformation stabilizes the cluster structure and makes it energetically much more favourable.

The situation with the $`LJ_{68}`$ cluster is similar. Again, we consider a highly asymmetric cluster isomer (see figure 15), which deviates significantly from the global energy minimum cluster structure based on the 55-atomic icosahedron shell (see figure 2). Fusion of a single atom to the boundary of the initial cluster configuration rearranges its structure completely forming a relatively compact $`LJ_{69}`$ cluster isomer (see figure 15). The $`LJ_{69}`$ cluster configuration possesses the evident elements of the 55-atomic icosahedron shell, which are absent in the initial $`LJ_{68}`$ cluster isomer. The compact configuration of the found $`LJ_{69}`$ cluster isomer brings its binding energy per atom close to the binding energy per atom calculated for the global energy minimum cluster structure.

### III.6 Symmetrical cluster fusion

Using the SY1-SY3 methods, we now consider the formation of symmetrical clusters and perform their energy analysis. We have generated the cluster configurations with $`N309`$ possessing the icosahedral and octahedral symmetry.

In figure 16, we plot the binding energies per atom calculated for clusters possessing the icosahedral and octahedral symmetry and compare them with the binding energies per atom calculated for the global energy minimum cluster structures. Figure 16 demonstrates the irregular behaviour of the binding energies of the calculated symmetrical cluster configurations. Note that this irregularities are not of the physical nature, because the SY1-SY3 methods do not model any concrete physical scenario of the cluster fusion process. These methods rather represent an efficient mathematical algorithm for the generation of symmetrical cluster configurations (see section II for details).

It is interesting that some of the found symmetrical cluster configurations within the size range $`N150`$ appear to be the global energy minimum cluster structures. This is the case for the clusters with $`N=`$13, 55, 135, 147, possessing the icosahedral point symmetry group, $`I_H`$. The binding energy of the highly symmetrical cluster isomer $`LJ_{115}`$ possessing also the $`I_H`$ point symmetry group is only about $`0.0033`$ units smaller than the energy of the corresponding global energy minimum. It is interesting that the $`LJ_{115}`$ global energy minimum cluster structure possesses the $`C_{5V}`$ point symmetry group, being lower than the $`I_H`$ point symmetry group. The fact that the cluster isomer with the lower symmetry becomes energetically more favourable can be understood after counting the total number of bonds in the system (see section III.1 for details). In the $`LJ_{115}`$ icosahedral isomer the total number of bonds is equal to 524, while in the global energy minimum isomer it is equal to 534. Here we have used the cut-off value equal to $`d=1.22^{1/6}\sigma `$. Comparison of the two number shows that the large number of bonds can be created in a system with the lower symmetry and that the larger number of bonds means the higher stability of the system.

Among the octahedral cluster configurations found with the use of the same technique there are also a few ones with the energy being very close to the cluster global energy minimum. Thus, the truncated octahedron $`LJ_{38}`$ appears to be the only global energy minimum cluster structure within the size range of $`N150`$ possessing the octahedral symmetry (see section III.1 for more details). The $`LJ_{116}`$ cluster isomer has also the form of truncated octahedron. Its binding energy is smaller than the binding energy of the global energy minimum cluster structure of the same size on $`0.0024`$ unit.

Note that some highly symmetric clusters, like the icosahedral $`LJ_{63}`$, $`LJ_{145}`$, $`LJ_{255}`$ can have relatively low binding energy as compared to their neighbours (see figure 16). The geometries of these clusters are shown in figure 17. They appear to be relatively low stable, because the total number of bonds in these clusters (222, 654 and 1194 respectively) is significantly lower than it is in the corresponding global energy minimum cluster structures (269 for $`LJ_{63}`$, 684 for $`LJ_{145}`$, for $`LJ_{255}`$ the total number of bonds in the global energy minimum is not known).

The calculation of the symmetrical cluster configuration and their energies turns out to be very useful for the preliminary analysis of the binding energy curve in a wide range. Thus, fitting the points $`N=13,55,115,135,147,267,297,309`$ with use of the liquid drop model equation (5), one derives the following values of $`\lambda _V=0.72948`$, $`\lambda _S=1.40469`$, $`\lambda _R=0.84143`$, which are in a good agreement with the results following from the global energy minimum calculations, see section III.3. We note that the performed analysis of symmetrical cluster configurations is essentially easier than the complete analysis of the global energy minima and thus it can be used for the efficient calculation of cluster binding energies in a wide range of $`N`$ at relatively high level of accuracy.

## IV Conclusion

In this paper we have discussed the classical models for the cluster growing process, but our ideas can be easily generalized on the quantum case and be applied to the cluster systems with different than LJ type of the inter-atomic interaction. It would be interesting to see to which extent the parameters of inter-atomic interaction can influence the cluster growing process and the corresponding sequence of magic numbers or whether the crystallization in the nuclear matter consisting of alpha particles and/or nucleons is possible. Studying cluster thermodynamic characteristics with the use of the developed technique is another interesting problem which is left open for future considerations.

## V Acknowledgments

The authors acknowledge support of this work by the Studienstiftung des deutschen Volkes, Alexander von Humbolt Foundation, DFG, Russian Foundation of Basic Research, Russian Academy of Sciences and the Royal Society of London.

## Appendix A Algorithm details

In this section we provide some details of the algorithms used in the computations NanoMeeting .

In the described algorithm, the potential energy of selected atom decreases during its motion. Therefore, the total potential energy of the cluster decreases at each step of the algorithm and it finally converges to a stable cluster configuration. In most cases, this convergence is rather fast. The time required for the obtaining of a stable configuration with the error $`\mathrm{\Delta }E`$ of the potential energy is usually proportional to $`ln\mathrm{\Delta }E`$. To speed up the process of conversion the procedure for the determination of the LJ force acting on the atom was written in Assembler x86. The performance of the calculations increased approximately 1.5 times comparing to the version written entirely in C.

Note that the described algorithm of the kinetic energy absorption is much more efficient as compared to those reducing the kinetic energy of moving atoms at each step of the Runge-Kutta integration procedure. Firstly, the following method was used: during the motion on each Runge-Kutta step the speed of an atom was decreased $`k`$ times. Therefore, the whole energy of oscillations that occur near the equilibrium position decreases because of the decreasing of the kinetic energy. So far, when the speed becomes very small, atoms start to oscillate very close to the potential energy minimum position. This method is rather simple, but too much time is required for atom coordinates to converge to their equilibrium position. Therefore we introduced another sufficient energy absorption algorithm, that was described in section II.

To prevent the fragmentation of the cluster caused by the unfavourable initial conditions, the adaptive choice of time-step in the Runge-Kutta integration procedure has been implemented. The fragmentation of the system might happen, for example, in the case when initially a pair of atoms are placed at a small distance one from another so that they experience a very strong repulsive force. In such a situation the moving atom gains high acceleration and during the first step of the Runge-Kutta procedure might go far away from the cluster. If the force acting on this atom in its final position is smaller than a given value $`F_{min}`$, it will be assumed in equilibrium and it will never be attracted back to the cluster, although the final configuration won’t be stable. The adaptive choice of the integration time-step prevents the atoms to fly apart. The idea of this procedure is rather simple. Each time step $`\mathrm{\Delta }t`$ in the Runge-Kutta integration procedure is repeated with the time step $`\mathrm{\Delta }t/2`$. If the coordinate or velocity change on the interval $`\mathrm{\Delta }t/2`$ is found too large then the integration step $`\mathrm{\Delta }t`$ is reduced and the procedure repeats. If the coordinate change is too small than the larger value of the interval $`\mathrm{\Delta }t`$ is set and the calculation is repeated again. The adaptive choice of the integration time step prevent the atoms to fly away, because in the unfavourable initial conditions situation the integration time-steps will become small, though the velocity will be high. Nevertheless, when the atom reaches the boundary of the cluster, it will stop: at that moment its velocity starts to decrease because of the attraction to the cluster.

## Appendix B Tables

## References


format: text

text:
# Chapter 1 Evaluation of Models for Diffuse Continuum Gamma Rays in EGRET Range

## Chapter 1 Evaluation of Models for Diffuse Continuum Gamma Rays in EGRET Range

### Abstract

The GALPROP model for cosmic-ray (CR) propagation produces explicit predictions for the angular distribution of Galactic diffuse $`\gamma `$-rays. We compare our current models with EGRET spectra for various regions of the sky. This allows a critical test of alternative hypotheses for the observed GeV excess. We show that a population of hard-spectrum $`\gamma `$-ray sources cannot be solely responsible for the excess since it also appears at high latitudes; on the other hand a hard CR electron spectrum model cannot explain the $`\gamma `$-ray excess in the inner Galaxy. Hence some combination of these explanations is suggested.

### 1. Introduction

Diffuse continuum $`\gamma `$-rays from the interstellar medium are potentially able to reveal much about the sources and propagation of CR, but in practice the exploitation of this well-known connection is problematic. We have previously  compared a range of models, based on our cosmic-ray propagation code GALPROP, with data from the Compton Gamma Ray Observatory. While it is rather easy to get agreement within a factor $``$ 2 from a few MeV to 10 GeV with a “conventional” set of parameters, the data quality warrant considerably better fits. Specifically, in order to reproduce the excess at GeV energies observed by EGRET, we adopted a model with a rather hard electron injection spectrum and a proton spectrum deviating moderately from that measured locally. An essential feature of our approach was that the locally-measured CR spectrum of electrons is not a good constraint because of the spatial fluctuations due to energy losses and the stochastic nature of the sources in space and time; the average interstellar electron spectrum responsible for $`\gamma `$-rays via IC emission can therefore be quite different from that measured locally. The hard electron spectrum interpretation is not entirely convincing however, since the electron fluctuations required are larger than expected. In addition, new accurate measurements of the local proton and Helium spectrum allow less freedom for deviations in the $`\pi ^0`$-decay component which were exploited in . Another suggestion which has been made  is that the $`\gamma `$-ray spectrum contains a $`\pi ^0`$-decay component from CR protons close to their (SNR) sources, so that in $`\gamma `$-rays we see the injection spectrum which is much flatter than the propagated spectrum which we measure locally.

In a companion paper  we use our model to re-determine the extragalactic $`\gamma `$-ray background (EGRB).

### 2. Data and method

We use the EGRET counts and exposure all-sky maps; the sources of the 3EG catalogue have been removed by the procedure described in . The predicted model skymaps are convolved with the EGRET point-spread function. The spectra are compared in the sky regions summarized in Table 1. Region A corresponds to the inner Galactic radian, region B is the Galactic plane excluding the inner radian, region C is the “outer Galaxy”, regions D and E cover higher latitudes at all longitudes, region F covers the Galactic poles. The models all use the locally-observed proton and Helium spectra. This is done because the nucleon data are now more precise (see references in ) than those which were available in . The nucleon injection spectra and the propagation parameters are chosen to reproduce the most recent measurements of primary and secondary nuclei, as described in detail in . The halo height is taken as 4 kpc as in , in accordance with our analysis of CR secondary/primary ratios .

### 3. Conventional model

We start by repeating the test of the “conventional” model which is based on the locally-observed electron (as well as nucleon) spectrum. The model spectra are compared with EGRET data in Fig. 1. As found in previous work, for such a model the GeV region shows an excess relative to that predicted; what is now evident is that this excess appears in all latitudes/longitude ranges. This already shows that the GeV excess is not a feature restricted to the Galactic ridge or the gas-related emission. Further it is clear that a simple upward rescaling of the $`\pi ^0`$-decay component will not improve the fit in any region, since the observed peak is at higher energies than the $`\pi ^0`$-decay peak. In the “SNR source” scenario  the spectrum in the inner Galaxy is attributed to an additional population of unresolved SNR, but this component cannot explain the excess at high latitudes, and hardly in the outer Galaxy. This explanation is therefore by itself insufficient, although as we will show it may well be part of the solution.

### 4. Optimized model

We next consider a model with a hard electron injection spectrum, as presented in . Our approach is to concentrate on obtaining a fit in all regions apart from the inner Galaxy, since then we can be reasonably sure that no population of unresolved sources is distorting the spectrum, and individual sources would be nearby, rare and easily identified, and so the situation should be relatively simple. (A population of weak sources in the halo could also contribute, but this is beyond the scope of this investigation). The result of our optimization procedure is shown in Fig. 2. An electron injection spectral index of 1.9 (with no break) is found optimal, consistent with what was found in . Values differing by up to 0.1 from this would also be acceptable. In order to be consistent with EGRET data $`>`$10 GeV, a cutoff in the electron spectrum at 3 TeV is required. Overall the fit is satisfactory and much improved over the conventional model and the GeV excess is reproduced in all regions except region A. At low latitudes in the inner Galaxy (region A) the peak around 1 GeV is not reproduced. Evidently a successful model requires both an adjustment of the IC via the electron spectrum and a hard-spectrum source population in the inner Galaxy. The sources must be concentrated in the inner Galaxy since in the remainder of the Galactic plane ($`30^{}<l<330^{},|b|<5^{}`$) the fit is satisfactory. Candidates for the source population include $`\gamma `$-ray pulsars like Geminga which does exhibit the required hard spectrum.

I.V.M. acknowledges partial support from a NASA Astrophysics Theory Program grant.

### 5. References

1. Berezhko E.G., Völk H.J. 2000, ApJ 540, 923

2. Moskalenko I.V., Strong A.W., Ormes J.F., Potgieter M.S. 2002, ApJ 565, 280

3. Strong A.W., Moskalenko I.V., Reimer O. 2000, ApJ 537, 763 (err. 541, 1109)

4. Strong A.W., Moskalenko I.V., Reimer O. 2003, these Proceedings


text:
# Stripes on a 6-Leg Hubbard Ladder

## Abstract

While DMRG calculations find stripes on doped n-leg t-J ladders, little is known about the possible formation of stripes on n-leg Hubbard ladders. Here we report results for a 7$`\times 6`$ Hubbard model with 4 holes. We find that a stripe forms for values of $`U/t`$ ranging from 6 to 20. For $`U/t34`$, the system exhibits the domain wall feature of a stripe, but the hole density is very broadened.

The nature of the ground state of the two dimensional t-J model continues to be controversial, even if one restricts one’s attention to numerical simulations. Density matrix renormalization group (DMRG) calculationsdmrg  for n-leg t-J ladders (with n ranging from 2 to 8) with J/t in the physical region of interest for the cuprates, find that stripes are formed when the ladders are doped with holeswhitescalapino . In contrast, the Green’s function Monte Carlo (GFMC) results of Sorella, et al., on square, periodic t-J lattices find a $`d_{x^2y^2}`$-wave superconducting ground statesorella . The aspect ratio of an n-leg lattice and the open boundary conditions used in the DMRG calculations may be sufficient to favor a striped state. Alternatively, the choice of the guiding trial wave function for the GFMC may bias the system towards a superconducting stateleecomment . It is clear that the t-J model is a delicately balanced system. Thus it is interesting to go back one step and study the underlying Hubbard model in which the onsite Coulomb interaction U is a parameter.

Much less is known about the nature of the ground state of the 2D Hubbard model, which is computationally more demanding. Previous DMRG calculationsbonca  on 3-leg Hubbard ladders found that diagonal three-hole stripes with a linear filling density of unity formed when U/t was greater than 5. Early mean-field Hartree-Fock calculations for the 2D Hubbard model found that vertical stripes were favored for U/t less than of order 4 and diagonal stripes formed at larger values of U/t schulz . These Hartree-Fock stripes have a linear filling density of one and their width is set by a coherence length $`\xi _0t/\mathrm{\Delta }_0`$ with $`\mathrm{\Delta }_0`$ the gap of the half-filled system. As U increases, $`\xi _0`$ decreases, approaching a lattice spacing at larger values of U/t. Eventually, when U exceeds twice the bandwidth, the mean-field stripes disappearinui . Note, however, that the Hartree-Fock treatments do not include pairing and that their energy per hole is of order $`1t`$ too high. Here, we extend the DMRG study of stripe formation in the Hubbard model and in particular examine the question of whether stripes form on a 6-leg Hubbard ladder and how their structure depends upon U/t. Our results represent the first reliable ground state results for a Hubbard cluster large enough to exhibit an unambiguous stripe. From this single cluster we are not able to make conclusions about the infinite 2D lattice; on the other hand, our results strongly suggest that stripes are low-lying states in the Hubbard model, as they are in the t-J model.

The cluster geometry was chosen with some care. Stripe-like behavior involving only two holes can be suggestiveprelovsek ; eder , but may only represent a feature of hole pairs. One would like at least four holes in a stripe, which makes the system size beyond the current reach of exact diagonalization. Periodic boundary conditions, when the dimensions of the cluster are even, frustrate a single stripeinui ; eder . Making one of the dimensions odd tends to force a single stripe. We choose a 7$`\times `$6 cluster, with cylindrical boundary conditions, periodic in the $`y`$ direction and open in the $`x`$ direction, and with four holes. Our previous calculations on the t-J model show that cylindrical L$`\times `$6 systems have ground states with four-hole stripes wrapped around the cylinder. The open boundaries neither frustrate nor force the stripes, although if they form, the boundaries serve to pin them. One may worry that Friedel oscillations in the charge density, induced by the open boundaries, mimic stripes, but we find, for a number of reasons, that our stripes in the t-J model are inconsistent with the Friedel oscillation scenariofriedel . Sorella et al. compared GFMC and DMRG on a 6$`\times `$6 cylindrical t-J cluster with 6 holessorella ; however, we expect 6 holes to split into a stripe and a pair, with resonance between them, obscuring any obvious signs of a stripe in the charge or spin densities. (A stripe having an odd number of holes would be frustrated in the $`y`$ direction.)

Here we begin by looking at the t-J model and then turn to the Hubbard model. A typical result for a 17$`\times 6`$ t-J lattice with J/t $`=`$ 0.35 and 12 holes, with cylindrical boundary conditions, is shown in Fig. 1. A maximum of $`m=5000`$ states were kept, for a discarded weight of $`7\times 10^5`$. One can see in Figure 1 that three charged stripes have formed and the spins form $`\pi `$-phase-shifted regions between the stripes. For an L$`\times `$6 ladder, the linear charge density along a stripe is 2/3. Studies of longer stripes find that the preferred density on a long stripe is 1/2. Competing with the four-hole stripe, but higher in energy per hole, is a single pair, which also exhibits some stripe-like features. As seen in Fig. 1, there are two types of four-hole stripes: the stripes on the ends are (mostly) bond-centered, and the stripe in the middle is site-centered. The energy per unit length of these two configurations is extremely close. The charge and spin density profiles along the x-direction are plotted in a more conventional fashion in Fig. 1(b). Since later we will be considering a 7$`\times `$6 ladder for the Hubbard model, also shown is the local hole density on a 7$`\times `$6, t-J ladder with 4 holes and J/t $`=`$ 0.35. One can see that the open boundaries on the small system mimic the presence of the other stripes on the larger system.

In Fig. 1, the presence of static magnetic moments makes the $`\pi `$-phase shift in the antiferromagnetic spin density as one crosses a stripe clearly visible. No external magnetic fields were applied to the system, and these static magnetic moments are artifacts of the DMRG calculation, for which no spin symmetry has been used. However, we argue that the presence of these static moments has very little effect on the validity of the results. To see this, we show DMRG results for an undoped 17$`\times `$6 system in Fig. 2. The undoped system converges very rapidly with the number of states $`m`$, allowing us to study the convergence in detail. In particular, we consider $`S_z(l)`$ versus the error in energy for a particular number of states kept $`m`$. The exact ground state of this system, by the Lieb-Mattis theorem, is a spin singlet, which is spin-rotationally invariant, and therefore $`S_z(l)=0`$ for all sites $`l`$. However, as measured by the spin-spin correlation function, this state has long-range antiferromagnetic order. The ground state is a superposition of all the static antiferromagnetic states with all possible orientations of the order parameter. The reduction in energy due to the macroscopic superposition of the states is very small. We see from Fig. 2 that there exist states with static moments of order 0.2 with energies only 0.0002$`t`$ per site above the ground state. In a doped system with some form of magnetic order, the behavior of the holes in each of the states with different spin orientations would be identical. It is reasonable to assume that the effect of the superposition on the hole behavior is slight.

We can also analyze this effect in the undoped system by applying a staggered magnetic field. We assume the energy per site of the system, for small applied fields, varies as

$$E(s,h)=E_0+as^2sh$$
(1)
where $`h`$ is the magnitude of the applied staggered field in the $`z`$ direction and $`s`$ is the average magnitude of $`S_z(l)`$ in response. We minimize $`E(s,h)`$ over $`s`$, keeping $`h`$ fixed, to find $`s(h)=h/(2a)`$, and hence $`E(s)=E_0as^2`$. From several very accurate DMRG simulations with various values of $`h`$, we find that for $`s<0.10.2`$, this energy dependence is accurate, and for the 17$`\times `$6 system $`a0.0034`$. To describe the case where there is no applied field, and a finite value of $`s`$ is considered an error, we define $`\mathrm{\Delta }E=as^2`$, yielding $`s(\mathrm{\Delta }E)=(\mathrm{\Delta }E/a)^{1/2}`$, which is shown as the solid line in Fig. 2. This should be considered to be an upper bound on the $`s(\mathrm{\Delta }E)`$ obtained from DMRG, for which there are other sources of error in the wavefunction besides a finite value of $`s`$.

We turn now to the Hubbard model. Our results for the Hubbard systems were performed using a new “single-site” DMRG methodsinglesite , which performs better than the standard DMRG algorithm having two sites in the center when the number of states per site is more than two or three. Using this method, we have been able to keep up to 7500 states per block in some cases. Unfortunately, discarded weights are not informative and we do not report them. When performing DMRG calculations on 2D clusters, one must deal with the possibility that the calculation will get stuck in a metastable state. For example, a striped state may be lower in energy than a state with two widely separated pairs, but the calculation may only be able to tunnel between these configurations when keeping very large numbers of states per block. In this case, one must repeat the calculations with constrained initial configurations, and compare final energies. On the other hand, a calculation is particularly robust if one sees that it does tunnel between very different states. Fig. 3(a) shows a plot of the ground state energy for a 7$`\times `$6 lattice with U/t $`=`$ 12 and 4 holes as a function of the number of basis states $`m`$. The hole density distribution obtained at a number of sweeps, labeled by the value of $`m`$, is shown in Fig. 3(b). The initial configuration consisted of two separate pairs. As the number of basis states increases and the ground state energy converges, one clearly sees the stripe develop, with the “tunneling” occuring for about $`m=1200`$ states. From Fig. 4 we see that, just as for the t-J model, there is a $`\pi `$-phase shift in the magnetization density across a stripe. One can see that as the stripe develops, the DMRG ground state energy decreases and that just as for the t-J model, doped holes on 6-leg Hubbard ladders can form striped ground state structures for U/t $`=`$ 12.

We have performed a limited study of the behavior of this system as a function of U/t. Fig. 5 shows the charge and spin densities for U/t ranging from 3 to 20. For U/t$`=`$8, we started the system with the holes as two separate pairs. We observed tunneling to the stripe state near $`m=3600`$, considerably later than for U/t$`=`$12. We let this calculation continue until it used all the available memory (3 Gb), taking about 1 week of computer time, on an Athlon MP 1800+ processor. In this case, further sweeps reached a maximum of $`m=7500`$, with little change in the charge and spin distributions. For U/t$`=`$6, we observed tunneling to the stripe state near $`m=4000`$, again with a maximum of $`m=7500`$. For U/t$`=`$4 and U/t$`=`$3, we did not try to observe the tunneling, instead starting with the four holes together in the center. Here, we found that a broadened striped configuration was stable, up to the maximum number of states kept ($`m=6000`$ and $`m=4000`$, respectively). For U/t$`=`$3, the charge distribution is perhaps too broad to consider it a stripe, but the domain wall nature of the spin configurations configurations is still fairly robust. However, since, as noted above, the spin expectation values are (useful) artifacts of the DMRG procedure, we do not attach much significance to the specific magnitudes shown in Fig. 5(b). Fig. 5 also shows results for U/t$`=`$20. For this run, we started it in a state with the four holes together, but during the first several sweeps, keeping only $`m200`$, the holes partially split apart. Subsequently, near $`m=1000`$, a definite stripe formed along with the antiferromagnetic domain wall. A maximum of $`m=6000`$ states were kept. For this large value of U/t, the stripe appears somewhat less stable than at U/t$`=`$12. We observe a very similar weakening in the stripe in the t-J model, as the value of J/t is reduced from 0.35 to 0.2 (the peak height in $`n_h(l)`$ drops from 0.21 to 0.18).

In conclusion, we have found that on a moderately sized Hubbard cluster, with cylindrical boundary conditions and doped with four holes, the ground state has a stripe. The stripe is narrow and well-defined for U/t$`=`$8-12. For smaller values of U, starting at U/t$``$6, the stripe broadens, until at U/t$`=`$3 the width of the stripe is the size of the system. At U/t$`=`$20, the stripe is somewhat broadened compared to U/t$`=`$12. The overall behavior is very similar to that seen in the t-J model. Although the open ends in one direction may encourage the stripe to form, the same roll may be played by neighboring stripes in larger systems. We find nearly the same hole density profile for a stripe in a 7$`\times `$6 t-J system as in the central portion of a t-J 17$`\times `$6 t-J system. For small U/t, the broad stripes we obtain are probably strongly influenced by the open boundaries.

The would like to that A. Chernyshev for useful comments. We acknowledge the support of the NSF under grants DMR98-70930 and DMR98-17242.