**1. Introduction**

Numerical models are often used for solving complex problems in many fields [1–5]. They allow researchers to observe complex systems in artificially created scenarios. This makes it possible to identify the major governing parameters and study the underlying physics. Using these observations, one can later develop sophisticated quantitative models. To be efficient, a heuristic model must account for important features of the phenomenon and its main governing factors.

We often consider it useful to construct mechanical systems from some interacting blocks (or 'particles'). The interactions between them must fulfill certain generic properties, e.g., for solids they should enable the system to sustain its shape, and this should be implemented in the model. The atomic lattice of crystalline solids is responsible for the long-range order in a solid. Its existence is owing to the physics of interatomic interactions. Such a description has general validity and also applies to amorphous solids and even fluids. It also applies to fracture processes, which occur via the subdivision of a solid into fragments before separation. The analogous approach can be applied to describe various phenomena and systems, such as pattern formation in nature, the mixing of powders, and the behavior of dynamic systems (see, for example [6,7]).

At present, it is still not possible to advance to satisfactory length and time scales by using classical molecular dynamics that operate with atoms or molecules. That is why representative objects are often used as particles, whose motion, in its entirety, reflects the pertinent aspects of the behavior of the system. For example, the moveable cellular automata [1], smooth particle hydrodynamics [2], and discrete automata methods [3,6,7] can be mentioned. Note that visualization and processing of the results of the computations are very involved [7]. In this context, the reduction in computation costs is essentially important if one needs a relatively long-term simulation to obtain physically interesting results.

Below we will present one of them, which applies (virtually or naturally) cylindrical and spherical symmetries to reduce the calculation time and memory consumption. Our concern in the present work is primarily the extraction of qualitative information from numerical simulations, with convenient access to direct visualization for various scenarios and over long durations. To do this, we used the isotropic 'moveable automata' technique (cf., e.g., [3,7]), which involves a weak, long-range interaction between particles in addition to a short-range repulsive one.

The form of the interparticle interactions used ensures the existence of a minimum overall potential and the attendant equilibrium state. Such a system does not require 'walls' created by the boundary conditions to form compact fragments of a crystalline lattice. Even in the case when the average density of the material is lower than that required for contiguous space-filling, the particles show a tendency to aggregate [7,8]. With these considerations in mind, we use a model whose details are given below. In solving the equations of motion, we confined ourselves to problems where the thermal energy *kBT*, where *T* is the absolute temperature and *kB* is the Boltzmann constant, is negligibly small compared to other energy terms involved. However, as is usual in dynamics, mutual transformations between potential and kinetic energy still exist. It can be proven [9] that in the limit of high kinetic energy density, *Ekin* >> *Cij*, where *Cij* is a stiffness parameter, the system behaves similar to a gas of nearly free-flying particles. In the opposite limit case, *Ekin* << *Cij*, a crystal lattice is formed spontaneously.

The nodes of the lattice oscillate around the bottoms of the potential valleys formed by the neighboring particles. It can be described in terms of a harmonic Hamiltonian, slightly perturbed by non-linear terms. In this sense, the system deforms elastically. While in the gaseous state, the system is isotropic on average, and in the solid state, its lattice has a hexagonal symmetry if the interaction potential between the particles is isotropic. In both limited cases, the system exhibits regular dynamic behavior. However, due to a difference in symmetry, there is a transition between the two states via a mixed disordered state, which, in the presence of dissipation, represents a viscous fluid. Thus, the minimalist model considered accounts for the occurrence of all three aggregate states of the system.

#### **2. Numerical Model**

The model describes a system of *N* particles represented by the vector radius **r***i*, the momentum **p***i*, and the interaction potential *U*( + +**r***i* − **r***j* + +) corresponding to the following Hamiltonian:

$$H(\mathbf{r}\_i, \mathbf{p}\_i) = \sum\_{i=1}^{N} \mathbf{p}\_i^2 / 2m\_i + \sum\_{i,j=1}^{N} \mathcal{U}(|\mathbf{r}\_i - \mathbf{r}\_j|) / 2. \tag{1}$$

