A NUMERICAL STUDY ON POLAR RING GALAXY (PRG) FORMATION

The association between galaxies is a phenomenon frequently occurring in the universe. Systems called Polar Ring Galaxies (PRGs) can be built for certain interactions. These structures consist of a galaxy known as a host and a ring of gas and stars surrounding a plane nearly perpendicular to its main axis. This research has contributed to the creation of a first library of N-body simulations, to study the parameters required to construct PRGs. A library of this kind is continuously expanding with new simulations based on these first findings. The inherent secular processes of isolated galaxies cannot give rise to this PRG geometry. Such systems are (i) created by galaxy interactions or (ii) by the introduction of cold gas from interstellar filaments within a galaxy. The most popular mechanisms are thought to be (i), and we are interested in this project, particularly where gas from one galaxy to another is introduced, without the need for a fusion. For formation of interaction, orbits with particular geometries are required in order to attach matter to form a polar ring. We review various mass ratios, inclinations and orbital energies in our society. The database which will be developed with the simulation results will contribute to future research on the interaction between galaxies, especially polar rings.


Introduction
The study of interacting galaxies is extremely important for the understanding of the history of the Universe, as such encounters modify the cosmic structures in a very striking way throughout their evolution. These phenomena are relentlessly determined by the purely 2 attractive character of the gravitational force, which in turn induces the interaction of large systems, the tidal force, and dynamic friction 1 . The process of modelling a real system is done by choosing the appropriate equations and identifying the correct parameters that characterize the system under study well. In the case of galaxies, very complex systems, it is impossible to create equations that represent them perfectly 2 . Thus, it should be borne in mind that modelling is a representation of an ideal system, that is, it contains only some of the characteristics of the real system. Logically, the model can be made more precise, increasing its complexity. However, a good correspondence must be established between the simplicity of the model and the precision of the results obtained.
Alongside theoretical physics and experimental physics, computational physics has become an indispensable tool in the treatment of astrophysical problems 3 . The numerical simulation of gravitational systems of N-bodies, that is, systems of N particles interacting gravitationally, is one of the most powerful techniques for the study of astronomical systems such as star clusters, galaxies and structures in large scale 4 . This computational treatment is only possible due to the previous theoretical study through the mathematical modelling of the dynamics of these systems of interest. Galaxies, in clusters, are relatively close to each other 5 .
As such, they are in constant interaction with each other, including collisions. Numerical simulations of colliding galaxies can be performed looking for a study of their dynamics based on the morphology and measurements of stellar velocities and interstellar gas of the galaxies, obtained observationally. From the best adjustments between simulation and real data, it is possible to infer the age of the interaction, which can be compared with the age of the star populations, redoing the complex history of the system and predicting its future. This work aims to develop a library of N-body simulations to investigate the space of parameters necessary for the formation of Polar Ring Galaxies, which are systems made up of a host galaxy and a ring of gas and stars orbiting an almost perpendicular plane with in relation to your master plan. Due to the geometry of the system, its formation cannot be due to the intrinsic secular processes of isolated galaxies. Through observational data and numerical simulations made by several authors, cited during this work, it is believed that the most common mechanism of formation of these systems is that of interactions between two galaxies. The study of such mechanisms requires an exploration of different mass ratios, inclinations and orbital energies between two galaxies in order to form polar ring galaxies.
Such exploration is done in this work. According Contopoulos and Efthymiopoulos (2011) 6 , the movement of stars and dark matter elements is governed purely by the gravitational force. The study of these movements and their combinations to form self-consistent statistical mechanical configurations, has an extremely important role in the study of Stellar Dynamics, and constitutes the central approach to the study of Galaxy Dynamics.
1.2 Galaxies with Polar Ring: Polar Ring Galaxies (PRGs) are one of the most striking objects of the large family of unbarred ringed galaxies 7 . A PRG is formed by an early, lenticular or elliptical host galaxy, surrounded by a ring of gas and stars orbiting a plane almost perpendicular to the main plane of the host galaxy. There is a rare case of a spiral host galaxy 8 , however they are generally of the 0 type.
Although less than one percent of the 0 galaxies have polar rings, these characteristics are usually seen only at the edge, since they have very low surface brightness, so it is estimated that the real proportion of galaxies with these characteristics is higher 9 . A lower limit for the local universe is that 4.5% of galaxies with host characteristics have a polar ring 10 . It is conjectured that the rings are newer than the host galaxy itself because they have a higher concentration of interstellar gas. After the process of forming a PRG, when the ring has already stabilized, much of the "memory" of the formation process is lost, making it difficult to study its dynamic history. The observational results obtained for a PRG under study can be compared with the database simulations and, based on the analysis of the best match, its current and past dynamics can be investigated.
The high rotational speeds observed in galaxy discs suggest the presence of a halo of dark matter around it. The presence of two components (galaxy and polar ring) in perpendicular dynamic planes in systems such as PRGs allows for a more accurate study of the shape of these halos, thus being of great importance for understanding the behavior of dark matter.
More specifically, when comparing simultaneous measurements of rotation curves in the plane of the main galaxy and in the plane of the polar ring, information about the shape of the gravitational potential is obtained. Still, statistical comparisons between maximum rotational speeds in the polar rings and luminosities of the main galaxy point to a significant flattening in the dark matter halos 11 . Numerical simulations of systems like this, carried out under different initial conditions, help to understand their formation and their dynamic history, giving the possibility to understand structures previously observed as bizarre.    15 . There is also a compilation of GADGET-2 that uses another integration scheme, called Particle Mesh. This scheme was not used in this work, so it will not be detailed.
Stars and dark matter are modelled as self-gravitational non-collisional fluids, that is, their fdistribution functions in the phase space satisfy the non-collisional Boltzmann Equation  These particles are integrated along these characteristic curves. This is essentially a Monte Carlo method 16 , approximating continuous f for one sample. In self-consistent simulations, the gravitational field Φ is generated by N bodies. The bodies in a simulation model should not be confused with the stars of the real physical system, as each body in a simulation can represent a set of many stars in a real star system.

