2.1 Formulation of the statistical model

As an example we take the transitional disk DoAr44, which is a completely symmetrical dust ring around its star ^[7]. To model it we ideally consider that it is constituted by a spherical gas star and a toroid, made up of planetesimals with equal masses. Through an accretion process, planetesimals will give rise to massive objects (planets, dwarf planets, asteroids or comets).

Initially the particle distribution, R(r), must be obtained throughout 0 < r < ∞. For the proto-star the number of particles grows exponentially as r → 0, which explains why it finally collapses into a point object with respect to the initial volume of the gas. The particles in the toroid are around the ave- rage distance \r) with a positive asymmetry, that is, the number of particles is greater when r < 〈r〉 and less when r > 〈r〉. Statistically, this corresponds to the Gamma density function:

In the expression (1) you have the parameter α > λ called form of distribution (if a = 1 a star is obtained and if α > 1 a toroid) and a parameter λ > 0 called scale, which determines the scope of positive asymmetry from the right, as illustrated in Figure 2. The term Γ (α) is the Gamma function (for this reason, this name is given to the density function), in this case it is an extension of the factorial concept to real numbers and is responsible for the convergence of the distribution.

To determine the parameters a and X we make use of the concept: Complete Spatial Randomness (CSR). It is assumed that the particles are randomly found in space (this is called an event), therefore, from an average occurrence frequency of events λ, the probability that a certain number of events will occur (k) in a certain position r, obeys a Poisson distribution:

In the region of the considered space, the probability of the encounter between two particles is very small, that is, statistically it is a "rare event". Consequently, the probability is zero, if the encounter of more than two particles is considered, in a moment of time. The encounter of two particles in a region is independent of the encounter of other particles in any other disjunct region of the previous one. And the encounter of the particles occurs with regularity and independence, being the number of encounters theoretically unlimited. This behavior validates the correspondence with the Poisson distribution (2).

If the toroids give rise to the planets orbiting the star, we can assign a number n to each orbit. Therefore, the semimajor axis of the orbit can be ex- pressed as α = α _s n, where as is a free parameter (in astronomical units). It is not required that the orbit numbers n for the planets, be consecutive. As the dynamics is determined by the mass of the proto-star regardless of the mass of the object that the orbit, the parameter as it must be proportional to the mass of the proto-star, at the time when the toroid was formed. The average occurrence of events is the encounter of two particles at a distance a, i.e., λ = 2/ α = 2/ α _s n. If in the distance a, two are given events, it is expected that at the distance nα, 2n events occur; therefore, k = 2n.

The Poisson distribution (2), 0 ≤ P(k, r, λ) ≤ 1, by multiplying it by a factor λ it becomes the Gamma density function (1), where α = 2n + 1 and λ = 2/ α _s n. As now a it is a positive integer (α ε ℕ) the Gamma function Γ (α) = (α - 1)!, i.e., Γ (2n + 1) = 2n!. The Gamma distribution, when α is a positive integer, is known as the Erlang distribution. Consequently, the particle distribution for the toroid generated by a planet in the orbit n is,

since this expression is a density function, the area under the curve is the unit: The density function (3) is reduced to a Gaussian distribution for n > 9.

As seen in the Figure 3 the maximum points are in α = α _s n ² , consequently, it is expected that one or more massive objects (planets or dwarf planets) will be formed in this position and smaller objects farther away from this position (satellites, asteroids or comets).

In the Solar System, the planets are located at certain distances from the Sun. A succession has been found that relates the distance of the planets from the Sun known as the law of Titius-Bode r = 0,4 + 0,3 X 2ⁿ, with n =- ∞ 0, 1,... , 7 ^[8]. However, a better agreement is obtained considering the distance proportional to the square of a whole number: α = α _sn^{2 [9]}.

For the symmetry of the toroid the angular part is given by the spherical harmonics:

This provides a Gaussian density function with respect to θ and a rotational invariance O2. By integrating in the entire range of coordinates the unit is obtained:

For the star we have the distribution function:

where ∫0 ∞ Sn(r)dr = 1, and

like wise

Therefore, to describe a protoplanetary disk like Hl-Tauri, where it is evidenced that it is constituted by a series of concentric toroids [¹], with a density τ (r, θ, ϕ ). The star is made up of a series of concentric layers with density σ (r θ, ϕ):

where mn the mass of each of the Nτ toroids that make up the protoplanetary disk, and Mn is the mass of each of the Nσ layers that will give rise to the proto-star. Each of the mass terms depends on the distribution of matter, the absence of matter makes some of these terms zero, so that not all orbit numbers contribute in the summation. By integrating the expression (7) and (8) we will obtain the mass of the protoplanetary disk, ,and the mass of the star,

2.2. Simulation process

To simulate the protoplanetary disk Hl-Tauri, we start from the values of the major semi-axes of the bright rings: (αo = 20.4, 38.1, 47, 55, 68.8, 81.3, 97) AU [¹]. In the equation α = αn2 we do not know the free parameter as and the orbit numbers, n. As the value for the Solar System is α☉≃ 0-031 AU and it is considered that the Hl-Tauri star has a mass of Ms = 1-7 M☉ [¹⁰], we propose a series of possible values αs (0.04 ≤ ka☉ ≤ 0.06) AU, with a descending sweep of 10-5 AU, where k = Ms/M☉.