It is opportune to represent the interaction potential by a pair of Gauss potentials:

$$\mathcal{U}(\left|\mathbf{r}\_{\bar{i}}-\mathbf{r}\_{\bar{j}}\right|) = \mathbb{C}\_{\bar{i}\bar{j}} \exp\left\{-\left[\left(\mathbf{r}\_{\bar{i}}-\mathbf{r}\_{\bar{j}}\right)/c\_{\bar{i}\bar{j}}\right]^2\right\} - D\_{\bar{i}\bar{j}} \exp\left\{-\left[\left(\mathbf{r}\_{\bar{i}}-\mathbf{r}\_{\bar{j}}\right)/d\_{\bar{i}\bar{j}}\right]^2\right\},\tag{2}$$

where *Cij* and *Dij* define the magnitude, while *cij* and *dij* are the radii of attraction and repulsion, respectively. The minimization condition corresponding to an equilibrium reads as follows: *Cij* >> *Dij*, *cij* < *dij*. Since the parameters *Cij*, *Dij*, *cij*, and *dij* form arrays covering all possible combinations, it can be said, without loss of generality, that the equations of motion

$$
\partial \mathbf{m}\_i \partial \mathbf{v}\_i / \partial t = -\partial H(\mathbf{x}\_i, \mathbf{p}\_i) / \partial \mathbf{p}\_i = \mathbf{f}\_i. \tag{3}
$$

describe systems with any number of components, including one- and two-component systems, which are of interest in this paper. There is a requirement to specify the parameters *C*11, *C*22, *C*12 = *C*21, etc.

In standard numerical experiments, the particles are initially placed at random on a rectangle with the side lengths of [0, *Lx*] and [0, *Ly*]. In addition, in the case of contact with a thermal bath, the boundary conditions need to be accompanied by additional ones. If the temperature *T* of the thermal bath at a given boundary (or within the system) is non-zero, an array of random *δ*-correlated Langevin forces need to be included. These are to satisfy the fluctuation-dissipation theorem in the form of

$$
\langle \langle \mathfrak{f}(t, \mathbf{x}\_i, y\_i) \mathfrak{f}(tt, \mathbf{x}\_j, y\_j) \rangle \rangle = D \delta(t - t') \delta\_{ij} \langle \mathfrak{f}(t, \mathbf{x}\_i, y\_i) \rangle = 0,\tag{4}
$$

where *D* = 2*γkBTm*/*dt*, *kB* is the Boltzmann constant, *γ* is a dissipation constant, and *dt* is the discrete time step.

Independently of the thermal source, the interaction between the particles produces kinetic energy and causes effective thermalization of the system. So, one has to introduce a dissipation channel acting to equilibrate the relative velocities of particles that happen to be within the distance *cv* close to the equilibrium. This can be performed by additive force

$$f\_i^{\mathbf{v}} \sim \sum\_{j=1}^{N} \left(\mathbf{v}\_i - \mathbf{v}\_j\right) \exp\left[-\left[\left(\mathbf{r}\_i - \mathbf{r}\_j\right)/c\_v\right]^2\right] \tag{5}$$

acting on the *i*-th particle, with ye<sup>t</sup> another dissipation constant *η*. As a result, the equations of motion assume the following form:

$$m\_i \frac{\partial \mathbf{v}\_i}{\partial t} = f\_i^{\mathbf{r}} - \eta f\_i^{\mathbf{v}} - \gamma \mathbf{v}\_i + \xi(t, \mathbf{r}\_i) \tag{6}$$

These equations can be integrated using Verlet's method, which conserves the energy of the system at each time step, provided no energy is supplied externally through mechanical work or temperature variation.

#### **3. Artificial Cylindrical Symmetry**

The simplicity of the equations of motion is deceptive. Summation occurs at every time step over all possible neighbors, including very distant ones. At the same time, to obtain a more or less realistic model of the system and exclude the effect of boundary conditions that are too close, one normally needs a sufficiently large number of "particles". It becomes especially important when a long-duration run of the numerical simulation (which often reflects the long duration of the real physical experiment) is needed.

For example, this problem appears when one models a process of strong shear deformation, which causes plastic transformations of the material, which results in so-called soli-state turbulence [9]. If the numerically generated system is too short, movable automata quickly reach a boundary of the specimen in a direction of shear deformation and repulse from it back to the internal regions of the system. It completely modifies the numerically obtained results, which strongly deviate from the picture of the real process. To avoid such an artifact, one has to apply as remote a boundary in the direction of shear as possible. To conserve a reasonably limited number of the automata in the numerical procedure, one has to reduce the width of the system in the orthogonal direction. A system that is too narrow with close boundaries in the orthogonal direction also gives some uncontrollable corrections.

In such a case, it is convenient to employ periodic boundary conditions on the horizontal axis. Accordingly, a particle leaving the interval [0, *Lx*] is returned to it at the opposite end, and the vectors **r***i* − **r***j* connect particles located within the interval [0, *Lx*] or the images of escaped particles at the opposite side of the system. On the vertical axis, the system is limited by plates at *y* = 0 and *y* = *Ly* with reflecting boundary conditions: *Uup* = *C* exp[(*y* − *Ly*)/*c*] and *Udown* = *C* exp[−*y*/*c*]. The stronger the inequalities *C* >> *Cij* and *c* << *cij*, the more rigid and sharp the walls are in relation to other forces and lengths relevant to the system.

