Joseph J. Smulsky
Institute of Earth’s Cryosphere, Tyum. SC of SB RAS, Federal Research Center, Tyumen, Russia.
DOI: 10.4236/jmp.2024.158051   PDF    HTML   XML   192 Downloads   699 Views   Citations

Abstract

Usually, models of globular star clusters are created by analyzing their luminosity and other observation parameters. The goal of this work is to create stable models of globular clusters based on the laws of mechanics. It is necessary to set the coordinates, velocities and masses of the stars so that as a result of their gravitational interaction the globular cluster is not destroyed. This is not an easy task, and it has been solved in this paper. Using an exact solution of the axisymmetric gravitational interaction of N-bodies, single-layer spherical structures were created. They are combined into multilayer models of globular clusters. An algorithm and a program for their creation is described. As a result of solving the problem of gravitational interaction of N bodies, evolution of 5-, 10-, and 15-layer structures was studied. During the inter-body interaction, there proceeds a transition from the initial specially organized structure to a structure with bodies, uniformly distributed in space. The number of inter-body collisions decreases, and the globular cluster model passes into the stable form of its existence. The collisions of bodies and the acquisition of rotational motion and thermal energy by them are considered. As a result of the passage to scaled dimensions, the results were recalculated to the conditions of globular star clusters. The periods of rotation and the temperatures of merged stars are calculated. Attention is paid to a decreased central-body mass in the analyzed models of globular star clusters.

Keywords

N-Body Problem, Solution, Globular Star Clusters, Properties

Share and Cite:

1. Introduction

Globular star clusters are spherical or somewhat flattened formations that contain from tens of thousands to a million stars [1]-[6]. The average density of substance in globular clusters is 0.4 stars per cubic parsec (pc) [7]. In the center, the density increases up to 100 - 1000 stars/pc3. The diameter of clusters can reach tens of parsecs [5]. Globular star clusters are common objects in the Universe: there are about 150 of them in our Galaxy alone [2]. They are old formations. The age of many of them exceeds 10 billion years (Gyr) [2] [5] and even exceeds 13.8 Gyr [7], i.e., the age of the Universe assumed in the “Big Bang” theory. Therefore, in recent decades, astronomers have reduced the age of globular clusters. However, after the James Webb Space Telescope discovered distant galaxies with an age exceeding 13.8 Gyr, the age of very old globular clusters is now raised to 26.7 Gyr [8]. It should be noted that relatively young globular clusters are observed in the Magellanic clouds, satellites of our Galaxy, and in other galaxies [5] [9]-[11].

There also exist dwarf spheroidal galaxies. Such galaxies are predominantly satellites of other galaxies, but they also occur as isolated objects [12]. The nuclei of galaxies and their surrounding halos also have a spherical shape. Thus, the spherical shape of stellar associations is often found in the Universe. Therefore, researchers believe that uncovering the mechanisms behind the existence of globular clusters will provide insight into processes in denser systems such as galactic nuclei [3].

Various methods are used to model stellar associations such as globular star clusters and galaxies [5] [3] [13]. In some of them, the entire cluster zone is considered as a continuous medium, and in others, as a set of objects with random kinematic characteristics. The latter models are not deterministic. In deterministic models, each body has its own size, mass, coordinates, and velocity. The gravitational interaction of each such body with any other body is investigated. Therefore, the position and velocity of any body are known at any time.

Depending on the method of solving the N-body problem, such models can differ significantly [4] [5]. For example, computer programs implementing such methods generally fall into four groups [6]. The most common programs for solving the N-body problem are the programs NBODY 2 - NBODY 6 [10] [14] developed by S.J. Aarseth [13]. Those programs, such as NBODY 6 [6], use the Hermite method based on the expansion of coordinates and velocity of bodies in Taylor series. This approach uses three derivatives [6]. In the process of solving the N-body problem, bodies can approach each other. At such moments, a step division algorithm is activated, and at small distances between bodies the calculation is terminated. Therefore, various regularization methods are introduced [4] [6] [13]. When two bodies approach together, as a rule, switching to the algorithm for the exact solution of the two-body problem is used. When three bodies approach each other, other algorithms are employed. Such regularization avoids slowing down the calculation process, but simultaneously leads to a loss of overall solution accuracy. That is why an additional effect is introduced, namely the tidal one [6] [13]. This effect is also added to take into account the influence due to objects located outside the globular cluster.

In these problems, one of the main tasks is the assignment of initial conditions, in particular, coordinates and velocities of bodies. For this purpose, many different methods were used [13]. In the most common of them, based on the brightness distribution and the number of stars in a cluster, the initial mass function (IMF) of the cluster [5] [14] or the initial concentration function [15] [16] are calculated. With their help, the velocity of stars at a certain radius from the cluster center can be found. The stars are then made statistically distributed in space. It should be noted that this leads to a violation of the determinism of the problem.

Since the positions and velocities of bodies obtained in this way can lead to the destruction of the globular cluster model, they must be further corrected and adapted [4] to avoid undesirable effects.

In many works on solving the N-body problem it was noted that there is only one exact solution to this problem, that for the two-body problem [4] [13]. An exact solution of the N-body problem is required for testing computer programs for numerical solution [13].

In fact, such solutions do exist. We have obtained some of them [17], for example, one solution for N bodies located axisymmetrically on a plane [18]. We should also mention a second solution, in which there are N2 concentric layers around a central body, each layer containing N3 bodies [19]. In this case, the entire structure rotates around its axis as a whole and the masses of one body in each layer are different.

In the present paper, globular cluster models are created based on the solution of the first problem [20] [21]. The solution to this problem exists for any number N of bodies. In this case, the central body may be absent. It should be noted that this planar axisymmetric N-body problem was also solved for Coulomb interaction [22]. Its solution exists only in the presence of a central particle, i.e. positively charged nucleus, and with the number of peripheral particles, i.e. electrons, N3 < 473.

Thus, by solving the planar axisymmetric N-body problem, the masses of the central and peripheral bodies, and their coordinates and velocities, are determined. As a result of research [20] [21], it was found that with a uniform arrangement of bodies in space, they move along the same orbits as on the plane. It should be noted here that those orbits can be a circle, an ellipse, a parabola, or a hyperbola. In the present work, globular clusters with circular orbits of bodies were treated. In such clusters, all bodies have certain masses, coordinates and velocities, that is, the cluster is strictly deterministic, with its motion being strictly pre-defined in the future.

For solving the N-body problem, we have developed a system called Galactica. This system was used to solve a number of problems in celestial and cosmic dynamics [23] [24]. The guide for working with Galactica is available in [25], the text of the program for gravitational interaction is published in [24], and the entire system is freely available¹. The Galactica system is a universal program. It can be used for calculating Coulomb interactions in the micro-word [22]. A guide to working with the Galactica system in this case is published in [26].

The Galactica program uses a highly accurate solution method. Like in the NBODY 6 program, the coordinates and velocities of bodies in it are expanded into Taylor series, with derivatives up to the sixth order inclusive taken into account. Due to this, the accuracy in solving Solar-system dynamics problems is several orders of magnitude higher than, for example, the accuracy offered by NASA programs [27] [28].

When solving the problem of evolution of the Solar system, a number of authors have arrived at a conclusion about Solar-system instability, which will manifest itself after 20 million years [29] [30]. Using the Galactica system, the problem of evolution of the Solar system over a period of 100 million years was solved, and it was shown that, on the contrary, the Solar system will remain stable, with no tendency to its changes observed [31] [32]. This is due to the greater accuracy of the Galactica program.

The Galactica program uses no regularization when bodies approach each other. At each step, accuracy is controlled and, if necessary, the step size is changed. When the bodies come in contact, they combine into one body, whose kinematic parameters can be determined as governed by the mechanics laws. Therefore, the resultant bodies have angular momentum, i.e. spins, and thermal energies. From these quantities, one can calculate the period of rotation of the body and its temperature. It should be noted here that Galactica has a number of functions that allow one to analyze and study the results obtained when solving the N-body problem.

In this way, the interactions of bodies in a spherical single-layer structure were analyzed [20] [21]. The construction of such a structure was substantiated, its development in the process of inter-body interaction was demonstrated, and evolution was studied. Such a structure presents a model of a globular star cluster with a rarefied core. In the present work, multilayer structures with each layer created according to the single-layer structure construction algorithm have been considered. The algorithm for constructing such structures is substantiated, issues concerning the choice of parameter values for these structures are discussed, and the evolution of several models of multilayer structures is studied. During this work, various phenomena and properties of modeled structures were revealed, with considerable attention having been given to the study of modeled objects. Such phenomena also occur in globular star clusters.

The structure of the paper is as follows. The method for constructing globular clusters is described at the beginning. The evolution of 5, 10, and 15-layer globular cluster models is then considered. Then their general properties are described. These studies were carried out in dimensionless form. Further, the results obtained are presented in dimensional form in the scale of a globular cluster. At the end of the paper, models of the central body are studied, which makes it possible to reduce its mass by tens of times.

2. Basic Principles in Constructing Multilayer Structures

A layered structure is created from a number of structures distributed over a sphere [20] [21]. The single-layer structure is based on solving the problem of interaction of N₃ peripheral bodies with a mass m_i = m₁ located axisymmetrically on a plane around a central body of mass m₀ [18] [33] [34]. Such a structure unfolds in space as follows. The second peripheral body m₂ (Figure 1(a)), together with the other bodies m₃, m₄, …, m_i, …, m_N3, is rotated in the x₁y₁-plane through an angle Δψ counted from the position of the first body. Then, through the x₂-axis passing through the second body, the entire plane with the bodies m₂ - m_N3 is rotated through an angle Δθ from the initial plane. The velocity vectors of these bodies are also located in a new plane. Such a process of rotations through the angles Δψ and Δθ is performed sequentially for all bodies from m₃ to m_N3. As a result, a structure distributed in space is formed (Figure 2(a)), in which the coordinates and velocities of bodies during their interaction allow them to move in space with periods almost the same as those on the plane. By changing the increments Δψ and Δθ, one can create various spherical structures.

There are also a number of other possibilities that allow one to create various structures using rotations through angles Δψ and Δθ. In the works by Smulsky [20] [21], four such possibilities were analyzed. In the fourth case, not the bodies m_i themselves but the pericenters of their orbits P_c,i are located on the initial circle (Figure 1(b)). The orbit can be an ellipse, a parabola, or a hyperbola. Peripheral bodies have similar orbits with eccentricity e, each body being located at an angular distance φ_0,i from its pericenter P_c,i. In this case, the bodies in the layer will be located in a 3D region remote from the pericenter radius R_p to the apocenter radius R_a.

Figure 1. Geometric and kinematic characteristics of a single-layer axisymmetric structure comprising N₃ bodies, with a central body of mass m₀ and peripheral bodies of mass m_i = m₁: (a) rotations of bodies and their velocities through angles Δψ and Δθ; (b) polar coordinates r_i and φ_0,i of the peripheral body m_i over a trajectory section; the polar angle φ_0,i is reckoned from the orbit pericenter P_c,i.

This option was used for creating multilayer structures. However, in the study reported in this paper there was no need to use elliptical orbits, therefore all treated structures were created with zero eccentricity.

Let there be N₂ layers in a multilayer structure, enumerated with numbers j = 1, 2, …, N₂. Consider the coordinates and velocities of a body in layer j at its initial on-plane position. The bodies on the ring of layer j are evenly spaced with an interval $\Delta {\phi }_j=2\pi /N_{3,j}$ . Those bodies move along orbits with an eccentricity e. They are located at an angular distance from the pericenters of their orbits (Figure 1(b)), which are reckoned from the axis x_p,j,i

${\phi }_{0,j,i}=\left(i_j-1\right)\cdot \Delta {\phi }_j,i_j=1,2,\cdots ,N_{3,j}$ .(1)

Note that, in contrast to Figure 1(b), in formula (1) and in subsequent formulas, an additional subscript indicating the layer number j is introduced.

From the trajectory equation of the peripheral body [18] [33], one can determine the distance of bodies r_j,i from the origin O, their radial velocity v_r,j,i, and the transversal velocity v_t,j,i. Then, in the coordinate system x_p,j,iy_p,j,iz_p,j,j with the x_p,j,i-axes passing through the pericenter P_c,j,i (Figure 1(b)), the coordinates and velocities of the peripheral bodies can be written as

$x_{p,j,i}=r_{j,i}\cdot \cos {\phi }_{0,j,i};y_{p,j,i}=r_{j,i}\cdot \sin {\phi }_{0,j,i};z_{p,j,i}=0$ ;(2)

$\begin{array}{l}{v}_{xp,j,i}=v_{r,j,i}\cdot \cos {\phi }_{0,j,i}-v_{t,j,i}\cdot \sin {\phi }_{0,j,i};\\ v_{yp,j,i}=v_{r,j,i}\cdot \sin {\phi }_{0,j,i}+v_{t,j,i}\cdot \cos {\phi }_{0,j,i};v_{zp,j,i}=0.\end{array}$ (3)

As a result of solving the problem of gravitational interaction of bodies in an axisymmetric structure [18] [33] [34], the trajectory equation for a peripheral body in the polar coordinate system r_j,i(φ_0,j,i) is obtained in the form

$r_{j,i}=\frac{R_{p,j}}{\left({\alpha }_{1,j}+1\right)\cdot \text{cos}{\phi }_{0,j,i}-{\alpha }_{1,j}}$ ,(4)

where R_p,j is the pericenter radius, i.e. the radius of the point on the orbit with the smallest distance to the origin O in Figure 1(b);

${\alpha }_{1,j}={\mu }_{1,j}/\left(R_{p,j}\cdot {\upsilon }_{p,j}^2\right)$ ;(5)

${\mu }_{1,j}=-G\left(m_{j,0}+m_{j,1}\cdot f_{N3,j}\right)$ ;(6)

$f_{N3,j}=0.25{\displaystyle \sum _{i_j=2}^{N_{3,j}}\frac{1}{\sin \left[\pi \left(i_j-1\right)/N_{3,j}\right]}}$ .(7)

In formulas (5)-(7), the following designations are used: α_1,j is the trajectory parameter; _1,j is the interaction parameter; and f_N3,j is the contribution due to the action of N_3,j − 1 peripheral bodies on one of the bodies.

Depending on the value of trajectory parameter α₁, the orbits of peripheral bodies can be circles (α₁ = −1), ellipses (−1 < α₁ < −0.5), parabolas (−1 < α₁ < −0.5), or hyperbolas (−0.5 < α₁ < 0). The time of motion of a body along the trajectory also depends on α₁ [18] [33] [34].

Below, four other parameters of the orbit of peripheral bodies will be needed [33] [34]: the orbital period

$P_j=-\frac{2\pi {\alpha }_{1,j}\cdot R_{p,j}}{v_{p,j}{\left(-2{\alpha }_{1,j}-1\right)}^{3/2}},$ (8)

the velocity at the pericenter

$v_{p,j}=\sqrt{{\mu }_{1,j}/\left({\alpha }_{1,j}\cdot R_{p,j}\right)}$ ,(9)

the eccentricity of the orbit

$e_j=-\left(1+1/{\alpha }_{1,j}\right)$ ,(10)

and the semi-major axis of the orbit

$a_j=R_{p,j}\left(2{\alpha }_{1,j}+1\right)/{\alpha }_{1,j}.$ (11)

The radial velocity of a peripheral body is given by [33] [34]

$v_{r,j,i}=\pm v_{p,j}\sqrt{{\left({\alpha }_{1,j}+1\right)}^2-{\left({\alpha }_{1,j}+R_{p,j}/r_{j,i}\right)}^2}$ .(12)

The radial velocity is positive when the body moves from pericenter to apocenter, and negative when returning back. The transversal velocity is written as [33]:

$v_{t,j.i}=v_{p,j}\cdot R_{p,j}/r_{j,i}$ (13)

Providing that the masses of the bodies m_j,0 and m_j,1, the pericenter radii R_p,j, and the trajectory parameters α_1,j or the eccentricities e_j are given, expressions (1)-(13) determine the coordinates and velocities of peripheral bodies in all N₂ layers of the flat structure.