TreeCode -Tree Algorithm:
To simplify the numerical integration of the system of equation 1, the gravitational field is smoothed over a softening length. The expression of the gravitational potential for each particle ∈ {1, . . . , } is given by the Plummer Potential 17 : with ε being the so-called "softening parameter". In this way, we try to avoid obtaining extremely high values for the potential when two particles are at a distance | ri -r j | too small.
The simplest method of calculating the force exerted on each pair of particles is the obvious direct sum on each pair. Thus, the result would be, for a system of N particles, N (N -1) pairs, that is, an O (N 2 ) complexity algorithm would be obtained. A widely used technique to reduce computational cost in simulations is the TreeCode hierarchical tree algorithm, which explores the fact that the interaction of a particle with its neighbors is much more important than its interaction with distant individual particles, that is, the method takes into account the distance between the particles in the system. For the execution of the TreeCode, there are two distinct phases: the construction of the tree according to the structure of the system and the calculation of the force in each particle following the hierarchy defined by the composition of the tree.
In the process of building the tree, a cell that contains all the particles of the system is considered first. Such a cell is called a root cell. The intersection between cells is called a node. The cell is divided into 8 equal cubic cells. At this point, we count the number of particles in each cell. If there are no particles in the cell, the cell is ignored. If there is a particle in the cell, it is stored as a leaf node. If there is more than one particle in the cell, it is recorded as a branch node and we go back to the previous step (subdivide it until we have only one particle per cell). This process was initially proposed by Barnes and Hut 18 , and for this reason it will be denoted hereinafter by Árvore BH. Note that such a tree is an octal tree, given the division of each cell by 8 during the construction process. The complexity of the  19 , which represent a set of particles. This way, significant computing time is saved for systems with large N. To define proximity criteria between particles, a parameter θ, called a tolerance parameter, is used as follows: let d be the separation between a given particle and a group of particles contained in a cell of size s. If inequality ≤ is satisfied, then the internal particle distribution can be overlooked and the interaction computed using a low-order expansion of the potential of the particle group in relation to its center of mass. Thus, the force for each node is calculated following the algorithm below: • Step 1: Choose a particle as a root node.
• Step 2: Calculate the value of s for the node determined in the previous step.
• Step 3: If s ≤ θ, then the internal structure is calculated for the cell as a whole.
Otherwise, go up one level in the tree and go back to Step 1.