The conceptual picture of such a system is shown in Figure 1. It is seen directly that two ends of the stripe in the horizontal directions are glued without a break. The particles mostly smoothly pass any (imaginary) gap between two ends of the stripe and continue their motion along a cylindrical surface. In the formal numerical model, the most important

thing is that the particles from one side of the stripe interact with the "images" of the particles from another one. However, such a trick is not needed if the particles are placed on the surface of a real cylinder in 3D space.

**Figure 1.** Conceptual picture of a system with periodic boundary conditions in one-dimensional and repulsing boundaries in the orthogonal one. The origin of the particular system and meaning of the colors in different subplots are described in the main text.

In any case, such effective "infiniteness" of the system provides an advantage (sometimes huge) in computer time consumption, because one does not need an extremely long system to minimize the effect of the boundaries, whereas the number of operations grows and squares *N*<sup>2</sup> of the number of particles *N*.

Development of the instant structure presented in Figure 1 with the duration is reproduced in the Movie\_01.avi [10]. It is seen directly how the particles continue their motion "infinitely" along the cylindrical surface without any disturbance passing any vertical line in subplot (c), including formally connected left and right boundaries of the planar pictures in subplots (a) and (b).

This advantage of the "infiniteness" becomes extremely important when in the extremely long-duration run, one needs to obtain the formation of a nontrivial self-organized system. In particular, mutual collisions of the ordered domains accompanied by their consolidation during the shear of a mixture of two compounds between the plates moving in opposite directions is used in Figure 1 and the corresponding Movie 1.avi [10] as a representative example. This is interesting as a particular realization of solid-state turbulence is discussed in recent publications (see [9] and references therein).

Besides the mathematically formal procedure convenient for the numerical simulations, the described model can correspond to the physically realizable cylindrical system, which potentially gives figurative feedback to the natural experiment, which can reproduce the quasi-infinite run under limited conditions of the laboratory or in some real systems. We will discuss such possibilities below.

#### **4. Practical Example with Stick-Slip Behavior of a System under Shear**

For practical needs and a further description, it is useful to calculate a measurable quantity: The external force acting on the system. First, we calculate the work produced by the forces acting on each particle of the system or the corresponding power *pi* = **f***i***v***i*. Summation over all particles yields the total power

$$2VF\_{\text{ext}}(t) = P(t) = \sum\_{i=1}^{N} \mathbf{f}\_i \mathbf{v}\_i \tag{7}$$

expended by a pair of external forces to maintain the motion of the boundary plates with a constant velocity, *V* = *const*. This yields the total force, *Fext*(*t*). Calculating spatially distributed power *pi* = **f***i***v***i* enables one to observe localization of excitations and hence the fragments of the lattice practically move as a whole.