In order to be able to vary the structures, the coefficient k_φ of the angle between the bodies in the layer and the coefficient k_φv of the angle of rotation of the velocity vector are introduced, with the help of which these angles can be expressed as

$\Delta \psi =k_{\phi }\cdot \Delta \phi$ ; $\Delta \theta =k_{\phi v}\cdot \Delta \phi$ ,(14)

where $\Delta {\phi }_j=2\pi /N_{3,j}$ .

This algorithm will be used to create a multilayer structure consisting of N₂ layers. The orbit of the bodies in a layer is given by the eccentricity e and the major semi-axis a. Using expression (10), the trajectory parameter α_1,j is determined by the eccentricity e, and according to formula (11), the pericenter radius R_p,j can be expressed in terms of the semi-axis of the orbit a_j as follows:

$R_{p,j}=a_j\cdot {\alpha }_{1,j}/\left(2{\alpha }_{1,j}+1\right)$ .(15)

It was shown [33] that with a centrally symmetrical arrangement of bodies in space, the force of the outer layer on the mass inside is zero. On the other hand, the force of action of such a structure on a mass located outside is equal to the force of action of a material point located in the center of the structure and having a mass equal to the mass of this structure. Therefore, for each layer, starting from j = 2, we assume that its center has a central body with mass m_j,0 equal to the mass of all bodies inside the layer j, i.e.

$m_{j,0}=m_{j-1,0}+N_{3,j-1}{m}_{j-1,1}$ .(16)

As already noted, a flat multilayer structure described by expressions (1)-(16) turns into a spatial one following successive rotations through angles Δψ and Δθ. These rotations are to be executed for each layer. The algorithm for applying these rotations was presented in [20] [21] as variant 4. Therefore, we do not reproduce it here.

In addition to this algorithm, three more operations have been added. The first operation is related with the refinement of the distance between the bodies in one layer (Figure 2(a)). The first peripheral body 1 is located on the x-axis to be followed by bodies 2, 3, 4, and so on, located at almost identical distances from each other. The last body 99 is not at the same distance from the first body. An algorithm for refining the coefficient k_φ in formula (14) was developed. This algorithm allows one to calculate the angle Δψ between the bodies so that the distance between the first and the last body be the same as between all other bodies. For this purpose, the coordinates of the body with number N_3,j + 1 are calculated. If these coordinates coincide with the coordinates of the first body, then the distance between the bodies N_3,j and 1 will be the same as the distance between the other bodies. The calculation is performed by the method of successive approximations until a relative change in the coefficient k_φ of the specified accuracy EPS is reached. As a rule, the coefficient k_φ for six iterations is determined with a relative accuracy EPS = 1 × 10⁻¹⁰.

The second operation is related to the correction of the position of bodies at the points of self-intersection of their initial location line (Figure 2(b)). For this purpose, the number N₃ of bodies in the layer is calculated so that the body on the line crossing the gap between two bodies is approximately at the same distance from them. We denote the average distance between bodies as d_m, and the minimum distance between the bodies i₁ and i₂ as d_mn, where the bodies i₁ and i₂ are located on different parts of the initial location line. Then, the length of this line is

Figure 2. Arrangement of 99 bodies in one layer at k_φ= 1.7 and k_φv = 1 in the xyz coordinate system: (a) in the uncorrected layer, the distance between the 99^th body and body 1 differs from the distances between other bodies; the segments show the velocity vectors of the 1^st and 99^th body; (b) correction of distances in the regions of self-intersection of the layer; (c) rotation of the layer through an angle β.

$l_0=N_3{d}_m$ ,(17)

and the distance from the x-axis, on which the body with number i = 1 is located, to the body with number i₁ is

$l_1=\left(i_1-1\right)d_m$ .(18)

The distance to the body with number i₂ will be greater by d_mn:

$l_2=\left(i_1-1\right)d_m+d_{mn}$ .(19)

Suppose that, with a new number of bodies N₃^’, the body i₂ is located in the middle between the bodies i₁ and i₁ + 1, so that the distance to this body is

${l^{\prime }}_2=\left(i_1-1\right){d^{\prime }}_m+0.5{d^{\prime }}_m={d^{\prime }}_m\left(i_1-0.5\right)$ .(20)

Since the length of the initial location line l₀ remains unchanged, then ${l^{\prime }}_2$ equals l₂; from here, taking into account formula (17) for ${d^{\prime }}_m$ , we obtain a new number of bodies in the layer:

${N^{\prime }}_3=\frac{N_3\left(i_1-0.5\right)}{i_1-1+d_{mn}/d_m}$ .(21)

The new amount of bodies ${N^{\prime }}_3$ calculated by formula (21) must be rounded up to the nearest integer number. As a rule, the new number of bodies thus obtained makes it possible to increase the minimum distance d_mn between the bodies i₁ and i₂ (Figure 2(b)) to acceptable values.

The third operation consists in the uniform rotation of the layers, starting from the second layer, in clockwise direction around the x-axis through an angle β_j (Figure 2(c)). The coordinate axes of the rotated layer are denoted as y_rtz_rt,. The projections onto the axes of the yz coordinate system are

$y=y_{rt}\cos \beta +z_{rt}\sin \beta$ ; $z=-y_{rt}\sin \beta +z_{rt}\cos \beta$ .(22)

The angle of rotation of the j-th layer is given by the expression

${\beta }_j=\left(j-1\right)\cdot \Delta \beta$ ,(23)

where $\Delta \beta =2\pi /N_2$ .

3. Software for Creating Multilayer Structures

For performing multiple calculations when creating a multilayer spherical structure, an MLSpStr2.for program has been developed. This program consists of three parts implementing the following operations: 1) reading the initial parameters; 2) construction of the multilayer structure; 3) creating a file with initial conditions for the Galactica system.

In the present paper, a dimensional/non-dimensional method of treating data is used. When setting values of structural parameters, the Solar system is used as an analogue providing parameter ratios. After the structure is created, its dynamics and evolution are analyzed in dimensionless form. The results are discussed in dimensional form with parameter values inherent to globular star clusters.

The main initial parameters are read from the data file MLSpStr2.dat. In this file, the following structural parameters are specified: N₂ is the number of layers; N₃₀ is the initial number of peripheral bodies in the first layer; mi is the initial mass of the central body and all first-layer bodies; p_m0 is the fraction of the mass mi due to the central body; A_Sm is the semi-axis of the orbits of first-layer peripheral bodies in astronomical units (AU); e is the eccentricity of the orbits of peripheral bodies; k_a is the coefficient of the semi-axes of the layers, starting from the second layer; k_N3 is the coefficient of the number of bodies on these layers; k_φ0 and k_φv are the coefficients of the initial angles of bodies and their velocities when constructing the structure; EPS is the allowed error in calculating k_φ for the body’s angle; I_nx is the key for initiating the uniform rotation of layers around the x-axis; I_cm is the key for issuing the coordinates and velocities of peripheral bodies in the information file MlSpStr2Err; and ρ_b is the absolute density of bodies. Note that the density of bodies ρ_b expressed in kg/m³ is necessary for calculating their radii. The radii of bodies are used in the Galactica program when calculating the inter-body collisions. The MLSpStr2.dat file also specifies a number of other parameters required for Galactica.

The above-listed data completely determine the parameters of the central body and first-layer bodies. The semi-axis of the orbits of the rest layers is calculated as

$a_j=A_{sm}\cdot AU\cdot \left(1+k_a\left(j-1\right)\right)$ ,(24)

where AU is the astronomical unit, and the number of bodies in a layer is determined in proportion to the semi-axis of the orbit:

$N_{3,j}=N_{3,1}\cdot k_{N3}\cdot \left(a_j/a_1\right)$ .(25)

The masses of peripheral bodies are identical in each layer. This mass is calculated from the difference between the mass m_i and the central-body mass, and from the number of bodies N_3,1.

After reading the initial parameters of the structure, the algorithm presented above calculates the coordinates and velocities of all bodies involved. Calculation results are issued in the form of three output files: fN3fvout.dat, MlSpStr2Err and, for example, MS15c49b.dat. The file fN3fvout.dat contains the layer numbers j, the number N_3,j of bodies, the coefficients k_φ,j and k_φv,j, the semi-axes a_j expressed in meters, and the mass m_j of one peripheral body in kg.

In the MlSpStr2Err file, for each layer, the period P_j in sidereal years, the number of bodies N_3,j, the average distance between bodies d_m in meters, the minimum distance d_mn between the bodies i₁ and i₂ at the point of intersection of their location line with indication of their numbers i₁ and i₂ are output. The distance between bodies 2 and 3 is output as d_m. In addition, the number of iterations in the calculation of coefficient k_φ,j and the two last values of this coefficient are indicated.

Then, the summary information for all layers is given, including the minimum distance between the bodies in the entire structure with indication of the layer and body numbers; the dimensionless coordinates and velocities of the center of mass of the entire system, and the first body in the first layer. For the center of mass, these values should be zero, and their non-zero values indicate the level of error in creating the system. For example, for a fifteen-layer structure MS15c49b.dat with the number of bodies N = 5866, the values of interest are at the level of 2 × 10⁻¹⁷ when the program is compiled with double precision, i.e. when the number length is 16 significant digits. This level of error indicates that the error is in fact extremely low.

At the end of this file, the initial data specified in the MLSpStr2.dat file are output. Thus, the MlSpStr2Err file is a kind of a passport of the created structure: it contains all necessary information about it.

Additionally, in the MlSpStr2Err file the coordinates and velocities of all bodies in the layers in dimensional form are output when the key Icm in the initial data file MLSpStr2.dat is set to 1. The MlSpStr2Err file is also intended for issuing error messages when the MLSpStr2.for program is running. When an error occurs, its decryption is written to this file.

Based on the initial data specified in the MLSpStr2.dat file, a structure is created, with the desired dimensions of layers and the number of bodies in them, according to the algorithm (24)-(25). In creating structures with a different algorithm, an additional source data file fN3fvinp.dat is used. The latter file specifies the layer numbers j, the number of bodies N_3,j, the coefficients k_φ,j and k_φv,j, the semi-axes a_j in meters, and the mass of one peripheral body m_j in kg. Based on these data, a multilayer structure is created. In this case, the mass of the central body is calculated based on the initial data file MLSpStr2.dat. Providing that a source data file fN3fvinp.dat is available, the MLSpStr2 program creates a structure from the data contained in this file; otherwise, algorithm (24)-(25) is used for this purpose.

The file of initial conditions for the Galactica program, for example, MS15c49b.dat, contains the masses, the coordinates, the velocities, and the radii of bodies, as well as a number of other parameters necessary for calculating the system dynamics and evolution. The Galactica system [23] [25] makes it possible to calculate the dynamics of a multilayer structure and study its evolution in time. In addition, Galactica is used to accomplish the creation of the structure. According to the algorithm presented above, a structure is created in which bodies in the layers are organized in a certain order. After their interaction for some time, the bodies will become evenly distributed over space. For implementing this distribution, the Galactica system is used.

A file of initial conditions such as MS15c49b.dat uses dimensionless values [25]. All body masses in it are related to the total mass of the system m_ss. Time T is expressed in hundreds of periods of revolution P₁ of first-layer bodies, where the periods P₁ are determined from the initial data according to formula (8). For this, the time factor

$k_t=1/\left(100\cdot P_1\right)$ (26)

is introduced. The geometric dimensions in the Galactica program are related to the quantity

$A_m={\left(G\cdot m_{ss}/k_t^2\right)}^{1/3}$ ,(27)

where G is the gravitational constant.

The program Galactica integrates the differential equations for bodies that interact according to the Newton law of gravitation. For example, in dimensionless form these equations as projected onto the x-axis look as follows:

$\frac{{\text{d}}^2{x}_j}{\text{d}{T}^2}=-{\displaystyle \sum _{k\ne j}^N\frac{m_{o,k}\left(x_j-x_k\right)}{r_{jk}^3}}$ , $j=k=\text{1},\text{2},\cdots ,N$ ,(28)

where x_j = x_C,j/A_m is the dimensionless coordinate of the j-th body; x_C,j is the dimensional coordinate of the j-th body relative to the center of mass of the entire structure, and m_o,k = m_k/m_ss is the dimensionless mass of the k-th body.

As already noted, the Galactica system, with a set of necessary tools for solving problems, is available in free access. Its description is presented in the GalDiscrp.pdf file in Russian, and in the GalDiscrpE.pdf file in English. The MLSpStr2 program, the data file MLSpStr2.dat, and the structure files mentioned here are available².

4. First Five-Layer Structures

When creating a structure, one must decide on the choice of parameters specified in the MLSpStr2.dat file. Some of these parameters were identified during the creation and study of single-layer spherical structures [20] [21]. The structures considered below have the following dimensional parameters: the initial mass of the central body and first-layer bodies is equal to the Solar-system mass mi = 1.99179 × 10³⁰ kg, with the mass fraction due to the central body being p_m0 = 0.99; the semi-axis length of first-layer bodies is equal to one astronomical unit, i.e. a₁ = 149.595 million km. In this case, the period of revolution of first-layer bodies is P₁ = 1 sidereal year. The rest parameters are as follows: A_sm = 1; e = 0; k_a = 1; k_N3 = 1; kφ = 2.83; k_φv = 1; EPS = 1 × 10⁻¹⁰; and ρ_b = 5∙10³ kg/m³. Note that the body masses are in the following correspondence to the masses of Solar-system bodies: m₀ = 0.991m_S and m₁ = 0.354m_Sa, where m_S is the mass of the Sun, and m_Sa is the mass of Saturn. The five-layer structures with these parameters, shown in Figure 3, are created with a uniform rotation of the layers around the x-axis, that is, with the key I_nx = 1.

According to formula (24), the sizes a_j of the layers in the five-layer structure MS05c99d.dat increase with each layer by the semi-axis of the orbit of first-layer bodies, and, according to formula (8), their periods P_j are equal to 2.80; 5.07; 7.66; and 10.46 periods of first-layer bodies. According to formula (25), the numbers of bodies in the layers, 99, 198, 297, 396, and 495, also increase by the number of bodies in the first layer. The total number of bodies is N = 1486. In Figure 3(a), the segments show the velocity vectors of the 2^nd and 1486^th body. Unlike in Figure 1 and Figure 2, here the enumeration of bodies begins from the central body (m₁), with the first body in the first layer designated as m₂. At the point of self-intersection of the line of formation of layers 2 and 4, the distance between bodies was respectively 66 and 14 times shorter than the average distance between bodies. Therefore, when calculating the motion of bodies in the structure using the Galactica system, intense collisions of bodies and their merging began at these places. Figure 3(b) shows the configuration of the structure at time T = 1.01, i.e., after 101 revolutions of first-layer bodies. There happened 91 collisions in this structure, while 88 bodies have acquired double masses, one a triple mass, and one body experienced a collision with the central body.

Figure 3(c) shows the distribution of relative body distances r_r1 from the center of mass depending on the body numbers k = 1, 2, …, N at time T = 1.01. The lines show the distances r_r1 at T = 0. These distances are normalized by the semi-axis of the orbits of first-layer bodies a₁, i.e. r_r1 = r/a₁. As it is seen, the greatest distance of one of the bodies reaches 20a₁. The points on the horizontal axis indicate the numbers of bodies that have merged with other bodies. It is seen that there are empty regions, or voids, in the 2^nd and 4^th layers, i.e. r_r1 = 0. In these places, there were small distances between bodies at the intersections of their location line. Therefore, at these places intense collisions of bodies occurred during their interaction.

In order to exclude such collisions, an MS05c99c.dat structure was created, in which the numbers of bodies N_3,2 and N_3,4 in layers 2 and 4 were corrected according to formula (21); as a result, the layers have become incorporating one additional body. Therefore, the total number of bodies in this structure (Figure 3(d)) has become N = 1488. In general appearance, this structure differs little from the previous one (see Figure 3(a)).