The SPH method and the TreeSPH implementation:
The evolution of interstellar gas is modelled numerically using the SPH (Smoothed Particle Hydrodynamics) technique in a way that is completely compatible with the tree structure of the TreeCode algorithm. The gas is treated as a smoothed fluid, following the Navier-Stokes equation for a compressible gas, with atoms and molecules exerting pressure and viscous forces on each other, which are not represented by the CBE. The gas density field is represented by particles that carry information capable of describing the local thermodynamic and hydrodynamic properties of the fluid, including an artificial viscosity that treats high shock densities and sources of entropy generation to account for heating and gas cooling. The dynamic equations are obtained from the Lagrangian form of the conservation laws for compressible fluids: where ρ is the density, v is the speed, P is the pressure, u is the thermal energy per unit mass and ℒ is the energy loss function, which includes all non-adiabatic energy sources and sinks.
Due to its nature, the SPH method is completely Lagrangian and naturally implementable in  Katz (1989) 20 , makes this unification. TreeSPH is considerably more flexible than other similar numerical schemes, being used by GADGET-2 in the treatment of three-dimensional self-gravitating fluids, with and without non-collisional material.
The purpose of the method is to describe the characteristics of the fluid by interpolating a set of particles. Interpolation is performed using a smoothing kernel W (smoothing kernel), which is the weighted average of the particles within a volume specified by a smoothing length h (smoothing length). W is normalized by W (r, h) dr = 1, and there are several techniques to choose from: all of them assume that there is no significant interaction between distant particles. Thus, a better performance of the calculations is achieved in a set of disordered points (the particles). Being (f (r)) the average value of one of the forces.
Interpolation allows any function to be expressed in terms of its field values f (r), the interpolation integral that allows finding such average value is given by ( ( )) = (| − ′|, ℎ) ( ′) ′. Of fundamental importance for any SPH formulation is density estimation. To do this, GADGET-2 uses the following equation: With each hi being the adaptive smoothing length of the i th particle.

Peano-Hilbert curve:
The Peano-Hilbert curve is a continuous space-filling fractal curve (they cross all points in the region to be covered), first described by the German mathematician David Hilbert (1862-1943) in 1891, as a variant of the curves of space filling space discovered by the mathematician Giuseppe Peano (1858-1932) 21 in the previous year.
The Peano-Hilbert curve is used in the implementation of GADGET-2 as a domain decomposition scheme, severely self-governing CPUs used. Such decomposition is done as follows: the total number of particles in the model is divided by the number of available processors, making all processors have the same amount of particles to work with. This prevents one processor from being overloaded while another is waiting.

Initial Conditions
Numerical tests performed by Bournaud and Combes (2015) 22 show that internal parameters of the donor and host galaxies (masses of discs, bowls and halos, as well as the rays of discs) are not the most relevant when compared to gaseous content (mass and radial extension) of 1) The relative speed before the encounter.
2) The angle between the plane of the donor disk and the equatorial plane of the host, before the material accretion event.
3) The angle between the donor's orbital plane in relation to the host, and the host's equatorial plane.
4) The minimum radius between the center of the host and the center of the donor.
The models generated by the MakeNewDisk code were written in the GADGET-17 format.
The program was modified to write the initial condition files in the GADGET-2 format.
These models can be rotated in a three-dimensional Cartesian reference system using the code uns_rotate.cc. Finally, the uns_stack.cc code causes two models in the new format to be joined in the same file of initial conditions, with the distance between the models defined, as well as the initial relative speed of the donor galaxy in relation to the host. The three codes were written in C ++ with the UNSIO API, and uns_rotate.cc was developed.