Namely, these colorized spots are clearly seen in some places in Figure 1. They definitely appear in different places and time moments and travel as energy waves along the system following the "events" corresponding to strong changes of the lattice due to mutual collisions of the ordered domains, and so on. It certainly causes avalanches of the energy consumption after stick-slip behavior of the total system under external shear.

Watching the Movie\_01.avi, one can observe how (initially randomly placed) relatively small islands of the rigid component collide with one another and form larger ordered "domains", moving and even rotating as a whole. However, they also can split into smaller ones again. The long-duration run allows us to accumulate statistics and observe that the distribution of the domain sizes becomes stationary. Another opposite process is normally organized in the experimental studies [9], where originally, a stripe of the hard material is placed between two soft ones. In this case, it deforms (elastically and plastically) and then splits into a set of islands moving and rotating as a whole. Their distribution tends to be the same equilibrium (stationary, but not static) as the opposite limit. The early infant stage of this process is illustrated in Figure 2, accompanied by the Movie\_02.avi [11], reproducing the initial stage of the process in dynamics. One can see the plastic deformation (and destruction into the separate parts later on) of the rigid layer was definitely accompanied by local excitations and even avalanches of them.

**Figure 2.** Localized excitation moving along rigid stripe inserted between two soft layers. It is seen directly that velocity of the rigid stripe inserted between two soft ones (green color) is almost constant and much smaller than velocities of the outer layers moving in opposite directions (blue and red colors, respectively).

From academic and practical points of view, it is interesting to learn how this process of a solitary rigid plate (without soft outer ones) behaves itself under shear. This behavior is reproduced in the Movie\_3.avi [12].

All the mentioned processes are clearly seen in dynamics in the sequence of the movies Movie\_03.avi, Movie\_04.avi, and Movie\_05.avi [12–14] where direct observation is accompanied by the curves giving quantitative information during the long (potentially infinite) run. In particular, one can see how the time–space map of the horizontal velocity of a selected vertical cross-section recording is static. It records information about the history of the events (wave propagation and their avalanches) in the system. For a two-component system, one can numerically analyze images from every frame of the movie and calculate the time depending on the number of domains in the rigid phase, the histogram of the area occupied by the domains of a given size, as well as the size of the maximal and averaged domain shown to accumulate during the long-duration run. Everything is reproduced directly in the movies. Analogous information can be accumulated in the static form and reproduced, as the examples in Figures 3–5 show below.