It should be noted that algorithm (21) for correcting the number of bodies N_3,j in a layer is not included in the MLSpStr2.for program. Therefore, new numbers of bodies in layers must be entered using an additional file fn3fvinp.dat.

When calculating the motion of bodies in this structure for time T = 1.05, i.e., for almost the same period as in the case of the previous structure, the number of collisions was found to equal 38. Thus, the elimination of the minimum distances in the second and fourth layers has led to a reduction in collisions by 2.4 times. During the time T = 1.95 (Figure 3(e)), there were 44 collisions in total. In this case, 40 double-mass bodies were formed, and one body had acquired a quadruple mass. In addition, there was one collision with the central body. Thus, during the second time interval ∆T = 0.9, there occurred a total of six collisions. In the latter structure (Figure 3(e)), the bodies are more uniformly distributed over space than in the previous structure shown in Figure 3(b).

Figure 3(f) shows the distribution of relative body distances for the second structure. The spread of distances in this structure is much smaller than in the previous one (Figure 3(c)). The greatest distance of one of the bodies is 13a₁. There are also no voids in the region of the 2^nd and 4^th layers.

Figure 3. Two five-layer structures at the beginning ((a), (d)) and by the end of the interaction ((b), (e)) of their constituent bodies: (a)-(c) MS05c99d.dat; (d)-(f) MS05c99c.dat; the first body of the first layer m₂ is on the x-axis; the velocity vectors of bodies m₂ and m_N are shown as segments; the lines in the graphs (c) and (f) show the body distances r_r1 from the center of mass at time T = 0.

For the two structures considered above, the problems of interaction of bodies were calculated in the Galactica system with a step dT = 1 × 10⁻⁷. Files with kinematic parameters of the structures were issued after Kl3 =1 × 10⁵ steps. This number of steps corresponded to the time interval ΔT = 0.01, which was equal to the period of revolution of first-layer bodies. For the structure MS05c99d.dat, 101 files were issued, and for the structure MS05c99c.dat, 195 files.

5. Evolution of the Five-Layer Structure

5.1. General Changes

In the structures shown in Figure 3, the layers were successively rotated about the x-axis; as a result, these structures became asymmetric. In this connection, the MS05с99e.dat structure was created with the same parameters as the MS05с99с.dat structure, but without the rotation of layers, that is, with the parameter I_nx = 0. As evident from Figure 4(a), the former structure is more symmetrical. It has N = 1488 bodies, and the layer parameters in it are the same as in the MS05c99c.dat structure. After 196 revolutions of first-layer bodies, that is, at T = 1.96, the structure under consideration (see Figure 4(b)) also proved to be more symmetrical than the previous configuration (Figure 3(e)). However, the spread of body distances r_r1 in it (Figure 4(c)) is somewhat greater than in the previous structure (Figure 3(f)). The greatest distance of a body from the center of mass is r_r1max = 15.7 in comparison with r_r1max = 13.1 in the previous structure. It also exhibits slightly more collisions, namely, 46, compared with the previous structure with 44 collisions. Note that body 1434 exhibits the greatest distance from the center of mass r_r1max = 15.7.

The dynamics of inter-body collisions in the structure under consideration is shown in Figure 4(d). By the time T = 1.96, 40 bodies with a double mass and 3 bodies with a triple mass have formed in this structure. As it is seen from Figure 4(d), collisions occur more frequently during the initial time interval of T < 0.5, and less frequently during the interval T > 0.7. Initially, the rate of collisions is equal to v_imp = 64 collisions per 100 revolutions of first-layer bodies and, then, v_imp = 8. That is, the rate of collisions has decreased by eight times. If we relate the rate of collisions to the number of bodies, then in the last section it will be v_imp1 = v_imp/N = 5.4 × 10⁻³.

Figure 4. Layer structure MS05c99e.dat without layer rotations. (b) at time T = 1.96. Part (d) of the figure shows the dynamics of collisions: k_0m is the number of collisions; k_2m is the number of formed bodies with a mass of 2m₁; k_3m is the number of formed bodies with a mass of 3m₁. For the rest designations see Figure 3.

From the visual analysis of the structures, it follows that the remains of their initial organization shown in Figure 4(a) still persist by the time T = 0.1, whereas by the time T = 0.2 they are no longer observed. That is why we can assume that, from the latter moment, a structure with a disordered arrangement of bodies has been created, with its appearance changing little in the future. When the rate v_imp reaches the second stage of its change, we can assume that the dynamics of the structure has passed into the stable phase of its existence.

5.2. Trajectories of Individual Bodies

The trajectories of bodies in different layers were studied: that of first-layer body 2, that of third-layer body 401, and those of bodies 1434 and 1488 in the last layer. The study was carried out for two time intervals, initial T = 0 - 0.15 (Figure 5(a), Figure 5(b)) and final T = 1.81 - 2.11 (Figures 5(c)-(e)).

Over the initial time interval, body 2 first moves in a circle (see Figure 5(a)) with a period P₁ = 0.01. Its orbit lies in the xy-plane. Over time, the orbit becomes elliptical, and the period of revolution increases slightly. Over a finite time interval (Figure 3(c)), the orbital eccentricity of body 2 increases to e = 0.385, but the dimensions of the orbit show a decrease, and the period also decreases and becomes shorter than 0.01. In the 3D graph of Figure 5(e), it is seen that the orbit of body 2 rotates in space, and with each revolution of this body

Figure 5. Trajectories of bodies 2, 401, 1434, and 1488 over the initial ((a), (b)) and final (c)-(e) periods of evolution of the MS05c99e.dat structure: I and F are the initial and final points of the trajectories; (e) three-dimensional images of the trajectories; the circle in image (d) shows the position of body 1434 at time T = 1.96.

around the center of mass, this body departs farther and farther from the xy-plane.

The orbit of body 401 over the initial time interval (Figure 5(a)) during the first revolution of the body is approximately a circle with a period of P₂ = 0.028. Then, the trajectory becomes elliptical, and the period of revolution increases slightly. Over a finite time interval (Figure 5(c)), the eccentricity reaches a value e = 0.373, the orbit size increases, and the period reaches 1.25P₂. As it is seen from Figure 5(e), here the orbit rotates in space with each body revolution.

The orbit of body 1488, the last body in the fifth layer, over the initial time interval (Figure 5(b)) is a circle with a period of revolution P₅ = 0.1046. Over the final time interval (Figure 5(d)), the orbit is already elliptical with an eccentricity e = 0.349, and the period has increased by a factor of 1.6. From Figure 5(e), it is evident that the orbit of body 1488 changes in space little during this time.

Figure 5(b) shows that during the first revolution the orbit of the most distant body 1434 exhibits notable changes. Over the final time interval (Figure 5(d)), the orbit is already an ellipse with an eccentricity e = 0.491 and a period of 0.30, while the initial period of revolution of five-layer bodies is 0.1046. That is, here the period has increased three-fold.

It should be noted that the orbits are exactly ellipses, parabolas and hyperbolas only in the two cases: 1) the interaction involves two bodies, and 2) the interaction involves N bodies located axisymmetrically on one plane [18]. That is why under an elliptic trajectory we mean a trajectory bounded in space, while under hyperbolic and parabolic trajectories we understand infinite trajectories. In the latter case, the velocity of a body on a parabolic trajectory tends to zero at infinity.

5.3. Dynamics of Encounters of an Individual Body with Other Bodies

Changes in the orbit of body 1434 occur when this body approaches other bodies. The approaches of body 1434 over the initial time interval are shown in Figure 6. Plotted on the vertical axis are the numbers of bodies k to which the body i = 1434 comes in proximity. Such bodies are marked with points or with

Figure 6. Approach of body i = 1434 by bodies k to a distance r_i,k for the number of integration steps Kl3 = 3000 made with a time step of dT = 1 × 10⁻⁷ over the initial time interval T = 0 - 0.15. The dimensionless distance r_i,k in the graph is increased by f = 3 × 10⁴ times.

horizontal segments if the points overlap. Also indicated on the vertical axis is the distance R_ik at the moment of approach of body 1434 to body k. The distance R_ik is increased by the factor f. The distance R_ik is defined as the minimum distance of a body i from other bodies over the number of integration steps Kl3, where Kl3 is a parameter of the Galactica system. In the case under consideration, the integration was performed with a step dT = 1 × 10⁻⁷.

As it follows from Figure 6, at the initial time T = 0 body 1434 passes at a distance of 2.40 × 10⁻³ from body 1062, which is at the intersection of the initial location line. In this case, the distance between the bodies on this line is 3.17 × 10⁻³, that is, of the same order.

Then, body 1434 moves with an almost constant distance to the neighboring body 1435 until the time T = 2.16 × 10⁻²; and then at time T = 2.32 × 10⁻² it approaches body 1311 to a distance R_ik = 1.66 × 10⁻⁴. The latter distance corresponds to 90 radii of the body. Then, at time T = 3.44 × 10⁻² the body approaches body 1186 to a distance of R_ik = 5.84 × 10⁻⁴. These two approaches lead to a significant change in the trajectory of the body of interest, with its orbit becoming elliptical. The body moves away from all bodies to a distance R_ik = 3.83 × 10⁻² from the nearest body 1421 at time T = 0.12. This is the most distant point of its orbit, after which body 1434 starts approaching the center of the structure again. At this point, the distance of body 1434 from the center is r = 0.109.

As a result of subsequent interactions, the orbital eccentricity increases in magnitude, and at the apocenter by the time T = 1.96 (Figure 5(d)) the body moves away from the center to a distance r = 0.205, which value is 3.14 times greater than the initial size of the structure.

5.4. Determining the Trajectory Parameters of a Distant Body

When considering the results of calculations for the motion of bodies in the structures shown in Figure 3(b), Figure 3(d), and Figure 4(b), it becomes necessary to determine the type of motion of remote bodies. This allows one to determine whether such a body is a body of this structure or it is a body ejected out of it. For solving this problem, it is necessary to perform a special study of the trajectory of such body. Consider a method based on the results of the two-body problem, which will give an answer to the question of interest without performing special studies.

For a body remote from the structure, we can assume that it is affected by the entire structure with mass m_ss located in its center of mass. Then, similarly to (5), the parameter of the body trajectory will be

${\alpha }_1={\mu }_1/\left(R_p\cdot v_p\right)$ (29)

where R_p and v_p are the pericenter radius of the body and its velocity in pericenter, and the interaction parameter is ${\mu }_1=-G{m}_{ss}$ . As a result of calculations made using the Galactica program, we obtain the coordinates x, y, z and the velocities v_x, v_y, v_z, of the body in the center-of-mass system. From the definition of the scalar product of the radius vector r of the body and its velocity vector v, we can write

$x{v}_x+y{v}_y+z{v}_z=r\cdot v\cdot \cos {\beta }_1$ ,(30)

where β₁ is the angle between the vectors r and v;

$r=\sqrt{x^2+y^2+z^2}$ ; $v=\sqrt{v_x^2+v_y^2+v_z^2}$ .

Then, from formula (30) we obtain the following expression for the angle β₁:

$\cos {\beta }_1=\frac{x{v}_x+y{v}_y+z{v}_z}{r\cdot v}$ (31)

The radial velocity v_r is directed along the radius vector r, and the transversal velocity v_t, in the direction perpendicular to the latter velocity, so that

$v_r=v\cdot \cos {\beta }_1$ ; $v_t=v\cdot \sin {\beta }_1$ .(32)

On the other hand, according to the two-body problem [33] [34], the radial and transversal velocities, similarly to (12) and (13), can be written as follows:

$v_r=\pm v_p\sqrt{{\left({\alpha }_1+1\right)}^2-{\left({\alpha }_1+R_p/r\right)}^2}$ ; $v_t=v_p\cdot R_p/r$ .(33)

Equating the right-hand sides in the velocity v_t from formulas (32) and (33), we obtain the pericentric velocity

$v_p=\frac{v\cdot r\cdot \sin {\beta }_1}{R_p}$ ,(34)

and excluding the velocity v_r from formulas (32) and (33) with taking into account formula (34), we obtain the expression

$R_p\cos {\beta }_1=r\sin {\beta }_1\sqrt{{\left({\alpha }_1+1\right)}^2-{\left({\alpha }_1+R_p/r\right)}^2}$ .(35)

Three Equations (29), (34) and (35) include three parameters α₁, R_p, and v_p. As a result of successive substitutions and solutions of quadratic equations, the pericenter radius is obtained in the following form:

$R_p=\frac{{\mu }_v\pm \sqrt{{\mu }_v^2+r^2\sin {\beta }_1\left(1+2{\mu }_v/r\right)}}{1+2{\mu }_v/r}$ (36)

where the designation μ_v = μ₁/v² is introduced. The quantity μ_v, which is measured in meters, is negative, μ_v < 0. Expression (36) gives two values of R_p: with the “-” sign in the case of an elliptical orbit we obtain the apocenter radius R_a, and in the case of “+”, the pericenter radius R_p.

Given the pericenter radius R_p, the pericenter velocity v_p can be found using expression (34). In the case of v_p > 0, the orbit is passed counterclockwise. Then, formula (29) can be used to calculate the trajectory parameter α₁, and formula (10) yields the orbital eccentricity e. In accordance with formula (4), these parameters make it possible to determine the trajectory of the body, as well as the time of motion along it [33] [34].

In the case of a parabolic or hyperbolic orbit, the velocity at infinity can be calculated as follows:

$v_{\infty }=v_p\sqrt{2{\alpha }_1+1}$ (37)

For a hyperbolic orbit, the apocenter radius R_a calculated using expression (36) with the “-” sign turns out to be negative.

For body 1434, the parameters calculated by this algorithm for time T = 1.96 have the following values: R_p = 0.0694, v_p = 4.638, α₁ = –0.669, e = 0.494, R_a = 0.2049, and v_a = 1.572. Since for an elliptical orbit we have: –0.5 > α₁ > –1, the trajectory of body 1434 is an ellipse. At time T = 1.96, the body is at a distance r = 0.2049 and has a velocity v = 1.572. These parameters coincide with the parameters of the apocenter. As it is seen from Figure 5(d), body 1434 shown with a circle is indeed located at the most distant point of the trajectory from the center, that is, at its apocenter. Therefore, algorithm (29)-(37) presented above can be used to estimate the trajectory of a remote body from its coordinates and velocity at some point in time.

6. Evolution of the Ten-Layer Structure

6.1. General Changes

Due to the fact that in the structure with layers rotated around the x-axis (Figure 3(d)) there were slightly fewer collisions than in the structure without rotations (Figure 4(a)), we have decided to create 10- and 15-layer structures with rotated layers. Such a ten-layer structure is shown in Figure 7(a). The initial parameters of this structure in the MLSpStr2.dat file are the same as those of the structure in Figure 3(d), except for the number of layers, which is set to N₂ = 10. All close encounters at the self-crossings of the formation line are eliminated. The number of bodies in the layers varies from 99 in the 1st layer to 991 in the 10^th layer, the semi-axis varies from a₁ = 0.0118 to a₁₀ = 0.1179, and the period, from P₁ = 0.0100 to P₁₀ = 0.2536. The total number of bodies is N = 5451.

In the previous cases, the calculation of bodies’ motion implemented with the help of the Galactica program was carried out with a time step dT = 1 × 10⁻⁷. At this step, the calculation time for the interval ∆T = 0.01, which is equivalent to one revolution of first-layer bodies lasted for 7 hours. For a ten-layer structure with N = 5451, the calculation time of this interval took 97 hours. For solving the problem over the required time interval, the duration of calculation will exceed one year. Therefore, the solution of problems with such a number of bodies was performed with a step of dT = 10⁻⁶ in automatic step selection mode, which can be launched in the Galactica system using the key Kl4 = 3. Since the results of calculation are issued following a certain number of integration steps, in the latter case the time interval between these results may be different.

It should be noted that the Galactica system also automatically modifies the step when the bodies approach each other up to a distance of the order of their diameters. These circumstances must be taken into account when analyzing points on the graphs of parameter variation over time.