For each value, ak, the number of orbit is calculated from the rings, obtaining a series of values (r_o, α _k n ² ), and a least-squares adjustment is made to the equation, r _o = α _sn2, where a _s is the parameter to be determined. In this way, we obtain the value a _s =, where is the average and e _s is the dispersion (error). Additionally, the correlation coefficient 0 ≤ R ² ≤ 1 is obtained, and the chi-square value by degree of freedom and the goodness of fit (likelihood): 0 ≤ ≤ 1. The maximum likelihood is used as a criterion to select the best value of α _s . The results are shown in table 1.

For the transition disk DoAr44, M _s = 1.3 M☉, R = 1.75 Ra M, = 0.005 M☉ and the toroid is between 32 and 60 AU ^[7]. The same procedure as above is followed, using as criteria, to select the best value as, the lowest percentage error: e _s = 100 | α _s ,- α _o |/α _o The results obtained are: as = (3,94858± 0,00006)× 10⁻² AU, n = 33 and r _n = 43 AU, with these values the Figure 5 was made.

Table 1 shows the mass m _n of each of the rings ^[10]. The number of particles in each ring N _t is proportional to its percentage. The particles are randomly distributed in the area under the curve of the function R _n (r) in the range 0 ≤ r ≤ 130 AU, and each particle is assigned a random angle: 0 ≤ ϕ < 2π. A circular sector is obtained in the XY plane (See upper part of Figure 8), then each of these particles are randomly distributed in the area under the curve Yn(θ, 0) in the interval 0 < θ ≤ π, (See lower part of Figure 8). With these three values ??(r, θ, ϕ) the position of each particle in the Cartesian coordinate system is determined: x=rsinθ cosϕ, y=r sinθsinϕ, z=rcosθ. Figure 9 shows the radial function, the cross-sectional and zenithal view of the particle system obtained. In Figure 9 we can observe the absence of a ring in 28 AU, which is consistent with the collapse of the ring to planet.

The function (m _n /Md)R _n (r)Y _n (θ,ϕ) is computationally represented as e _n R _n (r) f _n Y _n (r) where en and f _n is a scale factor that corresponds to the maximum heights of each of the distributions; which, in turn, depend on the number of particles that are randomly distributed under the curve, as shown to the right of the Figure 5. The domain of the distribution function, R _n (r), is 0 < r < r _ext (AU), which we divide into Nr intervals: 1 ≤ i ≤ N, each division has a width, h _r = r _ext /N _r . The height of the i-th division is calculated by averaging the two heights of the interval R(i) = 1/2 (R _n [(i-1)h _r ] + R _n [ih _r ]). A random radius is generated between 0 < r ≤ r _ext and it is located in the i-th division to which it belongs, as long as the points contained in the division are less than R(i), in such a way that filling is completed under the curve of Figure 5 with a certain number of points. The points are located at an angle 0 ≤ ϕ ≤ 2π random, obtaining the toroid on the XY plane. The axial symmetry is given by the ϕ angle that we already had in mind when projecting the toroid in the XY plane. This distribution function can be represented with a system of particles distributed in space as illustrated in Figure 5.

The domain of the angular distribution function Y _n (θ, 0) is 0 < θ ≤ π, which we divide into N _θ = N _r π/r _ext intervals: 1 ≤ i ≤ N _θ , each division has a width, h _θ = π /N _θ . The height of the i-th division is calculated by: Y (i) = 1/2 (Y _n [(i - 1)h _θ ] + Y _n [ih _θ ]). A random angle is generated between 0 < θ ≤ π and it is located in the i-th division to which it belongs, as long as the points contained in the division are less than Y (i), in such a way that filling is completed under the curve of Figure 5 with a certain number of points. The points are located along the z-axis from its XY position obtained previously with the radial distribution.

Finally, we have the radius 0 < r ≤ r _ext , the angles 0 ≤ θ < π and 0 ≤ Ø < 2 π, with which the position of the particles in the XYZ Cartesian coordinate system is determined:

The magnitude of the velocity of the particles, v _o , is determined from the opposing forces: the gravitational force against the centripetal force (-m d Φ (r)/dr|_r=α = sin θ - mv ² /α sin θ = 0):

where Φ(r) it is the potential under which the particles are subjected. The velocities in space are obtained by deriving the positions of the particles, equation (11), they are:

To simulate the collision between two particles, each one is considered as a sphere whose radius is half the minimum distance between the particles of the distribution. This was done with the purpose of inducing a first collision to accelerate the simulation. Initially, it is examined that there is no overlap between the spokes and the directions of the speeds must be opposite for the collision to exist; if these conditions are met, the second particle is suppressed, decreasing the number of particles of the simulation in one, this causes the speed of the simulation to increase depending on the frequency of the collisions. The resulting particle complies with the law of conservation of the mass, volume and momentum of the particles that collided. The Hl-Tauri disk consists of a series of toroids with a mass that depends on the number of number of particles in each of the radial distributions. The masses of each toroid were obtained from ^[7] m _n = 47, 53, 30, 62, 7, 101, 123 M⊕, which affect the height of the distribution function. The radial distribution function and the corresponding particle distribution are illustrated in the Figure 6. Unlike other models, the formation of planetary bodies is faster and generates stable and closed orbits ^[11].