**Figure 3.** Square root *h* (|*pi*|) 1/2 of the histogram *h*(|*pi*|) (**a**) and the time dependence of the fraction of the points for which the inequality |*pi*|<sup>&</sup>gt; *p*∗ holds (**b**). Note the bursts in the magnitude of this fraction corresponding to the peaks in the time dependence of the total force *Fext*(*t*) = *<sup>P</sup>*(*t*)/*V*, which is presented in (**c**).

Let us specifically mention once more that sliding and jumping over a number of contacting "lattice planes", which have different orientations, results in a phenomenon that can be treated as a kind of "stick-slip" behavior in the system under consideration.

Besides direct visualization, it is useful to characterize the excited sites by the absolute value of the power, |*pi*|=|**<sup>f</sup>***i***v***i*|. For that, we plot a histogram *h*(|*pi*|) for the probability of finding a given value of |*pi*|. As expected, the function *h*(|*pi*|) has a pronounced maximum at low values of |*pi*| associated with small perturbations of the system almost everywhere within the ordered domains. The regions of large strain contribute to the total power 2*VFext*(*t*) = *P*(*t*) = *N* ∑ *i*=1 **f***i***v***i* most significantly, but the partial contribution of an individual region is not very significant. As a result, the distribution *h*(|*pi*|) has an extended 'tail' at large |*pi*|.

**Figure 4.** Time–space map of the horizontal velocity of selected vertical cross-section *vx*(*<sup>t</sup>*, *y*) recording static information about history of the events (wave propagation) in the system.

**Figure 5.** Time-dependent number of rigid phase domains, histogram of area occupied by the domains of the given size, shown accumulated during long-duration run, as well size of the maximal and averaged domain.

For both the locations of the maxima and the asymptotes to be well-resolved at the same time, a nonlinear function *h* (|*pi*|) 1/2 is plotted in Figure 3. It can be roughly divided into two regions—left and right of a vertical dashed line that marks a certain characteristic value *p*<sup>∗</sup>. The fraction of points where |*pi*|<sup>&</sup>gt; *p*∗ holds, i.e., those points where intensive processes occur, evolves with time. The results of this procedure for a particular realization of the deformation process are shown in Figure 3b. One can recognize pronounced bursts of activity, which can be associated with the periods of emergence and development of rebuilt structures localized at the edges of domains (microscopic analogues of the earthquakes) seen in Movie\_3.avi [12]. The same periods can also be identified in the time dependence of the total picture presented in Figure 3c. To facilitate the comparison, some of the bursts are marked with dash-dotted lines.

The function *Fext*(*t*) represents the time dependence of the force required to slide a solid on a surface. From this dependence, the yield force, or yield stress, can be identified. In our computer experiments, this dependence resembled that of the resistance force for sliding a solid on a surface with random roughness in the stick-slip regime.

#### **5. Realistic Quasi-Infinite Objects**

Now let us discuss some physical objects with cylindrical and spherical symmetries, which can serve as examples of effectively infinite systems. First of all, quite often, longrun shear-force deformations are organized by means of rotating cylindrical systems. The typical configuration of such a system is illustrated in Figure 6 below. Here, an experimental material is placed between two cylindrical boundaries. The central one is motionless and the external one is rotating. It causes a configuration where different layers of the material inside tend to rotate with the velocity, which gradually increases from *<sup>V</sup>*(*<sup>R</sup>*min) 0 to *<sup>V</sup>*(*<sup>R</sup>*max) *Vext*. If both radii are large enough (and in following the system of the coordinates), it is mechanically equivalent to the motion between two plates moving with the velocities *V*max = − *Vext*/2 and *V*max = *Vext*/2. Another interesting example of a naturally cylindrical system can also be seen in biology [9].

**Figure 6.** Instant configuration of the cylindrical quasi-infinite system presented in the Movie\_6.avi [15].

This fact is clearly seen in the left subplot of Figure 6 where the particles are colored depending on their velocity (from maximum red to minimum blue, respectively). Movie\_06.avi [15] reproduces the process of formation and growth of the rigid domains inside the two-component mixture and the tendency towards the dynamic equilibrium distribution of their sizes. The process strongly resembles the same process in the abstract cylindrical system studied above. It also can be continued, as far as necessary.

The only (but very important) difference that such a system can show exists in reality. It even has its own advantage that the numerical procedure in such a configuration does not need fictional images. So, it is potentially expected to require less computer memory and be less time-consuming. However, it still has the disadvantage that due to the geometry, the inner and outer layers of the system are not exactly equivalent. To obtain better equivalence, one has to use both relatively large external and internal radiuses. It actually demands a relatively large number of mutually interacting particles (movable automata).

Spikes of the interaction forces in the places where the domains collide or strongly deform are seen directly in the middle subplot of Figure 6. There is one more interesting characteristic of the system that has to be elucidated. The system definitely rotates as a whole. At the same time, different regions mutually interact and can rotate quasiindependently. This effect is practically the same as that which can be observed in the abstract quasi-infinite system above.

To separate integral and differential rotations, it is convenient to directly calculate the local circulation (below, for briefness, we call it 'rotor' and denote it as *rotj*) of the particles in the proximity of every particle of the system and colorize them accordingly. According to this local circulation, the particle in position **r***j* is given by the following equation:

$$rot\_{j} = \sum\_{k}^{N\_{\text{prov}}} \left[ (\mathbf{v}\_{k} - \mathbf{v}\_{j}) \times (\mathbf{r}\_{k} - \mathbf{r}\_{j}) \right] / (|\mathbf{r}\_{k} - \mathbf{r}\_{j}|). \tag{8}$$