After 81.5 revolutions of first-layer bodies, i.e. at T = 0.815, the structure is shown in Figure 7(b). By the time T = 0.1, elements of the initial organization of

Figure 7. The ten-layer structure MS10c99b.dat and its evolution. (b) data at T = 0.815. Part (d) of the figure shows the dynamics of collisions: k_0m is the number of collisions; k_4m is the number of formed bodies with a mass equal to 4m₁; k_5m is the number of formed bodies with a mass equal to 5m₁. For other designations, see Figure 3 and Figure 4.

the structure are still preserved, and by the time T = 0.2 they completely disappear. By the time T = 0.4, the structure acquires a form later showing almost no changes. The scatter of body distances in Figure 7(c) is limited to 40a₁. About 30 bodies were ejected over a large distance, the largest of which is r_r1 = 600 for body 1575.

By the time T = 0.815, there were 306 collisions, one of which was with the central body. During this period, 202 double-mass bodies, 34 triple-mass bodies, 9 quadruple-mass bodies, and 2 bodies with mass 5m₁ were formed. The dynamics of collisions is shown in Figure 7(d). Before the time T = 0.061, there were no collisions. Then frequent collisions began, with the rate which we will also consider in two regions. On the first interval ∆T = 0.0688 - 0.142 the rate of collisions is v_imp = 2158 collisions per 100 revolutions of first-layer bodies, and on the second interval ∆T = 0.339 - 0.815 the average rate of collisions is v_imp = 69.3. At the same time, over this time interval there is a tendency towards a further decrease of velocity, possibly to zero at infinity. In this structure, the rate of collisions in the second trajectory section compared with the first section has decreased by 31 times. In this case, the rate of collisions per one body is v_imp1 = 1.27 × 10⁻², which value is 2.4 times higher than that in the five-layer structure.

6.2. Trajectories of Individual Bodies

Figure 8 shows the trajectories of the central body 1, first-layer body 57, second-layer body 170, sixth-layer body 1563, and tenth-layer body 5301. These calculations were performed over the interval T = 0.663 - 1.021 with a step dT = 1 × 10⁻⁷. The trajectory of the central body 1 around the center of mass is irregular (Figure 8(a)). The three-dimensional appearance of the trajectory is shown in Figure 8(d). Body 1 does not move away from the center of mass to a distance greater than 0.003, which value is 0.23 of the radius of the inner layer.

The orbit of body 57 (Figure 8(b)) continuously rotates in space (Figure 8(e)). The period of revolution of body 57 fluctuates within small limits relative to the period 1.34 P₁, where P₁ is the initial period. The semi-axis of the orbit also oscillates around a value of 0.015. The orbital eccentricity reaches 0.7.

Figure 8. Trajectories of the central body 1 and peripheral bodies 57, 170, 1563, and 5301 during the final period (T = 0.663 - 1.021) of the evolution of the MS10c99b.dat structure (a)-(c): parts (d) and (e) of the figure show 3D images of the trajectories; and I and F are the initial and final points of the trajectories.

The orbit of body 170 is an ellipse with a large eccentricity reaching 0.9 (Figure 8(b)). This orbit also varies in space (Figure 8(e)). The orbital period fluctuates around 1.16P₂, where P₂ = 0.028 is the initial period.

It should be noted that at the pericenters of the orbits in Figure 8(b) and Figure 8(d) there are trajectory sections with straight segments. This is due to the fact that during the intervals of outputting the trajectory coordinates (with the Galactica system parameter Kli = 2000), the orbits change more substantially than in other regions.

The trajectory of body 1563 is almost a circle, its eccentricity being e = 0.01 (Figure 8(c)). The orbit slightly changes in space (Figure 8(e)). The orbital period is 1.2P₆, where P₆ = 0.134. The semi-axis of the orbit has also changed little compared to the initial one and is equal to 1.06a₆.

The trajectory of body 5301 is an ellipse with eccentricity e = 0.23 (Figure 8(c)). Its plane is almost perpendicular to the xy-plane (Figure 8(e)). That is why on the xy-plane this trajectory is depicted as a flattened ellipse (Figure 8(c)).

As already noted, the distances r_r1 of about 30 bodies were outside the graph in Figure 7(c). For some of these bodies, according to algorithm (29)-(36), the types of trajectories and their parameters were determined. At time T = 0.815, the most remote body 1575 had a distance r = 7.076 and a velocity v = 9.735. Calculations by algorithm (29)-(36) gave the following results: R_p = 0.0560, v_p = 11.41; α₁ = −0.1371; e = 6.293; and v_∞ = 9.720. With this value of α₁, the trajectory is a hyperbola. Thus, the body was ejected from the structure along a hyperbolic orbit. At this time, its velocity already approached the velocity at infinity.

For body 1575, the calculations were repeated for the time T = 0.385, when the distance and velocity were r = 2.883 and v = 9.766. The parameters of the hyperbolic trajectory were confirmed up to seven digits. These calculations were also performed for the time T = 0.1035, when the parameters of the body were r = 0.1907 and v = 10.450. The parameters of the hyperbolic orbit were confirmed to within 4 digits. At this time, the distance r is close to the radius of the outer layer. It follows from here that algorithm (29)-(36) can be used for determining the trajectories of all bodies outside their main cluster.

These studies were carried out for all remote bodies. It was found that 14 bodies with mass m₁ and 3 bodies with mass 2m₁ were ejected from the structure. Two bodies have strongly elongated elliptical orbits, whose apocenter radii R_a are equal to 0.5646 for body 4601 and 3.024 for body 3032.

6.3. Approaches, Collisions, and Merging of Bodies

As already noted, in this structure there are two bodies with a mass of 5m₁. Figure 9 shows the encounters of one of these bodies, namely body 1951, with other bodies over the interval 0.069 < T < 0.075. The calculations were performed in step correction mode with key Kl4 = 3 and key Kl3 = 300 steps for outputting results. From the time T = 0.069, body 1951 approaches body 1952, and at time T = 7.153 × 10⁻² the bodies merge together. Then, the approach to body 1950 begins and at time T = 7.153 × 10⁻² the two bodies merge together.

Then, body 1951 approaches body 1953 and at time T = 7.172 × 10⁻² to merge with it. Then, the approach to body 1954 begins, and at time T = 7.321 × 10⁻² there occurs merging with this body. As a result of these four events, the mass of 1951 became equal to 5m₁. Further, at time T = 7.321 × 10⁻² the body 1951 approaches body 1949 to a distance of 14 peripheral-body radii R₁ = 1.676 × 10⁻⁶, and the bodies diverge. Then, at time T = 7.456 × 10⁻² there occurs an encounter with body 1955 to a distance of 48R₁.

Figure 10(a) shows the trajectories of all bodies approaching and colliding with body 1951. These calculations were performed in mode Kl4 = 3 with an initial step dT = 0.5 × 10⁻⁶ (Figure 10(b)) over the interval of output steps Kli =

Figure 9. Approach of body i = 1951 by bodies k to a distance r_i,k for the number of integration steps Kl3 = 300 made with a step dT = 0.5 × 10⁻⁶ over the initial interval T = 0.0690 - 0.0750. The dimensionless distance r_i,k in the graph is increased by a factor of f = 5 × 10⁴, and the number of bodies k is reduced by 1940.

Figure 10. The trajectories during the approach and collision of body 1951 with bodies 1949, 1950, 1952, 1953, and 1954. Table 1 shows the values of y₁ − y₁₂, x₁ − x₁₂.

Table 1. The values of the bounds in the graphs b, c, d, e, f, and g of Figure 10.

b		c		d		e		f		g
y1	0.06106	y3	0.061054	y5	0.06236	y7	0.06274	y9	0.06483	y11	0.0650
y2	0.06104	y4	0.061046	y6	0.06232	y8	0.06270	y10	0.066480	y12	0.0645
x1	−0.03948	x3	−0.039470	x5	−0.03762	x7	−0.03715	x9	−0.03372	x11	−0.0338
x2	−0.03946	x4	−0.039462	x6	−0.03758	x8	−0.03713	x10	−0.03370	x12	−0.0333

300. The first approach of body 1951 occurs with body 1952. On the section of approach (Figure 10(b)), the calculations were carried out with the interval of output steps Kli = 3. Body 1952 moves ahead of body 1951 and, as the latter body approaches the former, it begins to go around it. On a larger scale (Figure 10(c)), the diameters of the bodies are marked with circles. The bodies are almost at the same distance from the xy-plane. The merging of the two bodies occurs when their surfaces come in contact. The condensation of points on the trajectory indicates the inclusion of step correction mode before the collision. The length of the condensation section is proportional to the body velocity. It is seen that the velocity of body 1951 is greater than the velocity of body 1952. This difference may be due to the interaction between these bodies. In this case, the body 1951 is accelerating, and the body 1952, decelerating.

It should be noted that all images in the graphs are drawn to scale, and the values of the limits y₁, y₂, x₁, x₁, etc. are presented in Table 1.

The second approach of body 1951 occurs with body 1950 (Figure 10(a)). Figure 10(d) shows that the body 1950 catches up with the body 1951, crosses its trajectory, and due to the attraction of this body, the trajectory of body 1950 bends. The merging of the two bodies occurs at the moment of contact. In this case, the z-coordinate of body 1951 is greater, so the contact occurs below the visible contour of body 1951. As a result of the merging with body 1952, the radius of body 1951 increases to 1.26R₁. Before the collision, the body 1950 was slowing down as a result of the interaction with body 1951, while the body 1951 experienced acceleration.

The third approach of body 1951 occurs with body 1953 (Figure 10(a)). It is seen from Figure 10(e) that body 1953 moves ahead of body 1951. That is why body 1953 goes around body 1951 and collides with it on the side opposite to the side of approach. In this case, the radius of body 1951 was already 1.44R₁.

As it is evident from Figure 10(f), the fourth approach to body 1954 occurs similarly, yet with the enveloping motion of body 1954 being more pronounced than that of body 1953. This is due to the greater mass of body 1951. This mass is 4m₁, and the radius of 1951 is 1.59R₁. After the merging of body 1951 with body 1954, its mass became 5m₁, and its radius, 1.71R₁. In these two cases (see Figure 10(e), and Figure 10(f),), body 1951 more closely approaches the xy-plane than body 1954; therefore, this plane is partially shaded by this body.

In the four collisions considered above, three bodies approached the target body from the left side and collided with it on the right side (Figure 10(c), Figure 10(e), Figure 10(f)). On the other hand, the body approaching the target body from the right side collided with it on the left side (Figure 10(d)). At the same time, the greater the mass of body 1951, the greater is the angle through which the approaching body envelopes body 1951.

The approach of body 1949 to body 1951 is shown in Figure 10(g). These calculations were performed with the output interval Kli = 30. As noted above, in this case the body 1949 approached body 1951 to a distance of 14R₁. At the moment of approach, the velocity of body 1949 was almost three times greater than the velocity of body 1951. That is why the body 1949 overtook body 1951 with a further increase in velocity and went around it by almost 90˚. After approach, the velocity of body 1949 began to decrease.

Despite the above collisions, body 1951 kept moving within the structure. By the end of the studied time interval T = 0.815, the distance of body 1951 from the center of mass was r_r1 = 7.17. Body 1949 was also inside the structure with a distance of r_r1 = 4.28.

6.4. Spin of the Body and Its Thermal Energy after Collision

When bodies collide, the body formed from them acquires its own angular momentum, which we call spin S_p, and thermal energy E_t. The algorithm for calculating these quantities was reported in [33] [34], and its software implementation was described in [24]. Figure 11(a) shows the variation of spin projections S_px, S_py, S_pz and thermal energy E_t of body 1951 during its collision with four bodies. As a result of each collision, these quantities changed. The circles in the graphs mark the moments T of issuing output files by the Galactica program; those files were used to determine the values of interest.

The smallest change in spins occurs during the first collision (point 2m₁) with body 1952 (Figure 10(c)). This collision is close to a frontal impact. As is seen from Figure 11(a), the thermal energy E_t changes significantly in this case. During the second collision (point 3m₁) with body 1950, the spin projection S_p has increased more significantly, and the thermal energy E_t changed three times less than it did so during the first collision. As it is seen from Figure 10(d), this collision occurred along a tangent line. Other collisions with bodies 1953 and 1954 also occurred in nearly tangential directions. But body 1951 became more

Figure 11. Variation of thermal energy E_t and spin projections S_px, S_py, S_pz formed during the merging of bodies 1951 (a) and 3685 (b). The value of E_t is multiplied by the coefficients f = 2 × 10⁻⁶ (a) and f = 2 × 10⁻⁷ (b).

massive, and the velocity of the approaching bodies increased. Therefore, the increment of spins and thermal energy were significant. As a result, body 1951 acquired thermal energy E_t = 2.87 × 10⁻³ and spin modulus S_p = 6.385 × 10⁻⁹. The spin vector makes an angle β₂ = −23.8˚ with the z-axis. Since the angle is negative, body 1951 rotates clockwise.

For comparison, Figure 11(b) shows the rotational and thermal characteristics of body 3685, which also has a mass of 5m₁. During the first collision at point 2m₁, the body has acquired a small thermal energy E_t. During the second collision (3m₁), the spin projections S_px and S_py have slightly decreased, while the thermal energy E_t increased significantly. The third collision at point 5m₁ occurred with a body of double mass. In this case, the spin projection S_py dominates. The spin value is S_p = 3.48 × 10⁻⁹, and its vector makes an angle β₂ = −78.1˚ with the z-axis, i.e. the axis of rotation of body 3685 is close to the xy-plane. Its thermal energy is E_t = 6.45 × 10⁻³. Thus, with the same masses of these bodies, body 1951 has a 1.8 times greater spin and 2.2 times lower thermal energy.

7. Evolution of the 15-Layer Structure

7.1. General Changes

Due to the fact that with an increase in the number of bodies, the time required for calculating the evolution of a structure of interest becomes excessively large, a 15-layer structure was created with a smaller initial number of bodies in the first layer, namely N_3,1 = 49. In this case, the absolute peripheral-body mass was 2.02 times greater than that in the previous structures. This structure, contained in the MS15c49b.dat file, is shown in Figure 12(a). All close encounters at self-crossings of the formation line were eliminated by changing the number of bodies in accord with (21). The number of bodies in the layers varies from 47 in the first layer to 733 in the 15^th layer; the semi-axes, from a₁ = 0.0118 - 0.0105 to a₁₅ = 0.1573; and the periods, from P₁ = 0.0100 to P₁₀ = 0.3916. The total number of bodies was N = 5866.

After 77.9 revolutions of first-layer bodies, i.e. at T = 0.779, the structure is shown in Figure 12(b). By the time T = 0.2, elements of the initial organization of the structure are still preserved, whereas by the time T = 0.3 they completely disappear. By the time T = 0.4, the structure acquires a form that shows almost no subsequent changes. The scatter of body distances in Figure 12(c) is limited to 50 radii of first-layer bodies. About 20 bodies were ejected over a greater distance, the largest of which is r_r1 = 168 for body 5336.

By the time T = 0.779, there occurred 164 collisions. In this case, 120 bodies of double mass, 17 bodies of triple mass and 2 bodies of quadruple mass were formed. In addition, there occurred four collisions with the central body.

The dynamics of collisions is illustrated in Figure 12(d). For the interval T = 0.085 - 0.286, the rate of collisions was v_imp = 547 collisions per 100 revolutions of first-layer bodies, and during the second interval T = 0.378 - 0.779, the average rate of collisions was v_imp = 85. Thus, over the second time interval the collision rate has decreased by 6.4 times.

Figure 12. The fifteen-layer structure MS15c49b.dat and its evolution. Image b shows the structure at time T = 0.779. For the rest designations, see Figure 7.

When compared with the 10-layer structure, the velocity over the second time interval of the 15-layer structure was 1.23 times higher. The velocity per one body, v_imp1 = 1.45 × 10⁻², is also 1.14 times higher.

7.2. Trajectories of Individual Bodies