Simulations
First four models of galaxies were made: • Model A -Discoidal galaxy, based on the general parameters of the Milky Way. This model was designed to become the host galaxy for an eventual PRG.
• Model B -Galaxy with 50% of the mass of model A. More extended gas disk, which facilitates the capture of this material by the host to form the system ring.
• Model C -Similar to model B, but with 20% of the mass of model A.
• Model D -Similar to model B, but with 10% of the mass of model A.
The simulations grid consists of interacting model A with all the others, which sets the mass ratio between galaxies the first parameter to be tested. Thus, four groups of simulations are obtained: AA, AB, AC and AD. In a Cartesian system, the orbital plane in these four groups is the xy plane, and the pericenter of the orbit occurs over the positive section of the y axis. Model A is tilted 90 degrees with respect to the axis, with all models originally created with the plane on the  and as a donor (on the right). The green color represents the stars and the color bar represents the gas.
In all the figures referring to the simulations made, the reference legend located in the lower right corner must be understood as follows: the red axis is the axis, the green axis is the y and blue axis (with orientation coming out image) is the axis.
Two more parameters are tested in this grid: the pericenter distance q (for which the values of 8, 12 and 20 kpc were stipulated) and the initial relative speed v (three different values). The values of v were calculated for each group of simulations, the smallest being the one that leads to a zero energy orbit (which in the Keplerian case would be a parabola) that will be denoted by v0 (the escape velocity from the orbit of the donor galaxy with respect to the host), and the other two values v10 ≡ 1, 1v0 and v20 ≡ 1, 2v0. The grid therefore has 36 simulations. where X and Y are the interacting models, with pericenter distance n kpc and initial relative speed vj with j ∈ {0, 10, 20}.. Table 2 contains information about the formation or not of rings in each simulation of the grid.
The changes in speed in the initial conditions in each model of the grid mean changes in orbital energy: for speed v0, the orbital energy is zero, and the orbit is a perfect parabola. For strictly positive orbital energy (velocities greater than v0), the tidal force and dynamic friction can be expected to drag matter into the ring formation; for orbit with total negative energy, the orbit decay is faster and more violent, making it difficult for a possible ring to form.
The orbits are non-keplerian. When calculating the escape velocities at the pericenter for each simulation, it is possible to obtain the initial coordinates of velocity and position of the secondary galaxy with respect to the main. In this way, the initial distances between the two models in each pair of simulated galaxies were obtained. The speed at the beginning of each simulation was deduced, neglecting the effects caused by dynamic friction, through reverse temporal integration.

simAC-q8 simulations:
The simulations with initial relative speeds equal to v0 and v10 showed the formation of ring structures. The simAC-q8-v0 simulation showed, after approximately 2.4 Gyr, a very weak ring (Figure 3.2 (a)), which is quickly destroyed with the shock of the donor galaxy, resulting in a fusion of the two models, without the formation of rings. Here the ring formed is called "very weak" because it does not last long. However, the ring structure is well outlined. A weak ring structure is formed in the simAC-q8-v10 simulation. In this case, the transfer of gas from the donor to the host is more intense than in the previous case ( Figure 3.2 (b)), but also the ring is quickly destroyed with a shock between the two galaxies. The structure that forms after merging the models does not result in rings.

simAD-q8 simulations:
For the three speeds (v0, v10 and v20), this configuration of masses and distance from the pericentre presented a median ringed structure. In the three cases, the donor galaxy does not return after the initial shock for a new shock (and eventual fusion) until the 3 Gyr, which was the total simulated interaction time for all groups in the grid.