Differential rotation is visualized in the right subplot of Figure 6, and is especially clearly seen in the dynamics of the corresponding Movie\_06 [15]. One can distinguish plenty of places where the particles surrounding a given particle instantly rotate in the direction opposite to the rotation of the complete system. This effect corresponds to the "solid-state turbulence", which was already observed in a system confined between two planar plates moving in alternative directions.

The mentioned rotation is even more pronounced if the space between the internal and external cylinders is not so densely filled by the soft component. In this case, the "rocks" formed by the hard component aggregate faster and have more empty space to rotate as a whole. An example of such a system is reproduced in Movie\_07.avi [16]. Here, one can clearly observe the rotation of the solid "rocks" accompanied by the motion of their sides colored by blue in the left subplot, which marks the particles playing a role in the instant centers around which their neighbors locally rotate in the opposite direction. Besides, the collisions between the rocks are clearly marked by the flashes of force colors shown in the central subplot.

Let us note in advance that the visualization by means of the colorized rotor is also very useful for the study of rotating systems with spherical symmetry. In the context of the present paper, such symmetry is another realization of a quasi-infinite system existing in nature.

It appears when a spherical belt near the surface of the sphere can rotate more or less independently of the rotation of the whole object, for example, when it is an atmosphere or lithosphere of a planet. In both cases, the belts on the surface are partially involved in the rotation either by an effect (in some senses, just physical friction) on the solid surface of the planet or by the continental flotation on magma.

It is important to mention that it causes practically the same configuration of the velocities, which was studied above for the shear-induced turbulence. In reality, the velocity of the spherical belt monotonously depends on the latitude. It is maximal near the equator of the planet and has to decrease down to zero near the poles. In fact, one has a doubled picture of the shear-induced motion with two stripes mirroring one another.

Despite the complication related to the principally 3D configuration, it paradoxically gives an advantage from the computing point of view. It does not need phantom images of the particles, because the boundary does not exist at all. Certainly, it allows the simulation to run infinitely without any limitation from the boundary conditions. So, one can obtain rather realistic results using a relatively small number of movable automata.