Figure 13 shows the trajectories of the central body 1, first-layer body 25, seventh-layer body 1325, and 15^th-layer body 5660. These calculations were carried out over the interval T = 0.597 - 0.853 with a step dT = 1 × 10⁻⁷. The trajectory of the central body 1 (Figure 13(a)) is irregular and lies in the range of distances r < 0.003 from the center of mass.

The orbit of body 25 is continuously changing in space. Its eccentricity also changes: it decreases up to the time T = 0.77 and then increases. The average value of the eccentricity is e = 0.69. The period undergoes similar changes, and on average it equals 0.83P₁, where P₁ = 0.01 is the initial period. After T = 0.77, a significant change in the plane of the orbit occurs (Figure 13(c)), both its eccentricity and period increase, and its change in space also becomes more pronounced.

The orbit of body 1325 (Figure 13(b)) is an ellipse with eccentricity e = 0.408. The semi-axis of the orbit is 1.1 times greater than a₇ = 0.073, where a₇ is the semi-axis of the initial orbit, and the period is 0.95 times less than P₇ = 0.1645. As it follows from Figure 13(b) and Figure 13(d), the orbit changes significantly during one revolution.

The orbit of body 5600 in Figure 13(b) and Figure 13(e) is represented by its part smaller than half the body’s range of revolution about the center of mass. The distance of the body from the center is r = 0.238 and practically does not change. The slight ellipticity in Figure 13(b) is due to the inclination of the orbit to the xy-plane (see Figure 13(d)). The distance r is 1.5 a₁₅, and the period of revolution of body 5600 is therefore longer than the initial period P₁₅ = 0.3916.

As already noted, the distances r_r1 of approximately 20 bodies were beyond the bounds of Figure 12(c). For most of these bodies, algorithm (29)-(36) was used to determine their trajectory types. For bodies 5336, 543, 264, and 5146, remote to distances r_r1 = 168, 125, 115, and 107, respectively, the trajectories were hyperbolas, i.e. the bodies were ejected out of the structure. In total, six bodies were ejected. In this case, for two bodies, 5146 and 264, the eccentricities were 1.045 and 1.0098, respectively, i.e. they are close to unity, which value characterizes a parabolic orbit. For body 5594, remote to r_r1 = 94, the orbit is an ellipse with eccentricity e = 0.894. Other distant bodies also have elliptical orbits.

Figure 13. Trajectories of the central body 1 and peripheral bodies 25, 1325, and 5600 during the final period (T = 0.597 - 0.853) of the evolution of the MS15c99b.dat structure ((a), (b)): shown in parts (c)-(e) of the figure are the three-dimensional images of the trajectories; and I and F are the initial and final points of the trajectories.

7.3. Approach of a Remote Body

Encounters of the most distant body 5336 with other bodies were studied (Figure 14). In the first time interval T = 0.184 - 0.193 (see part a of Figure 14), this body first passes close to body 5335, and then its approach to body 5333 begins. At time T = 0.1854, the distance between the bodies R_ik becomes equal to 4.05 × 10⁻⁶, this value amounting to 2.15 radii of the body R₁, i.e. the distance between the surfaces of the bodies is 0.15∙R₁. In this case, the velocities of the bodies increase to v = 5.53 for body 5336 and v = 5.37 for body 5333. At the initial time T = 0, the velocities of both bodies were equal to 2.53. After the approach, the bodies move away from each other and their velocities at the final time T = 0.779 decrease to v = 2.77 and 1.87, respectively.

In the final time interval T = 0.75 - 0.85, it is seen that body 5336 is continuously moving away from body 5146, which is also moving away from the structure along a hyperbolic orbit. The distance between the bodies increases almost linearly with an average velocity of 1.38. The bodies move in approximately the same direction: the angle between their velocities is 8.27º.

As for body 5333, after moving away from body 5336, its velocity, as already noted, decreased more significantly, and its motion proceeded along an elliptical orbit with an eccentricity e = 0.597. At time T = 0.779, the body is at a distance of r_r1 = 13.87. Thus, the ejection of body 5336 has occurred due to its acceleration when approaching body 5333 to a distance r = 2.15R₁.

In the considered 5, 10, and 15-layer models of globular clusters, the arrangement of the layers relative to each other, the number of bodies in the first layer, and the number of layers changed. All of them are stable and do not destroy. Therefore, with the variations considered, it is possible to create models of globular clusters with any number of layers in them.

Figure 14. Graph of the approaches of body i = 5336 to bodies k to a distance r_i,k for the number of integration steps Kl3 = 400 with a variable step (Kl4 = 3) over the first interval T = 0.184 - 0.193 (a) and over the second interval T = 0.75 - 0.85 (b). The dimensionless distance r_i,k in the graph is increased by a factor of f = 1·10⁴, and the number of bodies k in part a of the figure is reduced by 5300.

8. General Characteristics of Structures

8.1. Velocity Profiles

In Figure 15(a) and Figure 15(b), the points show the velocity profiles of the five-layer MS05c99e.dat structure at the final time T = 1.96, and the lines show the same profiles at T = 0. In Figure 15(a), the velocities v are shown depending on the body numbers k, and in Figure 15(b), depending on the relative radius r_rN2 = r/r_N2. At the initial time T = 0, the line segments mark the velocities of bodies in five layers (Figure 15(a)). They decrease stepwise from the first layer to the last fifth layer. After the interaction of bodies proceeding during 196 revolutions of first-layer bodies, the distribution of velocities no longer shows the initial stepping behavior. The velocities were distributed around the initial velocity profile. At the same time, most of the velocities retain the initial trend of decreasing the velocity magnitude with increasing the body number k.

The zero velocities on the horizontal axis mark the bodies that have merged with other bodies, and their velocities were therefore assigned zero values. Body 76, which is at a distance r_rN2 = 0.082, or r_r1 = 0.41 from the center of mass, has the highest velocity v = 14.33; i.e. it is located near the central body. All in all, there are eight bodies whose velocities exceed the limit of the graph v = 12.

The distribution of velocities over body distances r_rN2 (Figure 15(b)) has a more compact form. In the trajectory section r_rN2 ≤ 1, the velocities of bodies are distributed uniformly around the initial profile represented by segments. It should be noted that the initial velocity profile is discrete. It is represented by points at which the line breaks.

Consider the reasons for the deviation of velocities from the initial profile. At first the orbits are circles, and then they become ellipses. In the pericenters of the orbits, the bodies have high velocities, and in the apocenters they have lower velocities. In addition, for some bodies the semi-axes of their orbits become smaller, and the velocities become greater. For the other part of the bodies, on the contrary, the semi-axes increase, and the velocities decrease. For these two reasons, there appears a velocity spread relative to the initial profile. Evidently, the action of these two factors is symmetrical with respect to this profile.

Figure 15. Distributions of body velocities v over their numbers k and over the relative body distances from the center of mass r_rN2 for structures formed by 5 ((a), (b)), 10 ((c), (d)), and 15 ((e), (f)) layers.

According to formulas (12)-(13), the velocity of the bodies at the initial time is equal to

$v=\sqrt{v_r^2+v_t^2}=v_p\sqrt{2{\alpha }_1\left(1-R_p/r\right)+1}$ , (38)

where the values of α₁, v_p, and R_p can be calculated using dependencies (4)-(9). For circular orbits with r = R_p, formula (38) gives the velocities of bodies at the initial time, and for arbitrary distances r it gives such velocities for all bodies. The above parameters depend on the mass of the bodies within the radius r, as well as on the number of bodies and their mass at the radius r (see formulas (6) and (7)). Therefore, expression (38) can be used to estimate the masses of stars in globular clusters from their velocities.

Outside the distances of the graph in Figure 15(b), there is only one body 1434 with r_rN2 = 3.14. The velocity of this body, v = 1.572, falls onto the same dependence v(r_rN2) as for other bodies at distances r_rN2 close to 3.

The velocity distributions of the ten-layer structure over body numbers (Figure 15(c)) are more compact than those for the five-layer structure. The velocities are distributed more uniformly about the initial profile for layers 1 to 8. The distribution of velocities over distances (Figure 15(d)) is also more compact compared to the distribution over body numbers. The highest velocity v = 17.7 is exhibited by the body located at a distance of r_rN2 = 0.0227, or r_r1 = 0.227, so that this body is almost twice as close to the center as in the five-layer structure. A total of four bodies exceed the velocity bounds on the graph. The bounds of the graph over distances are exceeded for 30 bodies, with the largest distance r_rN2 = 60 belonging to body 1575. The velocity v = 9.735 of this body does not fall on the dependence v(r_rN2) for distances close to 3. Starting from r_rN2 = 6, the bodies have an ascending velocity profile. This profile can be represented as an average dependence

$v=-0.526+0.171\cdot r_{rN2}$ for $r_{rN2}>6$ .(39)

These bodies leave the structure by moving along hyperbolic orbits.

The 15-layer structure has an even more compact distribution of bodies over their numbers (Figure 15(e)). It lasts up to the 14^th layer. The velocities here are distributed more evenly around the initial velocity profile.

The distribution of velocities over distances is also more compact (Figure 15(f)). However, the middle of this distribution noticeably deviates from the initial velocity profile in the region of the last layers. This difference is observed already for the 10-layer structure in Figure 15(d). In the 15-layer structure, the highest velocity v = 16.9 is exhibited by a body located at a distance of r_rN2 = 0.022, or r_r1 = 0.33. The velocities of only two bodies exceed the velocity bounds in the graph. The bounds of the graph in terms of body distances are exceeded for about 30 bodies with the largest distance r_rN2 = 11.2 found for body 5336. Its velocity is v = 2.763. Starting from r_rN2 = 4.4, an increasing velocity profile is set:

$v=-0.964+0.333\cdot r_{rN2}$ for $r_{rN2}>4.4$ (40)

The velocity profiles presented in Figure 15 for different structures depending on the numbers of bodies and on their distances are similar, but they differ in quantitative characteristics. These quantitative differences are seen from formulas (39)-(40) for the ascending sections of velocity profiles. From the analysis of velocity profiles over body numbers at the initial time (Figure 15(a), Figure 15(c), Figure 15(e)), it can be inferred that with an increase in the number of layers, the velocity difference between the layers decreases, and the v-profile asymptotically approaches a horizontal line. Moreover, in the 14^th layer (Figure 15(e)) the velocity becomes minimal, v = 2.524, and in the 15^th layer it already increases. This increase in body velocities can also be traced in the velocity profile over distances (Figure 15(f)) at r_rN2≈1 for the final time T = 0.779.

In order to elucidate the behavior of velocity in structures with a large number of layers, structures MS24c99b.dat with 24 layers (Figure 16(c)) and MS34c49b.dat with 34 layers (Figure 16(d)) were created. They have approximately the same number of bodies: N = 29701 and 29156, respectively. In the first structure, the number of bodies in the first layer is N_3,1 = 99, and in the second, N_3,1 = 49. That is why the second structure contains ten extra-layers at an approximately the same number of bodies N.

The images of the structures, obtained using the Galactica program, are shown in Figure 16. They are projections onto the xy-plane of the structures

Figure 16. Multilayer structures in projections onto the horizontal xy-plane as represented by the Galactica program: ((a), (b)) MS15c49b; (c) MS24c99b; (d) MS34c49b. Part b of the figure shows the 15-layer structure at the time T = 0.779 on the two times reduced scale, the lengths of segments are proportional to body velocities, and the point sizes are proportional to their masses.

pre-rotated around the x-axis through an angle of 20˚. The bodies in the structures are shown as points, and the velocity vectors, as segments. In order not to obscure the images of bodies in the 24-layer and 34-layer structures, the velocity vectors are reduced. Condensation of points in these structures occurs in layers whose mid-planes coincide with the line of sight. In the 15-layer structure at the final time T = 0.779 (Figure 16(b)), the enlarged points represent bodies with large masses. Since the graphs show projections of velocities, the small size of segments can be for bodies whose velocity is close to the perpendicular to the xy-plane. Sets of bodies with velocities directed similarly can be constellations.

All structures in Figure 16 were obtained with Inx = 1, i.e. all layers in them were uniformly rotated around the x-axis. As it is evident from the images of the structures, with an increase in the number of layers their shape approaches the spherical shape: the last structure (Figure 16(d)) shows no difference from a sphere.

Figure 17. Generalized velocity profile at the initial time for multilayer structures: 1—MS15c49b, 2—MS24c99b, and 3—MS34c49b.

When analyzing the velocity profiles of these structures, it turned out that those profiles are similar, and they can be generalized as shown in Figure 17. The velocity of bodies v_r1 is plotted along the vertical axis. This is the velocity normalized by the maximum velocity, or the body velocity in the first layer v₁. Plotted along the horizontal axis is the body distance r_r1 normalized by the semi-axis a₁ of the first layer. Numbers 1, 2, and 3 mark the velocity profiles of the 15-, 24-, and 36-layer structures. All those profiles form a general dependence. This dependence also determines the velocity profiles of the 5-layer and 10-layer structures. The latter velocity profiles are not shown in Figure 17.

The orbits of all bodies in the structures under consideration are circles; that is why their velocities can be calculated using expression (9) with r = R_p and v = v_p. From formula (9), it follows that the velocity v is inversely proportional to $\sqrt{r}$ and, according to formulas (6) and (16), they vary in proportion to $\sqrt{m}$ , where m is the mass of bodies inside the sphere of radius r. Therefore, with an increase in r, the velocity decreases unless the mass of the bodies increases so much that its influence becomes predominant.

In the ascending section of the graph in Figure 17, the velocity varies non-linearly, and its average change obeys the law

$v_{r1}=0.333+3.582\cdot {10}^{-3}\cdot r_{r1}$ for $r_{r1}>14$ .(41)

For structures 1, 2, and 3 shown in Figure 17 we have: v₁ = 6.590, 5.396, and 4.485, respectively. These velocities are given in dimensionless form. They change due to the fact that the dimensional velocities are multiplied by the velocity coefficient k_v, which, in accordance with formulas (26)-(27), can be determined as follows:

$k_v=1/\left(A_m{k}_t\right)$ .(42)

The presented values of v₁ vary in proportion to k_v. It follows from here that the generalized velocity profile in Figure 17 is also valid in dimensional form.

8.2. Angular Momenta

The angular momenta of five structures, from the 5-layer MS05c99b structure to the 34-layer MS34c49b structure, are given in Table 2 for the initial time T = 0. This quantity is also called the kinetic momentum. The table shows the projections M_x, M_y and M_z of the angular momentum onto the coordinate axes x, y, z, as well as the absolute value of M. Based on these projections, the angle β₃ between the angular momentum vector M and the z-axis is calculated as follows:

${\beta }_3=arctg\frac{\sqrt{M_x^2+M_y^2}}{M}$ . (43)

In addition, the total orbital momentum of all bodies M_us is given on the condition that the orbits of the bodies lie in the same plane. For circular orbits, the total orbital momentum is defined as follows:

$M_{us}={\displaystyle \sum _{i=1}^N{m}_i{v}_i{r}_i}$ .(44)

The orbital moments of bodies vary from layer to layer. For example, for the

Table 2. Orbital (M_x, M_y, M_z, M, M_us) and rotational (S_sx S_sy S_sz S_s) angular momenta in multilayer structures.

Structures	Orbital angular momentum
	T	Mx	My	Mz	M	β3	Mus	M/Mus
MS05c99b	0	8.690 × 10−3	2.295 × 10−4	0.02459	0.02608	19.5˚	0.02766	0.9428
MS10c99b	0	0.02944	−0.03525	−3.230 × 10−3	0.04603	94.0˚	0.09365	0.4915
MS15c49b	0	0.04884	−0.06214	4.176 × 10−3	0.07915	87.0˚	0.1554	0.5092
MS24c99b	0	0.07142	−0.09615	0.01974	0.12139	80.6	0.2273	0.5341
MS34c49b	0	0.08392	−0.1167	0.03297	0.14747	77.1	0.2670	0.5523

Structures	Rotational angular momentum (spins)
	T	Ssx	Ssy	Ssz	Ss	β3	Ss/M