simAC-q12 simulations:
The simAC-q12-v0 simulation generated the formation of a weak ring, as shown in Figure 3.3 (b), but it quickly dissipates with the encounter between the donor and host galaxy, culminating in a ringless fusion.
In the simAC-q12-v10 simulation, the formation of a ring was also observed (Figure 3.3 (c)).
However, such a ring dissipates quickly, before an eventual merger between the two galaxies, being categorized as a very weak ring. The same event occurs in the simAC-q12-v20 simulation, with the formed ring illustrated in Figure 3.3 (d).

SimAD-q12 Simulations:
The simAD-q12-v0 simulation shows the formation of a ring that has a short duration, dissipating before an eventual fusion between the two galaxies, being another case of very weak ring formation. Figure 3.4 (a) shows the moment when the ring is formed. The same occurs with the simAD-q12-v10 ( Figure 3.4 (b)) and simAD-q12-v20 ( Figure 3.4 (c)) simulations.

simAB-q20 simulations:
In this group of simulations, the only one that presented a considerable ring formation was the one with speed v10. However, the ring that forms around the host galaxy ( Figure 3.4 (d)) soon dissipates with the approach of the donor galaxy. After 3 Gyr, the encounter between the donor and the host culminates in a fusion, without the formation of rings.

Computational consumption:
The simulations were carried out on Hipercubo, an IP&D cluster, which has both older nodes, with 12 colors each, and newer nodes, with 20 colors each. Each simulation was performed on a single node, that is, the communication between  Table 3 shows the average run time for each model of pair of galaxies, regardless of the pericenter distance and initial relative speed of each simulation.   Gyr. (d) simAB-q20-v10 simulation, after 1,7 Gyr.
The simulations shown in Table 3 would consume approximately 104 processing days if the simulations were all run serially on a 20-color node (which are the fastest). If the Hypercube was used only for the simulations of this work 1, the approximate time for using the cluster would be 20 days. As the Hypercube was not exclusively used for these simulations, the approximate time of use was 50 days, distributed over the last 4 months. Before that, several test simulations were performed (but not computed in the grid and, therefore, were not commented on in this work).  14 with the original grid, and in this way to refine the most promising pair models. In most cases that resulted in rings, the donor galaxy remains close to the host after the ring is formed, which is not common in observations. As the interactions were simulated up to 3 Gyr of temporal evolution, continuing the simulation, it is possible to verify the stability of these rings with the distance from the donor galaxy (if such distance occurs). The simAD-q12 group of simulations, in particular, presented results that proved to be possibly more conducive to the formation of more robust rings if the amount of gas from the donor was greater (gas transfer is important in the formation of rings, as it -star formation is observed in polar rings). Also, if the host slope at the beginning of the simulation were changed, new results could be obtained.
Thus, another galaxy model was created, called Model E, similar to model D, but with a greater fraction of total mass associated with gas. In this way, three new models of pairs of galaxies were built, with the configurations as explained below: • simAD-q12-v20-yrot90: The donor is model D and the host is model A, but inclined 90 degrees with respect to the y axis, thus having the plane of the disk parallel to the xz plane (Figure 3.5). The pericenter distance is 12 kpc, and the initial relative speed is the same used in the simAD-q12-v20 simulation.
• simAE-q12-v20: It has the same configurations as simAD-q12-v20, with model E in place of model D.
• Figure 3.5 -Initial conditions for simAD-q12-v20-yrot90 simulation. The host galaxy is on the left. The green color represents the stars and the color bar represents the gas. degrees with respect to the y axis. The small increase in the gaseous mass fraction made in the simAE-q12-v20 model showed results very close to those found in the simAD-q12-v20 simulation. New simulations, starting from the simAE-q12-v20 model, with a lower fraction of the gas mass of the host will be tested looking for a greater difference between the fractions of gas mass between the host and the donor. In addition to these, other models of pairs are being designed, aiming at new simulations with better results. The simAC-q12-v0 and simAC-q12-v10 models showed the most robust rings, but with relatively little gas.
Therefore, models of pairs containing the donor with the highest gas fraction starting from these configurations of mass, speed and distance from the pericenter will be made. Also, new configurations will be tested, such as different slopes of the host galaxy at the beginning of the temporal evolution (as was done in the models simAD-q12-v20-yrot90 and simAE-q12-v20-yrot90) and new distances of pericenter q: according to the results obtained in the original grid with q = 20 kpc, if the inclination of the host galaxy is 90 degrees with respect to the x-axis no rings occur. However, the results for 90-degree inclination with respect to the y-axis can generate quite different data, according to results obtained in the simAD-q12-v20-yrot90 and simAE-q12-v20-yrot90 simulations when compared with the simAD-q12 models.
-v20 and simAE-q12-v20, respectively. In these cases, the original models showed formation of annular structures, while the new ones did not. Something reverse can happen with models with q = 20 kpc. According to Bournaud and Combes (2003), it is possible for a PRG to be formed by the addition scenario with the host galaxy in the position as the host is in Figure   3.5. Regarding the computational consumption data of the simulations, the use of nodes with 20 colors proved to be much more advantageous. The increase in colors does not necessarily represent a linear decrease in process execution time, since the communication time between internal colors in the same node must be considered. However, in this case, the increase from 12 to 20 internal colors in the same node was interesting. However, not always the gradual increase in the amount of colors, when it comes to different nodes used in parallel by the same process, leads to a decrease in the execution time, because the communication time between the nodes (which is slower than the communication between colors in the same node) must be taken into account. The predicted new simulations will be carried out between combinations of different nodes. In this way, communication between nodes (both only between new nodes and only between old nodes, as well as between old and new nodes) can be tested. It is important to mention that these simulations were all carried out using the compilation of the GADGET-2 program for the use of N-bodies with the TreeCode and SPH  Table-3 can be obtained by compiling the program using the Particle Mesh scheme, which integrates differently. Tests with this scheme will be done in some new simulations, so that the execution times with the two schemes can be compared. In this way, a table combining the times required with each compilation and quantities of colors and nodes used can be made, aiming to provide statistical data that allow to further optimize the processing times of future simulations.