The only thing one needs is to add an attraction (gravitation in fact) of the particles to the solid spherical surface. In principle, this force is a combination of the Newtonian gravitation potential *Uattract* ∼ 1/<sup>+</sup> +*rj* − *R* + + and exponential repulsion form the surface *Urepuls* ∼ exp[(*rj* − *<sup>R</sup>*)/*c*]. However, in the particular model, one needs only to "glue" the particles to the surface, so the cheapest way is to apply sufficiently strong confinement with the potential *Uconf*∼ (*rj* − *R*) 2 .

Generally speaking, in the case of a spherical surface, the definition *rotj* = *Nproxy*

∑ *k* [(**<sup>v</sup>***k* − **v***j*) × (**<sup>r</sup>***k* − **<sup>r</sup>***j*)]/(|**<sup>r</sup>***k* − **<sup>r</sup>***j*|) produces a vector in 3D space. To associate it with

the vortices and anti-vortices in the spherical belt (of atmo- or lithosphere) we visualize its projection perpendicular to the surface with the coordinate **r***j* and colorize it using the standard MatLab map. It is expected that, due to the general rotation of the 'planet', in the case of ideal laminar flow, the picture will have generally blue and red belts of color on the southern and northern hemispheres and a yellow–green stripe around the equator. However, if the vortices do exist, they will appear as spots of alternative colors (red in southern and blue in northern hemispheres, respectively).

For a trial, we start from the surface randomly covered by the relatively soft material (viscous "atmosphere") and allow it to self-organize due to the rotation. As expected, we obtain something similar to layered rotation with some vortices appearing between the layers. This picture reminds us of a Jupiter-like structure. Figure 7 reproduces a sequence of early configurations of the self-organization in a spherical quasi-infinite system. The first reproduces the very beginning of the process, starting from the randomly distributed particles on the sphere (a), followed by two intermediate ones (b)–(c), and the last subplot (d) reproduces a picture close to the stationary moving configuration. In this stage, the parallel belts of the vortices and anti-vortices in the "equatorial" region (shown by red and blue colors, respectively) are already seen.

**Figure 7.** Sequence of the instant configurations of self-organized spherical quasi-infinite system. Early (starting from the randomly distributed particles on the sphere), two intermediate, and close to the stationary moving configurations are shown in the consequent subplots (**<sup>a</sup>**–**d**), respectively. Parallel belts of the vortices and anti-vortices in the "equatorial" region corresponding to the positive and negative rotor projection perpendicular to the spherical surface (shown by red and blue colors, respectively) are seen directly.

Depending on the relations between the interaction and damping coefficients, as well as the speed of the rotation, one can obtain a different stationary picture. The vantage of the model is that, in principle, one can run it for as long as necessary to obtain rich stationary behavior. We performed a number of such experiments and obtained sometimes very realistic patterns, which closely resembled recently observed images of planetary atmospheres. See, for example, Movie\_08.avi [17].

It is important to note that (as it already was observed in the 2D case), depending on the values of interaction amplitudes *Cij*, *Dij* and distances *cij*, *dij* the equations of motion can describe either soft (liquid) materials, solid ones, or even a mixture of solid and liquid ones. In the context of a spherical system, it provides an intriguing opportunity to model a flotation of the "continental plates" involved in the motion of a rotating liquid "magma".

Certainly, we understand that the model cannot cover all the aspects of the extremely long and multifactorial planetary formation and tectonics. However, in the context of the present conceptual work, one can claim that the model has its own advantages. It is simple, very efficient, and allows for studying a large number of the "theoretical toys" with different relations between all the constants, randomly chosen initial conditions, and scenarios, including very exotic ones. Despite the wide variety of these scenarios, we already observed more or less typical invariants.

One of the typical intermediate stages created by the process reproduced in Movie\_09.avi [18] is shown in Figure 8, where the spatial distribution of the rigid and soft (brown and blue, respectively) components is reproduced in subplot (a). Corresponding to the present configuration, distributions of the local velocities and forces are shown in subplots (b) and (c), respectively. In particular, deep red spots of the local forces demonstrate some kind of "earthquakes", which normally appear when the "continents" either collide one with another or the "continents" are torn apart.

**Figure 8.** Typical intermediate configuration of the two-component spherical system reproduced in Movie\_09. Spatial distribution of the rigid and soft (brown and blue, respectively) components is reproduced in subplot (**a**). Corresponding to the present configuration, distributions of the local velocities and forces are shown in subplots (**b**,**<sup>c</sup>**), respectively. A tendency of the concentration of the "continents" to "equator" at fixed axis of the rotation is seen directly.

Our observations of the different scenarios testify to the general tendency of the concentration of the "continents" to the "equator" at a fixed axis of rotation. As a rule, after the long run, they finish in either a kind of "Pangea", or perhaps in two to three "continents" of very different sizes. It partially resembles the real geological history of the Earth, but partially differs from it, because it is known that "Pangea", over time, broke into several smaller parts, which drifted from the equatorial region.

One of the possible reasons for such catastrophes can be provoked by a periodic change in the direction of rotation of the planet. It can be caused by several celestial reasons, due to interplanetary perturbations or even dramatic collisions between planets. One can mention as an example Uranium, which, in contrast to the majority of planets, currently rotates around an axis lying almost in the plane of ecliptics.

We simulated such dramatic scenarios using periodic changes in the direction of the rotational axis. One of them is reproduced in Movie\_10.avi [19]. It shows a sequence of "geological" catastrophes, which are accompanied by huge surges of "geological activity". These catastrophes are separated by relatively long calm periods, during which the system normally relaxes to a small number of "continents", which form sometimes extremely realistic global patterns. Some of them strongly resemble the current pattern of our planet Earth, while others reproduce the pictures known from recorded geological history.