MS05c99b	1.96	3.741 × 10−9	−7.486 × 10−9	8.795 × 10−9	1.214 × 10−8	43.6˚	4.655 × 10−7
MS10c99b	0.815	−1.036 × 10−7	−1.587 × 10−7	1.277 × 10−7	2.285 × 10−7	56.0˚	4.964 × 10−6
MS15c49b	0.779	5.578 × 10−8	7.778 × 10−7	−4.633 × 10−7	9.070 × 10−7	120.7˚	1.146 × 10−5

34-layer structure, the angular momentum of the central body is 2.19 × 10⁻³³, for the first-layer bodies it is 9.41 × 10⁻⁷, and for the 34^th-layer bodies, 1.46 × 10⁻⁵. It should be noted that the presence of the kinetic moment of the central body is due to its deviation from the center of mass of the structure.

As evident from Table 2, the values of kinetic moments M and M_us increase with the number of layers in the structures. For example, on going from the 5-layer structure to the 34-layer structure the value of M_us increases by almost 10 times. The projections of moments M_x, M_y, and M_z also vary in a wide range, and they can be either positive or negative. The angular momentum vector M is least inclined to the z-axis for the MS05c99b structure: β₃ = 19.5º. In this structure, there are no rotations of layers about the x-axis. Therefore, the angular momentum of the structure M = 0.02608 differs little from the orbital momentum of all bodies M_us: their ratio is 0.9428.

In other structures, the layers are uniformly rotated around the x-axis. This has led to a significant change in the vector M. For instance, for the MS10c99b structure we have: β₃ = 94˚, i.e. the vector M lies practically in the xy-plane. In addition, it is located in the southern hemisphere. Therefore, the total rotation of bodies in the structure proceeds in clockwise direction. In other structures, the total rotation of bodies proceeds counterclockwise.

For the 15-layer structure, we have: β₃ = 87˚, i.e. the angular momentum vector lies even closer to the xy-plane, but it is located in the northern hemisphere. Therefore, the total rotation of bodies proceeds counterclockwise. In these two structures, the momentum ratio M/M_us is close to 0.5. In structures with a large number of layers, this ratio slightly increases, and the angle β₃ decreases, i.e. the vector M moves away from the xy-plane.

In Table 2, the values of the kinetic momentum are given for the initial time. The largest relative change in the angular momentum at the final time, δM_z = 5.8 × 10⁻⁵, was found to occur in the MS15c49b structure. This means that the magnitudes of the momentum remain unchanged to four significant digits. This means that the values of the angular momentum given in Table 2 are also valid for other moments of time.

8.3. Rotational Angular Momentum of Bodies and Their Thermal Energies

In Figure 11, the change in the spins of a body and its thermal energy was considered when other bodies were successively attached to it. Consider now the spins and thermal energies of all bodies at the final stage of structure formation.

Figure 18(a) shows the spin projections S_px, S_py, S_pz and thermal energy E_t for all merged bodies of the 5-layer structure at moment T = 1.96. The numbers of merged bodies i_u are plotted along the horizontal axis. The spin projections S_px, S_py, and S_pz are marked with a cross, a rhombus, and a square, respectively, and the thermal energies are marked with straight segments of a thin line that breaks at body numbers. The break points and the symbols of spin projections falling onto one and the same vertical line with number i_u, belong to the body i_u.

All in all, there are N_u = 43 merged bodies. The spin projections (Figure 18(a)) vary over wide ranges, and the plots show no dominant direction. The spin modulus S_p changes from a minimum value of 2.86 × 10⁻¹⁰ to a maximum value of 2.71 × 10⁻⁹, those values differing by one order of magnitude. The average value of spin modulus S_ppm = 1.21 × 10⁻⁹ in the graph is shown with a dashed line.

The total value of the angular momentum of the entire 5-layer structure, i.e. its spin, is given in Table 2 as projections S_sx, S_sy, and S_sz. The total spin modulus is S_s = 1.21 × 10⁻⁸, i.e. it is exactly ten times greater than the average spin of an individual body S_sm. However, there are N_u bodies in total, and the ratio of the total spin to the sum of the spins of individual bodies is S_s/(N_u∙S_ppm) = 0.23. Thus,

Figure 18. Variation of thermal energies E_t and spins S_p depending on the numbers of merged bodies in three structures, MS05c99e (a), MS10c99b (b), and MS15c49b (c): E_t is the thermal energy; E_tpm is the average thermal energy of the peripheral body; f = 1 × 10⁻⁶ is a multiplier; S_px, S_py, S_pz are the spin projections; and S_ppm is the mean spin modulus of the peripheral body.

the spins of individual bodies compensate each other to 77%.

As evident from Table 2, the total spin of the system is inclined to the z-axis by an angle β₃ = 120.7˚, whereas the kinetic momentum of the system M has an inclination angle β₃ = 19.5˚. The ratio of these momenta is 4.66 × 10⁻⁷, i.e. the rotational momentum is a small fraction of the orbital angular momentum.

The thermal energies E_t of bodies in Figure 18 are reduced by a factor f = 1 × 10⁻⁶. For the 5-layer structure, these energies vary from a minimum value of 5.85 × 10⁻⁵ to a maximum value of 3.24 × 10⁻³, these values differing by 55 times. The average value of the thermal energy of one body is E_tm = 1.08 × 10⁻³, and the thermal energy of all bodies is E_ts = 0.0466.

In the 10-layer structure (Figure 18(b)), there are N_u = 248 merged bodies in total. In this system, there was one merging with the central body, so we will consider the parameters of collisions of peripheral bodies without a central body. The spins of these bodies vary from

2.80 × 10⁻¹¹ to 6.39 × 10⁻⁹ at the mean value of S_ppm = 8.287 × 10⁻¹⁰. In Figure 18(b), the mean value is shown with a dashed line. Here again, the dominating spin direction cannot be identified on the spin projection graphs.

The spin of the central body, S_p1 = 2.13 × 10⁻⁷, is 257 times greater than the average spin of the peripheral body. Simultaneously, the total spin of the system (Table 2) is S_s = 2.29 × 10⁻⁷. Therefore, the projection of the spin of the central body determines the direction of the entire spin of the system. The spin makes an angle β₃ = 56˚ with the z-axis, while the orbital momentum vector is inclined to this axis by an angle β₃ = 94˚. In this case, if one look from the end of the z-axis onto the xy-plane, then the total orbital motion will proceed in clockwise direction, and the total rotational motion, counterclockwise.

The thermal energy of peripheral bodies varies from 5.46 × 10⁻⁶ to 6.45 × 10⁻³, i.e. by 1180 times, with the average thermal energy of one body being equal to E_tpm = 9.44 × 10⁻⁴. The thermal energy of the central body is E_t1 = 0.762, and that of all bodies, E_ts = 0.995. Thus, the central body makes the main contribution to the thermal energy of the system.

The number of merged bodies in the 15-layer structure is N_u = 140 (Figure 18(c)). Here, four bodies have merged with the central body. The spin modules of peripheral bodies vary from 1.73 × 10⁻¹⁰ to 3.75 × 10⁻⁹ at an average value of S_ppm = 1.29 × 10⁻⁹. The S_ppm-value is shown in Figure 18(c) with a dashed line.

The spin of the central body is S_p1 = 8.99 × 10⁻⁷, i.e. it exceeds the spin of the peripheral body by 696 times. Simultaneously, the total spin of the system (see Table 2) is S_s = 9.07 × 10⁻⁷, i.e. it is almost completely due to the central body. The total spin vector and, consequently, the spin of the central body makes an angle β₃ = 120.7˚ with the z-axis. Thus, the total rotational motion with respect to the z-axis proceeds in clockwise direction, and the total orbital motion with angle β₃ = 87 proceeds in an almost vertical plane. The ratio of the total spin to the total angular momentum is 1.15 × 10⁻⁵ (see Table 2).

The thermal energy of peripheral bodies varies from 2.25 × 10⁻⁵ to 6.46 × 10⁻³ with an average value of E_tpm = 1.56 × 10⁻³. The thermal energy of the central body is E_t1 = 0.149, and that of all bodies is E_ts = 0.366. The energy E_t1 is five times lower than the thermal energy of the central body in the 10-layer structure, despite the fact that there occurred four collisions with the central body. This fact indicates that the collision in the 10-layer structure has occurred at a higher velocity of the peripheral body. In this structure, the spin of the central body is four times less than that in the 15-layer structure. This indicates that the collision in the former structure was more frontal, while the collisions in the latter structure were closer to those occurring along tangential lines.

8.4. Constellations in the Structures

In the images of structures at the final time (Figure 4, Figure 7 and Figure 12), we observe close arrangements of bodies to each other involving two, three or more bodies. Some of these arrangements may be due to the coincidence of their visual lines at a sufficiently large distance between the bodies. Another fraction may be due to the short-term convergence of bodies in visual images. However, there may be cases when the bodies form stable associations, that is, constellations. In order to identify constellations in the output files of the Galactica system, a Constns.for program was developed. The output file of Galactica contains the coordinates of bodies at time T. The output files are generated after a set number of integration steps Kl3. The Constns.for program determines the number of bodies with numbers i for which there are bodies with numbers k located at a distance R_ik ≤ d_max from body i. For each body i, one can find up to six bodies numbered k with their distances R_ik. The total number of all nearby bodies k is determined without any restriction, but their numbers and distances R_ik are not memorized.

Figure 19(a) shows the results of calculations performed using this program for the ten-layer structure at the final time T = 0.815 for d_max = 0.003. This value of d_max corresponds to the distance between the bodies on the line of their location at time T = 0 (Figure 7(a)). For choosing this distance, the distances between closely spaced bodies in Figure 7(b) were measured. This distance turned out to be close to 0.003.

Along the horizontal axis, Figure 19(a) shows the numbers of bodies i, and along the vertical line, the numbers of bodies k. All in all, N_is = 1032 bodies i were found among the total number of bodies in the structure, N = 5451. If body i has several close bodies k, this number of bodies is marked in the graph by the number of points along the vertical line of body i. The largest number of bodies i had one close body each, and the smallest, five bodies. There were two such bodies.

It should be noted that the number N_is = 1032 includes all bodies that have a close body. Therefore, if each body has one close body, then there will be 1032/2 = 516 constellations, and if each body has two close bodies, then there will be 1032/3 = 344 constellations, etc.

From Figure 19(a), it is seen that for all bodies, starting from the body i = 1000, the number of close bodies k is concentrated on the diagonal line. In

Figure 19. Numbers of nearby bodies k for bodies i ((a), (b)) and the relative trajectories of bodies in constellations (c) for the ten-layer structure MS10c99b.

Figure 19(b), for bodies i = 4600 - 4800 the graph of nearby bodies is shown on an enlarged scale. The points show the location of bodies at time T = 0.815. It is seen from the graph that on the diagonal line the bodies k have numbers close to those of bodies i. As already noted, the close positions of the bodies determined at the time T = 0.815 may appear occasional. To make sure that these close positions are not occasional, but form constellations, we have to repeat such calculations for another point in time. The results of the calculation for the previous time T = 0.770 are shown with symbols “x”. Most of them coincide with the points for T = 0.815. Out of 41 cases, there are no matches in only three cases. Therefore, we can assume that in 38 cases of close locations, constellations appear. Of the 38 cases, three are different. Body 4746 has a nearby body 4750; body 4747 has two nearby bodies, 4750 and 5409; and body 4750 also has two nearby bodies, 4746 and 4747. This situation was also confirmed for T = 0.770. Therefore, the bodies 4746, 4747, 4750 and 5409 form a constellation of four bodies. It follows from here that in the region of bodies i = 4600 - 4800 there are 17 constellations formed by two bodies and one constellation formed by four bodies. In total, out of 200 bodies i in Figure 8(b) there are 18 constellations. Thus, the total number of constellations can be estimated as 18/200 = 0.09 of the total number of bodies. Then, there must be 0.09∙N = 491 constellations in the ten-layer structure. This estimate is close to the number 516 that was estimated previously from the number of encounters with one body. Therefore, we can assume that up to 9% of all bodies in the structure can form constellations.

The above-mentioned constellation 4746 - 5409 involving four bodies 4746, 4747, 4750 and 5409 is shown in Figure 20. Its position at time T = 0.770 is marked with the velocity vectors of bodies 4746 and 5409. The position of this constellation at the final time T = 0.815 is also shown in Figure 7(b). Evidently, during this time the constellation has moved on the sphere by about 90˚.

In this constellation, body 4747 has double mass, and the rest of the bodies have one mass m₁ each. When viewed on a larger scale, body 4747 is in the center, and the rest of the bodies surround it.

Bodies with a mass of 5m₁ do not form constellations. Most constellations are created with masses m₁. In Figure 19(a), there are four cases of bodies i with mass 4m₁, 13 cases with mass 3m₁, and 76 cases with mass 2m₁. The number of constellations with these bodies will be slightly less, since several bodies with increased masses can enter one constellation.

In the upper right corner of Figure 19(b), four bodies, 4776, 4781, 4785 and 4789, attract attention; one of those bodies, namely, body 4789, has a twice increased mass. It turned out that these bodies form two pairs: 4776 - 4789 and 4781 - 4785. As evident from Figure 20, at time T = 0.770 these pairs occupy position in different places of the 10-layer structure. Based on 20 points, starting from the time T = 0.7167, relative trajectories of the bodies of interest were analyzed: that of body 4776 relative to body 4789 and that of body 4781 relative to body 4785. Figure 19(c) shows the relative trajectories of these bodies, i.e. 4776 and 4781, as projected onto the xy-plane. Relative coordinates are denoted as x_r and y_r. From the initial point I to the final point F, body 4781 makes a little more than three revolutions around body 4789. The average interval between the points is ΔT = 5.195 × 10⁻³. The orbit of body 4781 is an ellipse with semi-axis a = 1.69 × 10⁻³ and eccentricity e = 0.255. The orbital period is P = 0.0325.

The orbit of body 4776, from the initial point I to the final point F, is represented with intersecting segments, since the period of revolution of this

Figure 20. Position of three constellations 4746 - 5409, 4776 - 4789, and 4781 - 4785 in the ten-layer structure MS10c99b at time T = 0.770.

body is comparable with the interval ΔT between the points. The orbit of the body is an ellipse with axis a = 3.12 × 10⁻⁴ and eccentricity e = 0.513. The period of revolution of body 4776 is given by formula (8). It equals P = 3.15 × 10⁻³, a value ten times shorter than the period of revolution of body 4781. During this time interval, body 4776 makes 31 revolutions about body 4789. The position of these two constellations at time T = 0.770 is shown in Figure 20, and at the final time T = 0.815, in Figure 7(b). During this time interval, the motion of these constellations over the sphere has occurred through an angle of not more than 45˚.

Calculations were performed using the Constns.for program with a doubled distance d_max = 6 × 10⁻³. For time T = 0.815, 2650 bodies i that approached other bodies k to a distance d ≤ d_max were identified. In this case, the number of bodies k reached 19 for body i = 93. For time T = 0.770, there were 2715 bodies i with the largest number of bodies k equal to 17 for body i = 84. For this case, the selected first few bodies were not confirmed for the time T = 0.815. Over the interval of change of bodies i = 4600 - 4800 (Figure 19(b)), the cases of discrepancy noted earlier have coincided in this case. Thus, a twofold increase in the distance d_max makes it possible to identify constellations with more distant bodies and with a large number of such bodies.

Similar calculations were also performed for the 5- and 15-layer structures. For the five-layer structure, at the final time T = 1.96, at d_max = 3.16 × 10⁻³, 367 bodies i approaching other bodies k to a distance R_ik ≤ d_max were identified. The number of bodies k reached 5. The distribution of bodies on a graph similar to that of Figure 19(a) was more uniform, and there was no concentration of bodies along the diagonal line.

For the 15-layer structure with d_max = 5.3 × 10⁻³, at the final time T = 0.780, 1843 bodies i that approached other bodies k to a distance R_ik ≤ d_max were identified. The largest number of bodies k was 7. The distribution of bodies in a graph similar to Figure 19(a) was also more uniform, but with lower concentrations of bodies near the diagonal line. The diagonal line could also be traced, but it started from i = 2000.