Conclusions
In this work, several numerical simulations were performed using Hipercubo. This work done by the researcher presented in this paper are useful for future simulations of N-bodies, helping to optimize the use of processors by researchers and students of the program.
Although the simulations performed did not result in the formation of robust polar rings, the parameter space for the addition scenario was restricted, facilitating the construction of new models of initial conditions for future simulations.
Due to the addition scenario, for the host tilted 90 degrees with respect to the x-axis, the formation of PRGs by interaction of galaxies of close masses seems to be less likely, given the low result of rings in models AA and AB (Table 2): of these, only the AB model with a pericenter distance equal to 20 kpc showed considerable results, which suggests that it may be possible to find formation of PRGs with this geometry and mass ratio only if the distance from the pericenter is slightly longer. The trend seems to be, for this geometry of the initial conditions, that polar rings are formed by the accretion scenario when the donor has a mass less than 50% of the host's mass, and the pericenter distance from the donor's orbit with respect to the host is less at 20 kpc. As many PRGs observed have the host galaxy being of type S0, new simulations with the host with little gas (or even without gas) are being made, lowering the gas mass fraction of model A and re-running the simulations that resulted in rings, even though they were weak. Few simulations with an initial host position other than that described above were discussed in this work. which will eventually contain even more data. Simulations for the galaxy models presented in this paper assembled with geometry conducive to the formation of PRGs by the fusion scenario will be carried out, although, this scenario is less likely to form PRGs than the accretion scenario.