9. Structure Scaling

9.1. Dimensionless Parameters of Structures

As noted above, in the MLSpStr2.for program, files of structures with initial conditions are created in dimensionless form using parameters m_ss, A_m, and k_t [25]. The dimensionless parameters of the structures are indicated in Table 3. Here, the number of bodies N varies from 1488 to 29,701; the mass of the central body m₀, from 0.8683 to 0.1427; the mass of peripheral bodies m₁, from 9.332 × 10⁻⁵ to 2.532 × 10⁻⁵; the radius of the central body Ra₀, from 2.171 × 10⁻⁵ to 3.965 × 10⁻⁵; the radii of peripheral bodies Ra₁, from 1.220 × 10⁻⁶ to 1.884 × 10⁻⁶; the semi-axis of the orbit of first-layer bodies a₁, from 0.0149 to 0.008587; and the semi-axis of the orbit of last-layer bodies a_N2, from 0.0652 to 0.2427. All structures have a dimensionless mass m_ss = 1 and the same period of revolution of first-layer bodies, P₁ = 0.01. The period of revolution of last-layer bodies P_N2 varies from 0.1046 to 0.7475.

As already noted, the dimensional semi-axes of these structures are identical, and the dimensional masses of peripheral bodies in the 5-, 10-, and 24-layer structures are also identical. In the 15- and 34-layer structures, the masses are twice as large. Yet, the dimensionless semi-axes and the masses of bodies are different in all structures. The passage to dimensionless quantities is related with the masses of structures, and those masses are different. Thus, the identical dimensional parameters in dimensionless form also become different.

The last column of Table 3 indicate the relative outer radius r_95rN2 of the structures whose evolution was studied. This radius was determined at the final time as the largest distance from the center of mass of 95% of the bodies, i.e. this is the size of the structure with 95% confidence. Evidently, for all structures the value of r_95rN2 is approximately the same, with an average value r_95rN2m = 1.662. Only 5% of all bodies are located at a distance from the center larger than r_95rN2m.

Table 3. Dimensionless parameters of the structures used in modeling spherical star clusters: the dimensionless mass of the structure is m_ss = 1; the rotation period of first-layer bodies P₁ = 0.01 is the same for all structures.

Structures	N	m0	m1, *10−5	Ra0, *10−5	Ra1, *10−6	a1, *100	aN2	PN2	r95rN2	Tkin	-U	Rv
MS05c99e	1488	0.8683	8.859	3.965	1.853	1.304	0.0652	0.1046	1.753	1.564	3.105	0.161
MS10c99b	5451	0.6426	6.557	3.586	1.676	1.179	0.1179	0.2535	1.522	2.169	4.260	0.117
MS15c49b	5866	0.4527	9.332	3.189	1.884	1.049	0.1573	0.3911	1.591	2.148	4.222	0.118
MS24c99b	29701	0.2481	2.532	2.611	1.220	0.8587	0.2061	0.5834	-	1.840	3.605	0.139
MS34c49b	29156	0.1427	2.941	2.171	1.282	0.7138	0.2427	0.7475	-	1.539	3.048	0.164

9.2. Scale Transition Algorithm

Based on the dimensionless parameters summarized in Table 3, we calculate the dimensional parameters of five models of globular star clusters. We set the peripheral-body mass to half the mass of the Sun: m_1,m = 0.5 m_S. In the center of globular constellations, the density of stars can reach 100 to 1000 stars per cubic parsec [7]. With this in mind, for all models in dimensional form, we set the semi-axis of the orbit of first-layer bodies to a_1,m = 0.2 pc, where one parsec (pc) is 3.0856776 × 10¹⁶ m.

Taking into account the values of the dimensionless masses in Table 3, the dimensional central-body mass will be

$m_{0,m}=m_{1,m}\cdot m_0/m_1$ ,(45)

and the dimensional mass of the entire structure,

$m_{ss,m}=m_{0,m}+\left(N-1\right)\cdot m_{1,m}$ .(46)

Given the semi-axis of the orbit of first-layer bodies a_1,m, the geometric scale is

$A_{m,m}=a_{1,m}/a_1$ .(47)

With the help of this scale, we find the semi-axis of the orbit of last-layer bodies a_N2,m = A_m,m·a_N2. Given the values of A_m,m and m_ss,m, formula (27) can be used to determine the time factor:

$k_{t.m}=\sqrt{G\cdot m_{ss,m}/A_{m,m}^3}$ .(48)

The velocity coefficient k_v,m is calculated by formula (42). Using these coefficients, converting the angular momentum and energy to dimensionless form, we obtain expressions for their scale factors:

$k_M=\frac{k_{v,m}}{A_{m,m}{m}_{ss,m}}$ ; $k_E=\frac{{\left(k_{v,m}\right)}^2}{m_{ss,m}}$ .(49)

With the help of scale factors, the transition from dimensionless quantities to dimensional ones is carried out. For example, the dimensional period of revolution of outer-layer bodies is calculated as follows: P_N2,m = P_N2/k_t,m.

9.3. Dimensional Parameters of Simulated Globular Clusters

Table 4 indicates the scale factors and dimensional values of globular clusters modeled with five multilayer structures. The body masses are given in solar masses; the semi-axes, in parsecs; and periods are in years. The last column gives the number of stars in a cubic parsec at the initial time. This value is defined as follows:

$n_{pk}=\frac{3N}{4\pi {\left(a_{N2,m}\right)}^3}$ .(50)

The density of stars n_pk in Table 4 changes from 355 for the 5-layer structure to 22 for the 34-layer structure. Since these are average densities, these values correspond to the observed densities of stars in the center of globular star clusters.

The revolution periods of stars in the first layer P_1,m vary from 120 to 170 thousand years (Table 4), and the revolution periods of stars in the last layer P_N2,m, from 1.2 to 12.7 million years. The periods of revolution of stars in constellations are much shorter than in the outer layer. For instance, in the ten-layer structure the period of revolution of body 4781 relative to body 4785 is 388 thousand years (Figure 19(c)), and that of body 4776 about body 4788 is 37.6 thousand years.

Taking into account the average value of r_95rN2m, the diameter of these globular star clusters D_g = 2 × 1.622 × r_95rN2m varies from 3.2 pc for the 5-layer structure to 22 pc for the 34-layer structure. The diameter of observed globular star clusters varies from 20 to 100 pc, and the number of stars in them, from 10⁴ to 10⁵ [7]. Evidently, the globular cluster models based on the 24- and 34-layer structures satisfy these parameters.

The shape of the observed globular star clusters may differ from spherical, so such clusters are characterized by ellipticity [2]. Ellipticity is defined as the difference between the largest and smallest diameters, divided by the largest diameter. In our Galaxy, the clusters NGC 7492, NGC 6144 and NGC 6273 have ellipticity of 0.24, 0.25 and 0.27, respectively [2]. The highest ellipticity of the simulated clusters, equal to 0.6, was exhibited by the 5-layer cluster. For other models, with an increase in the number of layers, the ellipticity approaches zero.

Table 4. Dimensional parameters of globular star clusters modeled with multilayer structures with the semi-axis of first-layer bodies a_1,m = 0.2 pc and star mass m_1,m = 0.5 m_S, where 1 parsec (pc) = 3.0856776 × 10¹⁶ m. The solar mass is m_S = 1.989118 × 10³⁰ kg.

Structures	m0,m, mS	mss,m, mS, *103	Am,m, m, *1017	kt,m, 1/s, *10−15	kv,m, s/m, *10−3	aN2,m, pc	P1,m, years, *105	PN2,m, years, *106	npk, 1/pc3
MS05c99e	4901	5.644	4.733	2.658	0.7948	1	1.192	1.247	355
MS10c99b	4900	7.625	5.234	2.656	0.7192	2	1.193	3.024	163
MS15c49b	2426	5.358	5.883	1.869	0.9095	3	1.696	6.631	52
MS24c99b	4899	19.75	7.187	2.657	0.5236	4.8	1.192	6.957	64
MS34c49b	2426	17.00	8.646	1.869	0.6189	6.8	1.696	12.68	22

It is also noted [7] that deformation tails were found in the shape of all globular clusters for which high-quality optical images were obtained; in other words, there are deviations from spherical symmetry. Therefore, we may assume that the asymmetry of the 10-layer and 15-layer models of globular clusters in Figure 7(b) and Figure 12(b) is their quite acceptable shape.

In order to verify the results of the scale transition, the period of revolution of first-layer bodies was calculated from the dimensional parameters according to formula (8). In formula (8) for a circular orbit we have α₁ = −1, the pericenter radius is R_p = a₁, the pericenter velocity v_p can be calculated by formula (9), and the interaction parameter µ₁, by formula (6). Thus, the periods of revolution calculated up to the units of the fourth decimal place proved to be coincident with the periods calculated using the scaling approach.

The presented parameters of globular cluster models are of interest to scientists who develop cluster models based on an analysis of their brightness. These models are based on the Virial Theorem -U/T_kin = 2, where U is the potential energy of the cluster and T_kin is its kinetic energy. They are defined in a known manner; see for example equation (1.2) in [13]. Dimensionless values of T_kin and -U are given in Table 3, from which it follows that the ratio -U/T_kin varies from 1.960 to 1.986, i.e. is close to 2.

The total mechanical energy E = T_kin + U in the process of evolution of the presented models towards the end slightly decreases in absolute value in the 5-layer structure by 0.02%, in the 10-layer structure by 1.2%, and in the 15-layer structure by 0.87%. In statistical modeling of clusters, one important scale is the virial radius $R_v=G\cdot m_{ss}^2/\left(\text{2}\left|U\right|\right)$ [5]. As can be seen from Table 3, the value of R_v almost completely coincides with the outer radius a_N2 for a 10-layer structure; in other cases, R_v is greater or less than a_N2.

9.4. The Rotation Periods of Stars and Their Temperatures

When considering the approaches and collisions of bodies in the three studied structures, the maximum spins and thermal energies were reported. In addition, the maximum values of velocities v_mx were indicated. With the help of scale factors, we find the dimensional values of the velocities, the periods of rotation of merged bodies P_rt, and the increase in their temperature Δt. For determining the radii of stars in globular star clusters, we specify their average density ρ_gk = ρ_S = 1.406 × 10³ kg/m³, where ρ_S is the average density of the Sun. Then, the dimensional radii of the bodies can be determined from their dimensional masses m_m:

$R{a}_m={\left(\frac{3m_m}{4\pi {\rho }_{gk}}\right)}^{1/3}$ .(51)

The moments of inertia of bodies are defined as for a homogeneous globe:

$J_m=0.4m_m\cdot R{a}_m^2$ .(52)

The dimensional spins and thermal energies are calculated using the scale factors

$S_{p,m}=S_p/k_M$ ; $E_{t,m}=E_t/k_E$ .(53)

Then, the dimensional periods of rotation of stars are defined as follows:

$P_{rt,m}=2\pi J_m/S_{p,m}$ .(54)

For calculating the heating temperature Δt of merged stars, we specify their specific heat equal to the specific heat of water C_t = 1.183 × 10³ J/(kg·deg). Then, the increase in the temperature of merged bodies will be:

$\Delta t=E_{t,m}/\left(m_m\cdot C_t\right)$ .(55)

The results of calculations performed according to these formulas are summarized in Table 5. The first five columns indicate the dimensionless spins, the thermal energies, and the maximum velocities. The values with subscript “0” refer to the central body, and those with subscript “1”, to peripheral bodies. For a peripheral body, the maximum values of spins and thermal energies are given. These bodies with the maximum parameters in the 5- and 15-layer structures have doubled masses, and in the 10-layer structure, the body mass is 5m₁. As already noted, bodies located near the center of structures have a maximum dimensionless velocity v_mx.

In the last five columns of Table 5, the dimensional values are given. The periods of rotation of the central body P_rt0,m vary from 2 to 17 days, while the period of rotation of the peripheral body (P_rt1,m) reach thousandths of a day, which is equivalent to several minutes. The heating temperature of the central body Δt₀ varies from several hundred to two thousand degrees, while that of the peripheral body (Δt₁) reaches several tens of thousands of degrees. In the central region of the clusters, the highest velocity of star motion is equal to two tens of kilometers per second.

With the masses of bodies used in the models of globular star clusters, their relative radii do not coincide with the relative radii at which the problem of interaction of bodies (28) was solved. Their rotation periods and temperatures Δt depend on the body radii. To make sure that the values in Table 5 are reliable, the periods of rotation and temperatures were calculated for the initial dimensional radii of the bodies.

Table 5. Star rotation periods P_rt in days, the collision-induced heating temperatures Δt of stars in Kelvins, and the maximum velocities v_mx in km/s of stars in globular star clusters modeled with multilayered structures (in km/s).

Structures	Dimensionless quantities					Dimensional quantities
	Sp0, *10−7	Sp1, *10−9	Et0	Et1, *10−3	vmx	Prt0,m, days	Prt1,m, days, *10−3	Δt0, K	Δt1, K, *10+4	vmx,m km/s
MS05c99e	-	2.71	-	3.24	14.3	-	1.549	-	2.447	17.99
MS10c99b	2.13	6.39	0.762	6.45	17.7	16.87	1.832	1938	3.215	24.61
MS15c49b	8.99	3.75	0.149	6.46	16.9	1.982	1.085	336.3	3.537	18.58

On average, the temperatures proved to be ten times higher. The temperatures are proportional to the squares of body velocities. At the initial parameters, the velocities are three times higher, which fact explains one order of magnitude higher temperatures.

At the initial parameters, the rotation periods were 60 times longer. The periods of rotation are inversely proportional to body velocities, and they vary in proportion to the squares of body radii. The radii of the bodies with the initial parameters are 15 times smaller. This explains the 60-fold increase in periods.

Thus, in order of magnitude, the increase in the temperature of stars and their periods of rotation in the models of globular star clusters presented in Table 5 adequately reflect these properties of stars.

The processes of merging of bodies depend on their radii. The connection of radii with body masses differs from the connection of distances with masses in the interaction Equation (28). That is why for more accurate modeling of body parameters during the merging of bodies, it is necessary to specify the radii of bodies corresponding to the simulated star system.

The rotation periods and temperatures of peripheral bodies presented in Table 5 are extreme. As it was shown in Section 8.3, the minimum values of spins and thermal energies are one or two orders of magnitude smaller. Therefore, on average, the temperatures Δt₁ for peripheral bodies can be expected to have several times lower, and the rotation periods P_rt1,m, several times larger values.

9.5. Central-Body Models

One of the main problems concerning such stellar associations as globular star clusters and galaxies is the central-body mass: how big should it be? For the problem of interaction of bodies of a plane axisymmetric structure [18] [33], there is an exact solution in the absence of a central body, that is, for m₀ = 0. However, such a structure is unstable. As a result of numerical experiments with a single-layer spherical structure [20] [21], it was found that such structures can be created with a relative central-body mass p_m0 < 0.95. Therefore, multilayer structures were created with certain margins at p_m0 = 0.99.

In a multilayer structure, the central-body mass, related to the mass of the entire structure, decreases with an increase in the number of layers from m₀ = 0.87 for the five-layer structure to m₀ = 0.14 for the 34-layer structure (Table 3). However, the dimensional central-body mass m_0,m remains equal to thousands of solar masses (Table 4). When modeling globular clusters, other researchers also obtain such masses of the central bodies [9]. However, in order to reduce the central-body mass, we will replace it with a model in the form of a multilayer structure similar to the above-considered structures. Let the mass of such a structure be equal to the central-body mass

$m_{ss3,m}=m_{0,m}$ ,(56)

and the outer semi-axis of the orbit is the fraction, equal to k_a1, of the semi-axis of the inner layer of the globular cluster:

$a_{N23,m}=k_{a1}\cdot a_{1,m}$ (57)

Here, the subscript “3” denotes the parameters of the model of the first central body. Subsequently, there will be more models of other central bodies. Two conditions (56)-(57) uniquely determine the parameters of this model.

Table 6. Models of central bodies in the form of multilayer structures under conditions (56)-(57) with k_a1 = 0.5.

Structures	Models of the 1st central body				Models of the 2nd central body
	m03,m, mS	m13,m, mS,	P13,m, years,	PN23,m, years, *104	m04,m, mS	m14,m, mS,	P14,m, years,	PN24,m, years
MS05c99e	4255	0.4341	4045	4.231	3695	0.3770	137.3	1436
MS10c99b	3149	0.3213	1664	4.217	2024	0.2065	23.20	588.2
MS15c49b	1073	0.2306	1534	5.998	474.7	0.1020	14.03	548.8
MS24c99b	1216	0.1240	719.9	4.200	301.6	0.0308	4.346	253.6
MS34c49b	346.1	0.0713	800.6	5.984	49.37	0.0102	3.780	282.6

Based on the parameters of the multilayer structures considered, central-body models were calculated. For a globular cluster model with N₂ layers, a central-body model with the same number of layers was created. Table 6 summarizes the parameters of the models of the 1st central body with k_a1 = 0.5, i.e. the outer size of the central-body model is equal to half the inner size of the multilayer structure. Given in Table 6 are the central-body mass m_03,m and the peripheral-body mass m_13,m (in solar masses), and the period of revolution of inner-layer bodies P_13,m and the period of revolution of outer-layer bodies P_N23,m (in sidereal years). For the five-layer model, the central-body mass m_0,m has decreased from 4901 (Table 4) down to 4255 solar masses (Table 6), i.e., this is an insignificant decrease. For the 34-layer structure, the central-body mass has decreased from 2426 to 346.1 solar masses, that is, by a factor of 7. In this structure, the peripheral bodies have a mass of m_1,m = 0.5 solar masses, and in the central-body model, the peripheral-body mass is m_13,m = 0.0713 solar masses. That is, also 7 times less. The period of revolution of bodies on the inner layer is P_13,m = 800.6 years, and on the outer layer, P_N23,m = 59,840 years. These periods are 212 times shorter than the corresponding periods in the globular cluster model (Table 4).

This central-body model models the central body of a multilayer structure, i.e. this is the model of the first central body. The central-body model also has a central body, i.e. this is a second central body. This body can also be substituted with a multilayer model according to conditions similar to (56)-(57). Table 6 summarizes the parameters of the model of the 2^nd central body; those parameters are marked with the subscript “4”. The mass of the 2^nd central body m_04,m for its 34-layer model has decreased most significantly, namely, also by 7 times with respect to m_03,m. The peripheral-body mass m_14,m has decreased by the same factor. The periods of revolution P_14,m and P_N24,m have also decreased by 212 times compared to the corresponding periods of the model of the first central body.

In these central-body models with two 34-layer structures, the central-body mass is 49.37 solar masses. If we additionally use the model of a 3rd central body, then the central-body mass will be further decreased to seven solar masses. In this case, the central body will be surrounded by 34 layers of revolving bodies with a mass equal to one and a half Jupiter masses. On the inner layer, the period of revolution of bodies will be 6.5 days, and on the 34^th layer, 1.3 years. There will be 29,155 peripheral bodies in total. Then, the same number of bodies with a mass equal to 0.0102 solar mass will be located in 34 layers with periods of revolution of 3.78 years in the 35^th layer and 282.6 years in the 68^th layer. The next 29,155 bodies with a mass of 0.0713 solar mass are located in 34 layers with periods of revolution of 800.6 years in the 69^th layer and 59.84 thousand years in the 102nd layer. The 34 outer layers of this globular cluster contain 29,155 stars with a mass of 0.5 solar mass and periods of revolution of 169.6 thousand years in layer 103 and 12.68 million years in layer 136. All in all, this cluster contains 116,620 peripheral bodies revolving around the central body with a mass equal to 7 solar masses. The radius of the first layer is 0.13 AU, and that of the 136^th layer, 6.8 pc, or 2.06 × 10⁵ AU. As it was already noted, the period of revolution of first-layer bodies here is 6.5 days, and that of bodies in the last layer, 12.68 million years.

Bodies in the inner layer with a mass equal to 1.5 Jupiter masses are not the sources of radiation, so they will shield the central body from light. In layers, starting from the 35^th layer, the mass of bodies will increase, and their luminosity will appear first in the infrared range, then in the red, and then in other spectral ranges.

The scaling of the results of the performed studies makes it possible to foresee such a hierarchically folded multilayer structure of a globular star cluster with a moderate central-body mass. The high velocities of bodies and their close arrangement on the inner layers can contribute to their more frequent collisions than in the structures considered; this circumstance can facilitate the destruction of such structures. Therefore, in order to make sure that such star clusters are possible, it is necessary to create their models and study their evolution. Such models can be created based on the multilayer structures whose evolution from the initially organized state to the steady state has already taken place, for example, based on the 15-layer MS15c49b.dat structure at the time T = 0.779. Then, the study of the evolution of such globular clusters will take much less time.

Unlike in the 34-layer structure, for the 15-layer structure it will be necessary to create not three but several more central-body models. For example, with seven such models, the peripheral-body mass on the inner layer will be equal to 2.2 Jupiter masses. In this case, the period of revolution of such a body will be eight orders of magnitude shorter, which raises doubts about the stability of such a structure.

With three central-body models, the peripheral-body mass in the inner layer will be 51.4 Jupiter masses, its semi-axis will be 1.96 × 10⁻² pc, or 4047 AU, and the period will be 9.13 × 10⁻² years, or 33.4 days. In this case, the central-body mass will be equal to 238 solar masses, i.e. a value ten times less than without the central-body model (Table 4). The globular cluster will contain 23,465 bodies, and the parameters of its 15 outer layers will be such as indicated in Table 4 for the 15-layer structure. Seemingly, such a hierarchical model of a globular cluster will be stable.

10. Discussion

Usually, at modeling of globular star clusters, for example using the NBODY 6 program, the evolution of the shape of the globular cluster is investigated and the change of its statistical characteristics are studied, for example, changes in the distribution of mass along the radius of the cluster. In this case, the internal dynamics of the globular cluster are not considered, the trajectories of the stars are not studied, the processes during their collision are not investigated, for example, the appearance of rotational motion of the stars and their thermal energy, etc.

In the present study, the N-body problem (28) was solved in dimensionless form. Therefore, its results can be applied to stellar associations of different scales, such as planetary systems, globular clusters and galaxies. However, the relative sizes of the bodies in these associations are different. Therefore, the characteristics of processes when bodies collide will be different. In further studies these circumstances will be taken into account.

11. Conclusions

A method has been developed for constructing models of globular star clusters in the form of multilayer spherical structures. In this case, the central body of the structure and the layer surrounding it are adopted as a new central body, around which the next layer is located.

As a result of the interaction of stars, they collide and merge together, and some stars get ejected from the structure. At the first stage, the rate of collisions is kept at a certain level, then it decreases by several times or, sometimes, by several tens of times, and the stage of steady cluster dynamics begins. The stars move in quasi-elliptical orbits around the center of mass of the cluster. The periods of revolution of stars in orbits increase with the increase in orbit sizes. With close passages of stars, the orbits change and their position in space also changes. Sometimes stars collide with each other or with the central body. Individual stars, when passing very close to other stars, can acquire a high velocity and leave the cluster.

Some of the stars moving around the center of mass of the cluster can unite forming constellations. In a constellation, there is a relative motion of stars around its center of mass.

The orbits of stars are located in different planes. Their orbital angular momentum changes due to the interaction, but the angular momentum of the entire cluster remains unchanged. In magnitude, it can reach half the sum of the modules of the orbital momenta of all stars.

When stars merge together, they acquire additional rotation and temperature. The axes of revolution of stars have different directions. The total angular momentum of the rotational motion does not coincide in direction with the total orbital momentum, and in magnitude it is hundreds of thousands times smaller.

The existence of a globular cluster is due to the attraction of bodies to the central body and due to the mutual attraction of stars to each other. With an increase in the number of layers, the relative central-body mass decreases from 0.87 in clusters with five layers to 0.14 in clusters with 34 layers. However, in absolute terms this mass remains at the level of several thousand solar masses. The developed central-body models in the form of multilayer structures show that the central-body mass can be reduced by one or two orders of magnitude. However, the speed of movement of stars on the inner layers becomes large, and this can lead to the destruction of the whole cluster. Apparently, there must be a minimum mass of the central body below which a globular cluster can no longer exist.

Acknowledgements

The evolution of multilayer structures was studied on supercomputers of the Central Shared-Use Center of the Siberian Supercomputer Center, ICM&MG SB RAS, Novosibirsk, Russia.

This research project was carried out under Contract Agreement No. FWRZ-2021-0007 with RAS.

Data Availability

The Galactica system for solving the N-body problem and studying its results is available at . The program for creating multilayer structures and the structure files mentioned in this paper are available at http://wgalactica.ru/smul1/smulski//Data/MLSpStr/.

NOTES

¹http://wgalactica.ru/smul1/smulski/GalactcW/

²http://wgalactica.ru/smul1/smulski//Data/MLSpStr/

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	King, I.R., Hedemann, E.J., Hodge, S.M. and White, R.E. (1968) The Structure of Star Clusters. V. Star Counts in 54 Globular Clusters. The Astronomical Journal, 73, 456-491.[CrossRef]
[2]	Harris, W.E. (1996) A Catalog of Parameters for Globular Clusters in the Milky Way. The Astronomical Journal, 112, 1487-1488.[CrossRef]
[3]	Heggie, D. and Hut, P. (2003). The Gravitational Million-Body Problem. Cambridge University Press. [Google Scholar] [CrossRef]
[4]	Orlov, V.V., Rubinov, A.V. (2008) N-Body Problem in Stellar Dynamics: Textbook. St. Petersburg State University.
[5]	Portegies Zwart, S.F., McMillan, S.L.W. and Gieles, M. (2010) Young Massive Star Clusters. Annual Review of Astronomy and Astrophysics, 48, 431-493.[CrossRef]
[6]	Heggie, D.C. (2014) Modelling and Understanding Globular Clusters. ISIMA.
[7]	Loktin, A.V. and Marsakov, V.A. (2009) Lectures on Stellar Astronomy. Educational and Scientific Monograph. Southern Federal University. (In Russian)
[8]	Llorente de Andrés, F. (2024) Some Old Globular Clusters (and Stars) Inferring That the Universe Is Older than Commonly Accepted. American Journal of Astronomy and Astrophysics, 11, 1-13.[CrossRef]
[9]	Baumgardt, H., Makino, J. and Hut, P. (2005) Which Globular Clusters Contain Intermediate-Mass Black Holes? The Astrophysical Journal, 620, 238-243.[CrossRef]
[10]	Mackey, A.D., Wilkinson, M.I., Davies, M.B. and Gilmore, G.F. (2008) Black Holes and Core Expansion in Massive Star Clusters. Monthly Notices of the Royal Astronomical Society, 386, 65-95.[CrossRef]
[11]	Cuevas-Otahola, B., Mayya, Y.D., Puerari, I. and Rosa-González, D. (2020) Mass-Radius Relation of Intermediate-Age Disc Super Star Clusters of M82. Monthly Notices of the Royal Astronomical Society, 500, 4422-4438.[CrossRef]
[12]	Tully, R.B., Rizzi, L., Dolphin, A.E., Karachentsev, I.D., Karachentseva, V.E., Makarov, D.I., et al. (2006) Associations of Dwarf Galaxies. The Astronomical Journal, 132, 729-748.[CrossRef]
[13]	Aarseth, S.J. (2003). Gravitational N-Body Simulations. Cambridge University Press. [Google Scholar] [CrossRef]
[14]	Baumgardt, H. (2016) n-Body Modelling of Globular Clusters: Masses, Mass-to-Light Ratios and Intermediate-Mass Black Holes. Monthly Notices of the Royal Astronomical Society, 464, 2174-2202.[CrossRef]
[15]	King, I.R. (1966) The Structure of Star Clusters. III. Some Simple Dvriamical Models. The Astronomical Journal, 71, 64-75.[CrossRef]
[16]	McLaughlin, D.E. (2000) Binding Energy and the Fundamental Plane of Globular Clusters. The Astrophysical Journal, 539, 618-640.[CrossRef]
[17]	Smulsky, J.J. (2019) The Upcoming Tasks of Fundamental Science. Sputnik+ Publishing House. (In Russian) http://www.ikz.ru/~smulski/Papers/InfPrZaFN.pdf
[18]	Smulsky, J.J. (2003) The Axisymmetric Problem of Gravitational Interaction of N-Bodies. Math Modeling, 15, 27-36. (In Russian) http://www.ikz.ru/~smulski/smul1/Russian1/IntSunSyst/Osvnb4.doc
[19]	Smulsky, J.J. (2015) Exact Solution to the Problem of N Bodies Forming a Multi-Layer Rotating Structure. SpringerPlus, 4, Article No. 361.[CrossRef] [PubMed]
[20]	Smulsky, J.J. (2016) Distributed Structures on the Sphere. Deposited in VINITI 22.08.2016, No. 112-V2016. (In Russian) http://www.ikz.ru/~smulski/Papers/SphDsSt2.pdf
[21]	Smulsky, J.J. (2019) Periodic Orbits of N Bodies on a Sphere. Cosmic Research, 57, 459-470.[CrossRef]
[22]	Smulsky, J.J. (2015) Multilayer Coulomb Structures: Mathematical Principia of Microcosm Mechanics. Open Access Library Journal, 2, 1-47.[CrossRef]
[23]	Smulsky, J.J. (2012) Galactica Software for Solving Gravitational Interaction Problems. Applied Physics Research, 4, 110-123.[CrossRef]
[24]	Smulsky, J.J. (2018) Future Space Problems and Their Solutions. Nova Science Publishers. http://www.ikz.ru/~smulski/Papers/InfFSPS.pdf
[25]	Smulsky, J.J. (2012) The System of Free Access Galactica to Compute Interactions of N-Bodies. International Journal of Modern Education and Computer Science, 4, 1-20.[CrossRef]
[26]	Smulsky, J. (2014) Module of System Galactica with Coulomb’s Interaction. International Journal of Modern Education and Computer Science, 6, 1-13.[CrossRef]
[27]	Smul’skii, I.I. and Krotov, O.I. (2015) Change of Angular Momentum in the Dynamics of the Solar System. Cosmic Research, 53, 237-245.[CrossRef]
[28]	Smulsky, J.J. (2019) Angular Momentum Due to Solar System Interactions. In: Gordon, O., Ed., A Comprehensive Guide to Angular Momentum, Nova Science Publishers, 1-40. http://www.ikz.ru/~smulski/Papers/CGAngMom1_2Cv.pdf
[29]	Laskar, J. (1996) Marginal Stability and Chaos in the Solar System. Symposium—International Astronomical Union, 172, 75-88.[CrossRef]
[30]	Laskar, J., Correia, A.C.M., Gastineau, M., Joutel, F., Levrard, B. and Robutel, P. (2004) Long Term Evolution and Chaotic Diffusion of the Insolation Quantities of Mars. Icarus, 170, 343-364.[CrossRef]
[31]	Melnikov, V.P. and Smulsky, J.J. (2009) Astronomical Theory of Ice Ages: New Approximations. Solutions and Challenges. Academic Publishing House “Geo”. http://www.ikz.ru/~smulski/Papers/AsThAnE.pdf
[32]	Smulsky, J.J. and Smulsky, Y.J. (2012) Dynamic Problems of the Planets and Asteroids, and Their Discussion. International Journal of Astronomy and Astrophysics, 2, 129-155.[CrossRef]
[33]	Smulsky, J.J. (1999) The Theory of Interaction. Publishing House of Novosibirsk University. (In Russian) http://www.ikz.ru/~smulski/TVfulA5_2.pdf
[34]	Smulsky, J.J. (2004) The Theory of Interaction. Cultural Information Bank. http://www.ikz.ru/~smulski/TVEnA5_2.pdf