© ESO, 2014

1. Introduction

The importance of stationary solutions to the Boltzmann or Fokker-Planck equations for stellar objects or other systems is that they provide a starting point for dynamical evolutions and a way to study, analytically or numerically, the stability of the corresponding system. Furthermore, from the distribution function of a system, we can obtain its thermodynamic variables by allowing a full understanding of the physics of each individual model. Many such studies have been conducted over the years (see e.g. Cohn 1979; Watanabe et al. 1981; Nishida et al. 1984; Yepes & Domínguez-Tenreiro 1992; Takahashi 1995; Theuns 1996; Takahashi et al. 1997; Joshi et al. 2001; Ashurov 2004; Fiestas & Spurzem 2010). Owing to the nature of the problem, the equations are usually solved numerically with approximative methods such as Monte Carlo (statistics) and moment equations (averages). On the other hand, direct numerical methods are not frequently considered, because of their computational complexity in two and three dimensions (Ujevic & Letelier 2005, 2006). Analytical solutions of the collisionless Boltzmann equation for thin discs were given in Mestel (1963) and Kalnajs (1972).

Regardless the importance of the Kalnajs and Mestel thin discs, we have to keep in mind that they do not model a realistic physical situation because they exhibit odd behaviours within the disc. Namely, in the Kalnajs discs, the distribution function goes to infinity for certain values of the energy density and the angular momentum; while in the Mestel discs, the gravitational potential Φ diverges at the centre of the disc. Thus, the problem of finding self-consistent stellar models for galaxies starting from known well-behaved gravitational potential models is of broad astrophysical interest. However, the procedure for finding reasonable physical models (either analytical or numerical) of the complete Fokker-Planck equation is not a simple task.

The main goal of our work is to use a direct numerical method to find self-consistent planar ring models that satisfy the Fokker-Planck equation and the Poisson equation and to study some of their properties. In many astrophysical systems, light is emitted by a thin stellar disc, thus our study of planar galaxies or planar systems is more than academically interesting. From the variety of gravitational ring models that have been proposed in the literature, we are going to focus on families of analytic ring models of Vogt & Letelier (2009) and the analytic homogeneous ring solution of Lass & Blitzer (1983). In the formulation of the Fokker-Planck equation, we include the collision term of the local approximation and calculate the diffusion coefficients as prescribed by Rosenbluth et al. (1957). Furthermore, based on the mass density profile of each particular model, a simple boundary condition for the Fokker-Planck equation has been successfully applied. To our knowledge, there has been no analytical or numerical solution of the Fokker-Planck equation for flat systems with rings.

The article is organized as follows. Section 2 introduces the Fokker-Planck equation for a thin disc with the collision term given by the local approximation and the diffusion coefficients. In Sect. 3, the computational code used to solve the Fokker-Planck equation is briefly explained. Section 4, presents the properties of the flattened ring models used in this article and the numerical results obtained. Finally, in Sect. 5, our results are summarized.

2. Fokker-Planck equation for flat systems with rings

To obtain a physical self-consistent model for an astrophysical system, its distribution function f must satisfy the Fokker-Planck equation and the Poisson equation simultaneously in a given time t, say, $\begin{array}{ccc}& & \frac{\mathit{\partial f}}{\mathit{\partial t}}\mathrm{+}{v}\mathrm{·}\mathrm{\nabla }\mathit{f}\mathrm{-}\mathrm{\nabla \Phi }\mathrm{·}{\mathrm{\nabla }}_{\mathit{v}}\mathit{f}\mathrm{=}\mathrm{\Gamma }\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathit{,}\\ & & {\mathrm{\nabla }}^{\mathrm{2}}\mathrm{\Phi }\mathrm{=}\mathrm{4}\mathit{\pi Gm}\mathrm{\int }\mathit{f}{\mathrm{d}}^{\mathrm{3}}\mathit{v,}\end{array}$where v represents the velocity of the particles (e.g. stars), ∇ is the usual gradient, ∇_v is the gradient in the velocity space (derivatives with respect to the velocities), G the gravitational constant, Φ is the gravitational potential, m the mass of each particle, and the symbol Γ [f] represents the Fokker-Planck collision term. The collision term in the local approximation has the form $\mathrm{\Gamma }\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{v}}^{\mathit{i}}}\mathrm{[}\mathit{fD}{\mathrm{\left(}\mathrm{\Delta }{\mathit{v}}^{\mathit{i}}\mathrm{\right)}}^{\mathrm{]}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i,j}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{v}}^{\mathit{i}}\mathit{\partial }{\mathit{v}}^{\mathit{j}}}\mathrm{[}\mathit{fD}{\mathrm{\left(}\mathrm{\Delta }{\mathit{v}}^{\mathit{i}}\mathrm{\Delta }{\mathit{v}}^{\mathit{j}}\mathrm{\right)}}^{\mathrm{]}}\mathit{,}$(3)where the functions D(Δv_i) and D(Δv_iΔv_j) are known as the diffusion coefficients.

These coefficients were calculated (Rosenbluth et al. 1957) in the three-dimensional case by assuming an inverse square force for the interaction, and each stellar encounter involves only a single pair of stars and is independent of all others. These coefficients can be simplified if the field-star distribution function is described by a Maxwellian distribution (Binney & Tremaine 2008). We recall that the diffusion coefficients are calculated considering a test star of mass m moving through an infinite homogeneous sea of field stars of mass m_a who has mean velocity equal to zero. The use of these coefficient in real actual systems, which are not homogeneous neither infinite, therefore have to be seen only as a first approximation. The reason for this simplification is that in general the diffusion coefficients also depend on the distribution function, so the Fokker-Planck equation is actually a more complicated integro-differential equation. Also, the diffusion coefficients are calculated by assuming encounters between the elements of the system and not collisions. (We use the word “encounter” to denote the gravitational perturbation of the orbit of one element of the system by another element, and the word “collision” to denote physical contact between elements. However, we keep “collision term” for the quantity Γ [f ].) Therefore, the results obtained in the following sections may not be valid in systems where collisions are important. Furthermore, the collision term is valid in the local approximation, which means that the encounter affects only the velocity and not the position of the interacting stars. For this reason strong deflections might be excluded.

Regarding the Fokker-Planck equation, it is important to mention that the individual stellar encounters gradually perturb stars away from their trajectories such that after many encounters, a star eventually loses track of its initial trajectory. The characteristic time interval in which these encounters are relevant to the evolution of the system corresponds to the relaxation time. The Fokker-Planck equation is appropriate for timescales exceeding the relaxation time, while the collisionless Boltzmann equation is appropriate for timescales lower than the relaxation time. It is known that encounters may have played an important part in the structure of stellar systems, such as globular clusters, galactic nuclei, open clusters nuclei, and clusters of galaxies.

Equations (1) and (2), together with the collision term (3) are valid for a general three-dimensional case. As a first approximation in our study we consider a case of a very thin disc with thickness approximately equal to δz ≈ 10^-8% of the size of its radius and stars with low velocities in the z direction. By taking these approximations into account, we can rewrite the Fokker-Planck equation in cylindrical coordinates as $\frac{\mathit{\partial f}}{\mathit{\partial t}}\mathrm{+}\mathit{ṙ}\frac{\mathit{\partial f}}{\mathit{\partial r}}\mathrm{+}\mathit{\varphi ̇}\frac{\mathit{\partial f}}{\mathit{\partial \varphi }}\mathrm{+}\left(\frac{{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{2}}}{\mathit{r}}\mathrm{-}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial r}}\right)\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{v}}_{\mathit{r}}}\mathrm{-}\left(\frac{{\mathit{v}}_{\mathit{r}}{\mathit{v}}_{\mathit{\varphi }}}{\mathit{r}}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial \varphi }}\right)\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{v}}_{\mathit{\varphi }}}\mathrm{=}\mathrm{\Gamma }\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathit{,}$(4)with $\mathrm{\Gamma }\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{v}}^{\mathit{i}}}\mathrm{[}\mathit{fD}{\mathrm{\left(}\mathrm{\Delta }{\mathit{v}}^{\mathit{i}}\mathrm{\right)}}^{\mathrm{]}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i,j}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{v}}^{\mathit{i}}\mathit{\partial }{\mathit{v}}^{\mathit{j}}}\mathrm{[}\mathit{fD}{\mathrm{\left(}\mathrm{\Delta }{\mathit{v}}^{\mathit{i}}\mathrm{\Delta }{\mathit{v}}^{\mathit{j}}\mathrm{\right)}}^{\mathrm{]}}\mathit{.}$(5)The same result can be obtain if we start with Eq. (1) and perform a projection on the z = 0 and v_z = 0 planes using Dirac deltas, i.e. doing the replacement f(x,v,t) → f(x,v,t)δ(z)δ(v_z) and integrating in the z and v_z variables (Ujevic & Letelier 2005). The explicit form of the collision term Γ [f] is presented in Appendix A.

3. Numerical method

The general procedure for finding the self-consistent distribution function is the following. First, we fix our system in a given time t, then we insert the gravitational potential (Φ) of our system into the Fokker-Planck equation. After that, we set the values of the physical parameters for the model and the boundary conditions of the system. Then, we solve the Fokker-Planck equation numerically to obtain the distribution function (f). Once we determine the distribution function, we calculate the mass density numerically with the relation $\mathrm{\Sigma }\mathrm{(}\mathit{r}{\mathrm{)}}_{\mathrm{Num}}\mathrm{=}\mathit{m}\mathrm{\int }\mathit{f}\mathrm{(}\mathit{r,}{\mathit{v}}_{\mathit{r}}\mathit{,}{\mathit{v}}_{\mathit{\varphi }}\mathrm{)}\mathrm{d}{\mathit{v}}_{\mathit{r}}\mathrm{d}{\mathit{v}}_{\mathit{\varphi }}\mathit{.}$(6)To simplify the integration in Eq. (6) we set the upper boundary velocity in each domain point to be the maximum escape velocity of the whole system, i.e. $\mathit{v}\mathrm{=}\sqrt{{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathit{\varphi }}^{\mathrm{2}}}\mathrm{\le }{\mathit{v}}_{\mathrm{esc}}$, where ${\mathit{v}}_{\mathrm{esc}}\mathrm{=}\sqrt{\mathrm{2}\mathrm{|}\mathrm{\Phi }{\mathrm{|}}_{\max }}$. Finally, we compare the numerical mass density obtained by Eq. (6) with the original analytical mass density, which is associated with the gravitational potential inserted in the Fokker-Planck equation. If both densities coincide, then we conclude that the distribution function is self-consistent, and if they do not coincide, we modify the values of the physical parameters and repeat the procedure until a self-consistent distribution function is reached. When the desired distribution function is obtained, we can calculate the average value of any observable as a function of r by using $\mathrm{⟨}\mathrm{𝒪}\mathrm{⟩}\mathrm{=}\mathrm{\int }\mathrm{𝒪}\mathit{f}\mathrm{\left(}\mathit{r,}{\mathit{v}}_{\mathit{r}}\mathit{,}{\mathit{v}}_{\mathit{\varphi }}\mathrm{\right)}\mathrm{d}{\mathit{v}}_{\mathit{r}}\mathrm{d}{\mathit{v}}_{\mathit{\varphi }}/\mathrm{\int }\mathit{f}\mathrm{\left(}\mathit{r,}{\mathit{v}}_{\mathit{r}}\mathit{,}{\mathit{v}}_{\mathit{\varphi }}\mathrm{\right)}\mathrm{d}{\mathit{v}}_{\mathit{r}}\mathrm{d}{\mathit{v}}_{\mathit{\varphi }}\mathit{,}$(7)where the integrals (6) and (7) are calculated numerically using the Simpson 1/3 method.

We discretized the Fokker-Planck equation by applying a finite differencing method, which uses a central-differences second-order scheme for the first and the second-order partial derivatives. In this work, the gravitational potentials are axially symmetric, therefore the distribution function f depends on three variables: one spatial variable (r) and two velocity variables (v_r and v_ϕ). We therefore used an eleven-point discretization molecule. (The mixed derivatives of Γ [f ] in Eq. (5) introduce a more complex discretization molecule.) We used a grid of 51 nodes for the r variable and 25 nodes for each velocity variable.

After the discretization, the resulting system of linear algebraic equations is solved using the robust left conjugate direction method (Catabriga et al. 2004). The boundary condition for our problem must have a spatial part and a velocity part. An appropriate boundary condition for the velocity part can be found if we recall that the velocity distribution in a star system is determined by two competing processes: there is the relaxation through the stellar encounters that drives the distribution toward a Maxwellian form, and on the other hand, there are stars that, by going beyond the finite escape velocity, continually disappear from the system. Since the escaping stars traverse the system in a small fraction of a relaxation time, the resulting velocity distribution should be zero beyond the escape velocity, but otherwise nearly Maxwellian. Thus, we propose a Maxwellian distribution for the velocities at the border. Now, for the spatial boundary condition we found that the best results are obtained when we set the spatial part equal to the analytic mass surface density of the particular model evaluated at the border. Therefore, we set the boundary condition at the border of the system as ${\mathit{f}}_{\mathrm{Boun}}\mathrm{\left(}\mathit{r,v}\mathrm{\right)}\mathrm{\propto }\exp \mathrm{[}\mathrm{-}{\mathit{v}}^{\mathrm{2}}\mathit{/}\mathrm{2}{{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}^{\mathrm{]}}\mathrm{·}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(8)where Σ(r) is the analytical mass surface density of the particular model, and σ_v is a constant with velocity dimension.

With this computational setup we have obtained self-consistent solutions when the collision term is of the same order of magnitude (or bigger) as the larger term on the left-hand-side of Eq. (1). For lower values of the collision term, we start having problems finding self-consistent models. One computational problem is that when we decrease the values of the collision terms the discretization matrix (representing the system of linear equations) that we need to solve has low values in the main diagonal, and most of the computational methods fail to give us a correct (or accurate) solution for this case. Perhaps changing the discretization method or the numerical method used to invert the matrix will solve this problem. Also, we may have to consider more realistic diffusion coeficients for our systems. At the moment we are studying different alternatives that will be discussed elsewhere.

4. Flattened galactic ring model

4.1. Flat Kuzmin-Toomre family of rings

Recently, Vogt & Letelier (2009) have found a family of flat rings system (characterized by zero density on r = 0) and a family of a flat discs with rings using the superposition of the members of the Kuzmin-Toomre family of discs. All the models have axial symmetry and infinite extension, thus, to make the numerical calculations, we limit the extension of the systems by setting a cutoff radius (r = r_max). It is interesting that this family of solutions also admits a relativistic representation (Ujevic et al. 2011).

Fig. 1 Results for the ring (1,1). We consider the parameters a = 4 kpc, σ = 48.87 km s-1, ℳ(1,1) = 109 M⊙, ρ = 1.35 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 4 kpc, δz = 6 × 10-10 kpc, σv = 0.8σ, and rmax = 4a. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface profile obtained from the numerical integration of Eq. (6) (full circles). In the region with the highest density, the difference between them is less than 0.8%. Graph b) shows the values of ⟨ v2 ⟩ /σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for case r/a ≈ 0.72, respectively. Graphs e) and f) show the spatial distribution function and the level curves for case vϕ = 0, respectively.

4.1.1. Family of one-ring systems

The family of one-ring systems is defined by the mass surface density ${\mathrm{\Sigma }}^{\mathrm{\left(}\mathit{i,j}\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{(}\mathit{r}\mathit{/}\mathit{a}{\mathrm{)}}^{\mathrm{2}\mathit{i}}}{{\left(\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\right)}^{\mathit{i}\mathrm{+}\mathit{j}\mathrm{+}\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}$(9)where Σ₀ is a constant with dimension of mass surface density, a is a parameter of the model, i = 1,2,... and j = 0,1,.... The total mass ℳ^(i,j) can be obtained by simply integrating Σ^(i,j) over the space. The gravitational potentials associated to the density (9) can be calculated with the same coefficients as used in the superposition of (9) (see Vogt & Letelier (2009) for details). For example, if we take the (1,1) model (this corresponds to i = 1 and j = 1), the mass density of the model is ${\mathrm{\Sigma }}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{\right)}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{(}\mathit{r}\mathit{/}\mathit{a}{\mathrm{)}}^{\mathrm{2}}}{{\left(\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\right)}^{\mathrm{7}\mathit{/}\mathrm{2}}}\mathit{,}$(10)and the respective gravitational potential is (in z = 0) ${\mathrm{\Phi }}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{\right)}}\mathrm{(}\mathit{r,z}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}\mathrm{-}\frac{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{0}}\mathit{G}{\mathit{a}}^{\mathrm{2}}}{\mathrm{15}{\left({\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\right)}^{\mathrm{5}\mathit{/}\mathrm{2}}}\mathrm{(}\mathrm{2}{\mathit{r}}^{\mathrm{4}}\mathrm{+}\mathrm{7}{\mathit{r}}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{+}\mathrm{2}{\mathit{a}}^{\mathrm{4}}\mathrm{)}\mathit{.}$(11)Now, with (11) in the Eq. (4) we can find the distribution function of the system, and then the numerical mass surface density using Eq. (6). For the outer boundary condition, we set Eq. (8) and for the inner boundary condition we set f = 0.

Figure 1 shows the numerical results for the single ring (1,1) with parameters a = 4 kpc, σ = 48.87 km s^-1, ℳ^(1,1) = 10⁹M_⊙, ρ = 1.35 × 10^-13 kg m^-3, m_a = M_⊙, v_typ = 0.8σ, b_m = 4 kpc, δz = 6 × 10^-10 kpc, σ_v = 0.8σ, and r_max = 4a. In Fig. 1a the mass surface density obtained after integrating (6) is compared with the exact mass density (10). The agreement between them is very good. Namely, in the densest part of the graph, near r/a ≈ 0.7, the biggest difference appears between them into account, which is less than 0.8%. If we take all the assumptions made so far and the fact that we are trying to match a real finite physical system to an infinite geometric potential that was developed without using any observational data, then the result is more than satisfactory.

With the help of the distribution function and Eq. (7) in Fig. 1b, the values of ⟨ v² ⟩/σ² within the ring are plotted, and this quantity is related with the temperature of the system. We see that the highest value of the temperature coincides with the densest region, the temperature drops to zero near the origin, and it is almost constant in the rest of the system. To get more complete information about the system, we present in Fig. 1c, d the velocity profile and the level curves of the distribution function for the value r/a ≈ 0.72, respectively, and in Fig. 1e, f the spatial profile and the level curves of the distribution function for the value v_ϕ = 0, respectively. The level curves of Fig. 1d move to the left for values of r/a< 0.72, and for values of r/a> 0.72, the level curves move to the right. In Fig. 1f the level curves are the same for different values of v_ϕ, changing only the magnitude of the distribution function. We performed several calculations of the single-ring model (1,1) with different values of ℳ^(1,1) and a, finding self-consistent models for different physical parameters, given that we have properly chosen the values of σ and ρ. For instance, for v_typ = 0.8σ, b_m ≈ a, r_max = 4a, and σ_v = 0.8σ, the self-consistency occurs when the parameters ℳ^(1,1), m_a, a, and δz satisfy the inequality ln [ 0.71(ℳ^(1,1)/m_a) ] ·(m_a/ ℳ^(1,1))·(a/δz) ≳ 102.19. After many calculations we have concluded that this inequality is only a limitation of the numerical method employed and not a physical limitation.

Besides the single-ring model (1,1) presented here, we also investigated different single-ring models (i,j) of Eq. (9) with good results.

Fig. 2 Results for the double ring (1, 2). We consider the parameters a = 4 kpc, σ = 24.80 km s-1, ${\mathrm{ℳ}}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, κ = 1.2, ρ = 0.12 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 4 kpc, δz = 3 × 10-9 kpc, σv = 0.8σ, and rmax = 4a. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface density profile obtained from the numerical integration of Eq. (6) (full circles). In the region with the highest density the difference is less than 1.7%. Graph b) shows the values of ⟨ v2 ⟩/σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for the case r/a ≈ 2.32, respectively. Graphs e) and f) show the spatial distribution function and the level curves for the case vϕ = 0, respectively.

4.1.2. Family of double-ring systems

In this case, the family of double-ring systems is defined by the mass surface density ${\mathrm{\Sigma }}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathit{i,j}\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{(}\mathit{r}\mathit{/}\mathit{a}{\mathrm{)}}^{\mathrm{\left(}\mathrm{2}\mathit{i}\mathrm{\right)}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\kappa }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}}{{\left(\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\right)}^{\mathit{i}\mathrm{+}\mathit{j}\mathrm{+}\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}$(12)where κ is a dimensionless nonzero constant. In these models a gap is placed at r/a = 1/κ, thus generating a family of two concentric flat rings. As a study case we take the (1, 2) model. For this model, the mass surface density is ${\mathrm{\Sigma }}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}\frac{\mathrm{(}\mathit{r}\mathit{/}\mathit{a}{\mathrm{)}}^{\mathrm{2}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\kappa }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}}{{\left(\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\right)}^{\mathrm{9}\mathit{/}\mathrm{2}}}\mathit{,}$(13)and the respective gravitational potential is (in z = 0) $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{(}\mathit{r,z}\mathrm{=}\mathrm{0}\mathrm{)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{0}}\mathit{G}{\mathit{a}}^{\mathrm{2}}}{\mathrm{105}{\mathit{\chi }}^{\mathrm{7}\mathit{/}\mathrm{2}}}\{\mathrm{6}{\mathit{\chi }}^{\mathrm{3}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{8}{\mathit{\kappa }}^{\mathrm{4}}\mathrm{)}\\ & & \mathrm{+}\hspace{0.17em}{\mathit{\chi }}^{\mathrm{2}}\mathrm{[}\mathrm{-}\mathrm{16}{\mathit{\kappa }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}{\mathrm{(}\mathrm{9}\mathrm{-}\mathrm{26}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{-}\mathrm{54}{\mathit{\kappa }}^{\mathrm{4}}\mathrm{)}}^{\mathrm{]}}\\ & & \mathrm{-}\hspace{0.17em}\mathrm{2}{\mathit{a}}^{\mathrm{2}}\mathit{\chi }\mathrm{[}{\mathit{r}}^{\mathrm{2}}{\mathrm{(}\mathrm{2}\mathrm{+}\mathrm{11}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{+}\mathrm{9}{\mathit{\kappa }}^{\mathrm{4}}\mathrm{)}}^{}\\ & & \mathrm{-}\hspace{0.17em}{\mathit{a}}^{\mathrm{2}}{\mathrm{(}\mathrm{7}\mathrm{+}\mathrm{28}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{+}\mathrm{21}{\mathit{\kappa }}^{\mathrm{4}}\mathrm{)}}^{\mathrm{]}}\\ & & \mathrm{-}\hspace{0.17em}\mathrm{3}{\mathit{a}}^{\mathrm{4}}{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}\mathrm{[}\mathrm{2}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathrm{7}{\mathit{a}}^{\mathrm{2}}\mathrm{]}\}\mathit{,}\end{array}$(14)where χ = r² + a². For this family of rings we consider the same boundary conditions used in the single-ring systems.

In Fig. 2 we present the numerical results for the double-ring model (1,2) with parameters a = 4 kpc, σ = 24.80 km s^-1, ${\mathrm{ℳ}}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, κ = 1.2, ρ = 0.12 × 10^-13 kg m^-3, m_a = M_⊙, v_typ = 0.8σ, b_m = 4 kpc, δz = 3 × 10^-9 kpc, σ_v = 0.8σ, and r_max = 4a. In Fig. 2a the mass surface density obtained after integrating (6) is compared with the exact mass surface density (13). The biggest difference between them is less than 1.7%, and it occurs in the densest region. This is again a good result for the same reasons as in the single-ring model. In Fig. 2b the quantity ⟨ v² ⟩/σ², which corresponds to the temperature of the system, is plotted.

The temperature drops to zero in the regions without particles, i.e. at the origin and at the intersection of the rings. The highest temperature is located in the inner ring, and in the outer ring the temperature is almost constant. For completeness, we present in Fig. 2c, d the velocity profile and the level curves of the distribution function for the value r/a ≈ 2.32, respectively. For values of r/a > 2.32, the level curves of Fig. 2d are nearly the same, for values of 0.37 < r/a < 0.80 the level curves move to the right and for other values of r the level curves move to the left. Moreover, in Fig. 2e, f we plot the spatial profile and the level curves of the distribution function for the value v_ϕ = 0, respectively. The level curves of Fig. 2f have the same profile for different values of v_ϕ, varying the magnitude of the distribution function.

We performed several calculations for the double-ring model (1, 2) with different values of ${\mathrm{ℳ}}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}$, a, and κ and find self-consistent models, given again (as in the single-ring case) that we have properly chosen the values of σ and ρ. As an example for the double – ring model, for v_typ = 0.8σ, b_m ≈ a, r_max = 4a, σ_v = 0.8σ, and κ = 1.2, the self-consistency occurs when the parameters ${\mathrm{ℳ}}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}$, m_a, a, and δz satisfy the inequality $\ln \mathrm{\left[}\mathrm{0.18}\mathrm{\right(}{\mathrm{ℳ}}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathit{/}{\mathit{m}}_{\mathit{a}}\mathrm{\left)}\mathrm{\right]}\mathrm{·}\mathrm{(}{\mathit{m}}_{\mathit{a}}\mathit{/}{\mathrm{ℳ}}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{)}\mathrm{·}\mathrm{(}\mathit{a}\mathit{/}\mathit{\delta z}\mathrm{)}\gtrsim \mathrm{19}$. For this case and as in the previous section, the conclusion is that this inequality is only a limitation of the numerical method employed and not a physical limitation.

Besides the double-ring model (1, 2) discussed above, we went on to investigate different double- ring models (i,j) of Eq. (12) with similarly good results.

Fig. 3 Results for the disc with a ring for the i = 2 model. We consider the parameters a = 4 kpc, σ = 29.74 km s-1, ${\mathrm{ℳ}}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, κ = 1.5, ρ = 0.21 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 4 kpc, δz = 2 × 10-9 kpc, σv = σ, and rmax = 4a. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface density profile obtained from the numerical integration of Eq. (6) (full circles). In the region with the highest density, the difference is less than 2%. Graph b) shows the values of ⟨ v2 ⟩/σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for the case r/a ≈ 0.24, respectively. Graphs e) and f) show the spatial distribution function and the level curves for the case vr = 0, respectively.

Fig. 4 Results for the homogeneous ring model. We consider the following parameters a = 10 kpc, b = 12 kpc, σ = 36.95 km s-1, ℳR = 109M⊙, ρ = 0.24 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 2 kpc, δz = 2 × 10-8 kpc, and σv = 0.9σ. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface density profile obtained from the numerical integration of Eq. (6) (full circles). The highest difference between them is less than 2.2%. Graph b) shows the values of ⟨ v2 ⟩ /σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for the case r/a ≈ 1.09, respectively. Graphs e) and f) show the spatial distribution function and the level curves for the case vr = 0, respectively.

4.1.3. Family of discs with flat rings

The family of discs with flat rings is defined by the mass surface density ${\mathrm{\Sigma }}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathit{i}\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\kappa }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}}{{\left(\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\right)}^{\mathit{i}\mathrm{+}\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}$(15)where κ is a dimensionless nonzero constant. This surface mass density represents a disc with radius r/a = 1/κ and a flat ring between 1/κ<r/a < ∞. As a study case, we take the i = 2 model. For this model, the mass surface density is ${\mathrm{\Sigma }}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\frac{{\mathrm{\Sigma }}_{\mathrm{0}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\kappa }}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{2}}}{{\left(\mathrm{1}\mathrm{+}{\mathit{r}}^{\mathrm{2}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\right)}^{\mathrm{7}\mathit{/}\mathrm{2}}}\mathit{,}$(16)and the respective gravitational potential is (in z = 0) $\begin{array}{ccc}{\mathrm{\Phi }}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{(}\mathit{r,z}\mathrm{=}\mathrm{0}\mathrm{)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{2}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{0}}\mathit{G}{\mathit{a}}^{\mathrm{2}}}{\mathrm{15}{\mathit{\chi }}^{\mathrm{5}\mathit{/}\mathrm{2}}}\{{\mathit{\chi }}^{\mathrm{2}}\mathrm{(}\mathrm{3}\mathrm{-}\mathrm{4}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathit{\chi }\mathrm{[}\mathrm{3}{\mathit{a}}^{\mathrm{2}}\mathrm{-}\mathrm{4}{\mathit{\kappa }}^{\mathrm{2}}{\mathit{a}}^{\mathrm{2}}\mathrm{+}\mathrm{8}{\mathit{\kappa }}^{\mathrm{4}}{\mathit{r}}^{\mathrm{2}}\mathrm{]}\\ & & \mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{[}\mathrm{-}{\mathit{r}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{a}}^{\mathrm{2}}{\mathrm{(}\mathrm{2}\mathrm{+}\mathrm{4}{\mathit{\kappa }}^{\mathrm{2}}\mathrm{+}\mathrm{3}{\mathit{\kappa }}^{\mathrm{4}}\mathrm{)}}^{\mathrm{]}}\}\mathit{,}\end{array}$(17)where again χ = r² + a².

In Fig. 3 we present the numerical results for the model i = 2 of a disc with a ring for the parameters a = 4 kpc, σ = 29.74 km s^-1, ${\mathrm{ℳ}}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, κ = 1.5, ρ = 0.21 × 10^-13 kg m^-3, m_a = M_⊙, v_typ = 0.8σ, b_m = 4 kpc, δz = 2 × 10^-9 kpc, σ_v = σ, and r_max = 4a. In Fig. 3a, the exact mass surface density (16) is compared with the numerical integrated mass density (6). We find good agreement between them, and the biggest difference is again in the densest region (less than 2.0%). In Fig. 3b we plot the quantity ⟨ v² ⟩/σ² to measure the temperature of the system. The temperature drops to zero in the intersection of the disc with the ring where there are no particles. The highest temperature is in the central disc, and the temperature is almost constant in the ring. Figure 3c, d present the velocity profile and the level curves of the distribution function for the value r/a ≈ 0.24, respectively, and in Fig. 3e, f the spatial profile and the level curves of the distribution function for the value v_r = 0 is plotted, respectively. In Fig. 3d, the level curves move to the right for values of r/a between 0.24 < r/a < 0.63 and to the left for other values of r. For different values of v_r, we obtain the same profile of Fig. 3f, changing only their magnitudes. We performed several calculations of the i = 2 model with different values of ${\mathrm{ℳ}}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$, a, and κ, again finding self-consistent models for properly chosen values of σ and ρ. This time our example shows that for v_typ = 0.8σ, b_m ≈ a, r_max = 4a, σ_v = σ, and κ = 1.5, the self-consistency occurs when the parameters ${\mathrm{ℳ}}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$, m_a, a, and δz satisfy the inequality $\ln \mathrm{\left[}\mathrm{0.26}\mathrm{\right(}{\mathrm{ℳ}}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{/}{\mathit{m}}_{\mathit{a}}\mathrm{\left)}\mathrm{\right]}\mathrm{·}\mathrm{(}{\mathit{m}}_{\mathit{a}}\mathit{/}{\mathrm{ℳ}}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{)}\mathrm{·}\mathrm{(}\mathit{a}\mathit{/}\mathit{\delta z}\mathrm{)}\gtrsim \mathrm{27}$. As in previous sections, this inequality should only be a limitation of the numerical method. Furthermore, good results were obtained for different i models of Eq. (15).

4.2. Homogeneous ring

In this section we consider an analytic homogeneous ring model written in terms of complete elliptic integrals (Lass & Blitzer 1983). In this model, the ring has inner radius a, outer radius b, and a constant density Σ₀, while the mass M_R generates the gravitational potential (written in cylindrical coordinates) $\begin{array}{ccc}\mathrm{\Phi }\mathrm{(}\mathit{r,z}\mathrm{=}\mathrm{0}\mathrm{)}& \mathrm{=}& \mathrm{2}\mathit{G}{\mathrm{\Sigma }}_{\mathrm{0}}[\mathrm{(}\mathit{a}\mathrm{+}\mathit{r}\mathrm{)}\mathit{E}\mathrm{\left(}{\mathit{\eta }}_{\mathit{a}}\mathrm{\right)}\mathrm{+}\mathrm{(}\mathit{a}\mathrm{-}\mathit{r}\mathrm{)}\mathit{K}\mathrm{\left(}{\mathit{\eta }}_{\mathit{a}}\mathrm{\right)}\\ & & \mathrm{-}\hspace{0.17em}\mathrm{(}\mathit{b}\mathrm{+}\mathit{r}\mathrm{)}\mathit{E}\mathrm{\left(}{\mathit{\eta }}_{\mathit{b}}\mathrm{\right)}\mathrm{-}\mathrm{(}\mathit{b}\mathrm{-}\mathit{r}\mathrm{)}\mathit{K}\mathrm{\left(}{\mathit{\eta }}_{\mathit{b}}\mathrm{\right)}]\mathit{,}\end{array}$(18)where η_a = 4ar/ (a + r)², η_b = 4br/ (b + r)², and K and E are complete elliptic integrals of the first and of the second order, respectively. The gravitational acceleration can be written as $\begin{array}{ccc}\mathrm{-}\frac{\mathit{\partial }{\mathrm{\Phi }}_{\mathit{ba}}}{\mathit{\partial r}}& \mathrm{=}& \mathrm{2}\mathit{G}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{[}\mathrm{(}\mathit{r}\mathrm{+}\mathit{b}\mathrm{)}\mathit{E}\mathrm{\left(}{\mathit{\eta }}_{\mathit{b}}\mathrm{\right)}\mathrm{􏼮}\mathit{r}\mathrm{-}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{b}}^{\mathrm{2}}\mathrm{)}}^{}\mathit{K}\mathrm{\left(}{\mathit{\eta }}_{\mathit{b}}\mathrm{\right)}{\mathrm{􏼮}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathit{b}{\mathit{r}}^{\mathrm{)}}}^{}}^{}\\ & & \mathrm{-}\hspace{0.17em}\mathrm{(}\mathit{r}\mathrm{+}\mathit{a}\mathrm{)}\mathit{E}\mathrm{\left(}{\mathit{\eta }}_{\mathit{a}}\mathrm{\right)}\mathrm{􏼮}\mathit{r}\mathrm{+}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}{\mathit{a}}^{\mathrm{2}}\mathrm{)}}^{}\mathit{K}\mathrm{\left(}{\mathit{\eta }}_{\mathit{a}}\mathrm{\right)}{\mathrm{􏼮}{\mathrm{(}{\mathit{r}}^{\mathrm{2}}\mathrm{+}\mathit{a}{\mathit{r}}^{\mathrm{)}}}^{}}^{\mathrm{]}}\mathit{,}\end{array}$(19)where we used the relations $\frac{\mathit{\partial K}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}}{\mathit{\partial \eta }}\mathrm{=}\frac{\mathit{E}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }\mathrm{)}\mathit{K}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}}{\mathrm{2}\mathit{\eta }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\eta }\mathrm{)}}\mathit{,} \frac{\mathit{\partial E}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}}{\mathit{\partial \eta }}\mathrm{=}\frac{\mathit{E}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{-}\mathit{K}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}}{\mathrm{2}\mathit{\eta }}\mathrm{·}$(20)The gravitational acceleration has divergences at the boundary of the ring, i.e. in r = a and r = b, because the function K(η) diverges when η = 1. To incorporate the radial acceleration in the Fokker-Planck equation is necessary to calculate all the values of the function E(η) and K(η). These functions will be evaluated using Romberg’s integration scheme with a Richardson type of extrapolation. For the inner ring and for the outer ring boundary condition, we again consider condition (8).

In Fig. 4 we present the numerical results for a homogeneous ring with parameters a = 10 kpc, b = 12 kpc, σ = 36.95 km s^-1, ℳ_R = 10⁹M_⊙, ρ = 0.24 × 10^-13 kg m^-3, m_a = M_⊙, v_typ = 0.8σ, b_m = 2 kpc, δz = 2 × 10^-8 kpc, and σ_v = 0.9σ. In Fig. 4a the exact mass surface density (a constant) is compared with the numerical integrated mass density. We see that in the proximities of the borders of the ring, the calculated numerical density is less accurate than in the rest of the ring, and they show oscillatory behaviour. The error at the borders is less than 2.2%, thus the numerical integrated mass is in good agreement with the exact one considering that the gravitational acceleration has discontinuities at the borders. In Fig. 4b we plot the quantity ⟨ v² ⟩/σ². The temperature is almost constant within the ring, with some fluctuations in the vicinity of the borders possibly due to the boundary condition. Figure 4c,d shows the velocity profile and the level curves of the distribution function for the value r/a ≈ 1.09, respectively. For different values of r and v_r, Fig. 4d and f have the same profiles, respectively. Finally, in Fig. 4e, f, we plot the spatial profile and the level curves of the distribution function for the value v_r = 0.

Besides the model presented in Fig. 4 we performed several calculations with different values of ℳ_R, a, and b, and found self-consistent models for proper values of σ and ρ. For example, for v_typ = 0.8σ, b_m ≈ a, b/a = 1.2, and σ_v = 0.9σ, the self-consistency occurs when the parameters ℳ_R, m_a, a, and δz satisfy the inequality ln [0.20(ℳ_R/m_a)] ·(m_a/ℳ_R)·(a/δz) ≳ 2. As in previous sections, this inequality is only attributed to the limitation of the numerical method.

5. Conclusion

To obtain initial self-consistent data for the Fokker-Planck equation is, in general, not an easy task. In this work we found several solutions for systems of flat rings that satisfy the Fokker-Planck equation and the Poisson equation simultaneously. The collision term of the Fokker-Planck equation was obtained by taking the local approximation into account, while the diffusion coefficients were calculated following the work of Rosenbluth et al. (1957), where only gravitational encounters are taken into account. In particular, we found the distribution function of the flattened ring models of Vogt & Letelier (2009) and the flat homogeneous ring model of Lass & Blitzer (1983). It is interesting that we managed to model physical systems starting from geometric potentials that were developed without any observational data. To our knowledge, these are the first numerical solutions of the Fokker-Planck equation for flat systems with rings. Furthermore, it appears that the simple boundary condition used in our work gives good results for a large set of physical parameters. It will be very important to test this condition in different two-dimensional models and three-dimensional structures. Our results provide the basis needed for further evolution studies of the Fokker-Planck equation when different kinds of perturbations are applied.

Acknowledgments

C.A. thanks the CAPES for financial support, R.C. thanks the UFABC and CAPES for financial support, and M.U. thanks the CNPq for financial support.

References

Appendix A: Explicit determination of the collision term

The collision term for the two-dimensional FP equation in the local approximation is given by Eq. (5), where the diffusion coefficients were determined by Rosenbluth in 1957 and have the following form $\begin{array}{ccc}& & \mathit{D}\mathrm{\left[}\mathrm{\Delta }{\mathit{v}}_{\mathit{i}}\mathrm{\right]}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{G}}^{\mathrm{2}}{\mathit{m}}_{\mathit{a}}\mathrm{(}\mathit{m}\mathrm{+}{\mathit{m}}_{\mathit{a}}\mathrm{)}\ln \mathrm{\Lambda }\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{v}}_{\mathit{i}}}\mathrm{\int }{\mathit{d}}^{\mathrm{3}}{{v}}_{\mathit{a}}\frac{{\mathit{f}}_{\mathit{a}}\mathrm{\left(}{{v}}_{\mathit{a}}\mathrm{\right)}}{\mathrm{|}{v}\mathrm{-}{{v}}_{\mathit{a}}\mathrm{|}}\mathit{,}\\ & & \mathit{D}\mathrm{\left[}\mathrm{\Delta }{\mathit{v}}_{\mathit{i}}\mathrm{\Delta }{\mathit{v}}_{\mathit{j}}\mathrm{\right]}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{G}}^{\mathrm{2}}{\mathit{m}}_{\mathit{a}}^{\mathrm{2}}\ln \mathrm{\Lambda }\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{v}}_{\mathit{i}}\mathit{\partial }{\mathit{v}}_{\mathit{j}}}\mathrm{\int }{\mathit{d}}^{\mathrm{3}}{{v}}_{\mathit{a}}{\mathit{f}}_{\mathit{a}}\mathrm{\left(}{{v}}_{\mathit{a}}\mathrm{\right)}\mathrm{|}{v}\mathrm{-}{{v}}_{\mathit{a}}\mathrm{|}\mathit{,}\end{array}$(A.1)where m and m_a are the mass of the test particle and the field particle, respectively, $\mathrm{\Lambda }\mathrm{=}\frac{{\mathit{b}}_{\mathit{m}}{\mathit{v}}_{\mathrm{typ}}^{\mathrm{2}}}{\mathit{G}\mathrm{(}\mathit{m}\mathrm{+}{\mathit{m}}_{\mathit{a}}\mathrm{)}}$, v_typ is the typical velocity of the particles in the system, and b_m is the maximum impact parameter. Numerical experiments show that the appropriate value for b_m is approximately the orbital radius (Theuns 1996).

The diffusion coefficients can be evaluated explicitly when the distribution function of the field particles is known. The most important case is when the distribution function of field particles is Maxwellian, i.e. ${\mathit{f}}_{\mathit{a}}\mathrm{\left(}{\mathit{v}}_{\mathit{a}}\mathrm{\right)}\mathrm{=}\frac{\mathit{n}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{e}}^{\mathrm{-}{\mathit{v}}_{\mathit{a}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\mathit{,}$(A.2)where n is the number density of particles, and σ the one-dimensional velocity dispersion of the field particles. For this particular case, it is possible to find explicit expressions for the diffusion coefficients $\mathit{D}\mathrm{\left[}\mathrm{\Delta }{\mathit{v}}_{\mathit{i}}\mathrm{\right]}\mathrm{=}\mathit{\gamma \sigma }\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\right]}\mathit{,} \mathit{D}\mathrm{\left[}\mathrm{\Delta }{\mathit{v}}_{\mathit{i}}\mathrm{\Delta }{\mathit{v}}_{\mathit{j}}\mathrm{\right]}\mathrm{=}\mathit{\gamma }{\mathit{\sigma }}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{j}}\mathrm{\right]}\mathit{,}$(A.3)where $\mathit{\gamma }\mathrm{=}\frac{\mathrm{4}\mathit{\pi }{\mathit{G}}^{\mathrm{2}}\mathit{\rho }{\mathit{m}}_{\mathit{a}}\ln \mathrm{\Lambda }}{{\mathit{\sigma }}^{\mathrm{3}}}$, ρ = nm_a, ϵ_i = v_i/σ. We define the new dimensionless diffusion coefficients $\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\right]}\mathrm{=}\mathrm{-}\mathrm{𝒜}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}{\mathit{ϵ}}_{\mathit{i}}\mathit{,} \begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{j}}\mathrm{\right]}\mathrm{=}\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }\mathrm{ℬ}}{\mathit{\partial ϵ}}{\mathit{ϵ}}_{\mathit{i}}{\mathit{ϵ}}_{\mathit{j}}\mathrm{+}\mathrm{ℬ}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}{\mathit{\delta }}_{\mathit{ij}}\mathit{,}$(A.4)where $\begin{array}{ccc}& & \mathrm{𝒜}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\mathit{\mu }\frac{\mathrm{Erf}\left(\frac{\mathit{ϵ}}{\sqrt{\mathrm{2}}}\right)\mathrm{-}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathit{ϵ}{\mathrm{e}}^{\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{3}}}\mathit{,}\\ & & \mathrm{ℬ}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}{\mathit{ϵ}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{Erf}\left(\frac{\mathit{ϵ}}{\sqrt{\mathrm{2}}}\right)\mathrm{+}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathit{ϵ}{\mathrm{e}}^{\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{3}}}\mathit{,}\\ & & \mathit{\mu }\mathrm{=}\frac{\mathit{m}\mathrm{+}{\mathit{m}}_{\mathit{a}}}{{\mathit{m}}_{\mathit{a}}}\mathit{,} \mathit{ϵ}\mathrm{=}\sqrt{{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{ϵ}}_{\mathrm{2}}^{\mathrm{2}}}\mathit{.}\end{array}$(A.5)Replacing the expressions (A.3) in the collision term (5), we find$\mathrm{\Gamma }\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathit{/}\mathit{\gamma }\mathrm{=}\mathrm{-}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{ϵ}}_{\mathit{i}}}\{\mathit{f}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\right]}\}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sum _{\mathit{i,j}\mathrm{=}\mathrm{1}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}}{\mathit{\partial }{\mathit{ϵ}}_{\mathit{i}}\mathit{\partial }{\mathit{ϵ}}_{\mathit{j}}}\{\mathit{f}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{j}}\mathrm{\right]}\}\mathit{.}$(A.6)Evaluating the derivatives we can write the previous equation as $\mathrm{\Gamma }\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathit{/}\mathit{\gamma }\mathrm{=}{\mathit{C}}_{\mathrm{0}}\mathit{f}\mathrm{+}{\mathit{C}}_{\mathrm{1}}\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{ϵ}}_{\mathrm{1}}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}\frac{\mathit{\partial f}}{\mathit{\partial }{\mathit{ϵ}}_{\mathrm{2}}}\mathrm{+}{\mathit{C}}_{\mathrm{11}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{+}{\mathit{C}}_{\mathrm{22}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{ϵ}}_{\mathrm{2}}^{\mathrm{2}}}\mathrm{+}{\mathit{C}}_{\mathrm{12}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{f}}{\mathit{\partial }{\mathit{ϵ}}_{\mathrm{1}}\mathit{\partial }{\mathit{ϵ}}_{\mathrm{2}}}\mathit{,}$(A.7)where $\begin{array}{ccc}{\mathit{C}}_{\mathrm{0}}& \mathrm{=}& \hspace{0.17em}\frac{\mathrm{1}}{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}{\mathrm{]}}_{\mathit{,}\mathrm{11}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}{\mathrm{]}}_{\mathit{,}\mathrm{22}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}{\mathrm{]}}_{\mathit{,}\mathrm{12}}\\ & & \mathrm{-}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}{\mathrm{]}}_{\mathit{,}\mathrm{1}}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}{\mathrm{]}}_{\mathit{,}\mathrm{2}}\mathit{,}\\ {\mathit{C}}_{\mathrm{1}}& \mathrm{=}& \hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}{\mathrm{]}}_{\mathit{,}\mathrm{1}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}{\mathrm{]}}_{\mathit{,}\mathrm{2}}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}\mathrm{\right]}\mathit{,}\\ {\mathit{C}}_{\mathrm{2}}& \mathrm{=}& \hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}{\mathrm{]}}_{\mathit{,}\mathrm{2}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{1}}{\mathrm{]}}_{\mathit{,}\mathrm{1}}\mathrm{-}\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathrm{2}}\mathrm{\right]}\mathit{,}\\ {\mathit{C}}_{\mathrm{11}}& \mathrm{=}& \end{array}$(A.8)Finally, the derivatives of the dimensionless diffusion coefficients are determined from the following relationships: $\begin{array}{ccc}& & \begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}{\mathrm{]}}_{\mathit{,i}}\mathrm{\equiv }\frac{\mathit{\partial }\begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{\left[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\right]}}{\mathit{\partial }{\mathit{ϵ}}_{\mathit{i}}}\mathrm{=}\mathrm{-}\mathrm{𝒜}\mathrm{-}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }\mathrm{𝒜}}{\mathit{\partial ϵ}}\mathit{,}\\ & & \begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{j}}{\mathrm{]}}_{\mathit{,i}}\mathrm{=}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}{\mathit{ϵ}}_{\mathit{j}}{\left[\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }}{\mathit{\partial ϵ}}\right]}^{\mathrm{2}}\mathrm{ℬ}\mathrm{+}\mathrm{(}{\mathit{ϵ}}_{\mathit{j}}\mathrm{+}\mathrm{2}{\mathit{\delta }}_{\mathit{ij}}{\mathit{ϵ}}_{\mathit{i}}\mathrm{)}\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }\mathrm{ℬ}}{\mathit{\partial ϵ}}\mathit{,}\\ & & \begin{array}{c}\mathrm{˜}\\ \mathit{D}\end{array}\mathrm{[}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{i}}\mathrm{\Delta }{\mathit{ϵ}}_{\mathit{j}}{\mathrm{]}}_{\mathit{,ij}}\mathrm{=}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}{\mathit{ϵ}}_{\mathit{j}}^{\mathrm{2}}{\left[\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }}{\mathit{\partial ϵ}}\right]}^{\mathrm{3}}\mathrm{ℬ}\mathrm{+}\mathrm{(}{\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{ϵ}}_{\mathit{j}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{\delta }}_{\mathit{ij}}{\mathit{ϵ}}_{\mathit{i}}{\mathit{ϵ}}_{\mathit{j}}\mathrm{)}{\left[\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }}{\mathit{\partial ϵ}}\right]}^{\mathrm{2}}\mathrm{ℬ}\\ & & \end{array}$(A.9)

where $\begin{array}{ccc}& & \frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }\mathrm{𝒜}}{\mathit{\partial ϵ}}\mathrm{=}\mathit{\mu }\frac{\mathrm{-}\mathrm{3}\mathrm{Erf}\left(\frac{\mathit{ϵ}}{\sqrt{\mathrm{2}}}\right)\mathrm{+}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathit{ϵ}\mathrm{(}{\mathit{ϵ}}^{\mathrm{2}}\mathrm{+}\mathrm{3}\mathrm{)}{\mathrm{e}}^{\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{5}}}\mathit{,}\\ & & \frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }\mathrm{ℬ}}{\mathit{\partial ϵ}}\mathrm{=}\mathrm{-}\frac{\mathrm{(}{\mathit{ϵ}}^{\mathrm{2}}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{Erf}\left(\frac{\mathit{ϵ}}{\sqrt{\mathrm{2}}}\right)\mathrm{+}\mathrm{3}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathit{ϵ}{\mathrm{e}}^{\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{5}}}\mathit{,}\\ & & {\left[\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }}{\mathit{\partial ϵ}}\right]}^{\mathrm{2}}\mathrm{ℬ}\mathrm{=}\frac{\mathrm{3}\mathrm{(}{\mathit{ϵ}}^{\mathrm{2}}\mathrm{-}\mathrm{5}\mathrm{)}\mathrm{Erf}\left(\frac{\mathit{ϵ}}{\sqrt{\mathrm{2}}}\right)\mathrm{+}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathit{ϵ}\mathrm{(}\mathrm{2}{\mathit{ϵ}}^{\mathrm{2}}\mathrm{+}\mathrm{15}\mathrm{)}{\mathrm{e}}^{\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{7}}}\mathit{,}\\ & & {\left[\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }}{\mathit{\partial ϵ}}\right]}^{\mathrm{3}}\mathrm{ℬ}\mathrm{=}\mathrm{-}\frac{\mathrm{15}\mathrm{(}{\mathit{ϵ}}^{\mathrm{2}}\mathrm{-}\mathrm{7}\mathrm{)}\mathrm{Erf}\left(\frac{\mathit{ϵ}}{\sqrt{\mathrm{2}}}\right)\mathrm{+}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}\mathit{ϵ}\mathrm{(}\mathrm{2}{\mathit{ϵ}}^{\mathrm{4}}\mathrm{+}\mathrm{20}{\mathit{ϵ}}^{\mathrm{2}}\mathrm{+}\mathrm{105}\mathrm{)}{\mathrm{e}}^{\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}\mathit{/}\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{9}}}\mathrm{·}\\ & & \end{array}$(A.10)In the case where ϵ → 0, we replace in Eq. (A.9) the following limits $\begin{array}{ccc}& & \underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\lim }\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }\mathrm{𝒜}}{\mathit{\partial ϵ}}\mathrm{\to }\mathrm{-}\mathit{\mu }\frac{\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}}{\mathrm{5}}\mathit{,} \underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\lim }\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }\mathrm{ℬ}}{\mathit{\partial ϵ}}\mathrm{\to }\mathrm{-}\frac{\mathrm{2}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}}{\mathrm{15}}\mathit{,}\\ & & \underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\lim }{\left[\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }}{\mathit{\partial ϵ}}\right]}^{\mathrm{2}}\mathrm{ℬ}\mathrm{\to }\frac{\mathrm{2}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}}{\mathrm{35}}\mathit{,}\underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\lim }{\left[\frac{\mathrm{1}}{\mathit{ϵ}}\frac{\mathit{\partial }}{\mathit{\partial ϵ}}\right]}^{\mathrm{3}}\mathrm{ℬ}\mathrm{\to }\mathrm{-}\frac{\mathrm{2}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}}{\mathrm{63}}\mathit{,}\\ & & \underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\lim }\mathrm{𝒜}\mathrm{\to }\mathit{\mu }\frac{\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}}{\mathrm{3}}\mathit{,} \underset{\mathit{ϵ}\mathrm{\to }\mathrm{0}}{\lim }\mathrm{ℬ}\mathrm{\to }\frac{\mathrm{2}\sqrt{\frac{\mathrm{2}}{\mathit{\pi }}}}{\mathrm{3}}\mathit{,}\end{array}$(A.11)in order to avoid numerical instabilities.

All Figures

Fig. 1 Results for the ring (1,1). We consider the parameters a = 4 kpc, σ = 48.87 km s-1, ℳ(1,1) = 109 M⊙, ρ = 1.35 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 4 kpc, δz = 6 × 10-10 kpc, σv = 0.8σ, and rmax = 4a. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface profile obtained from the numerical integration of Eq. (6) (full circles). In the region with the highest density, the difference between them is less than 0.8%. Graph b) shows the values of ⟨ v2 ⟩ /σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for case r/a ≈ 0.72, respectively. Graphs e) and f) show the spatial distribution function and the level curves for case vϕ = 0, respectively.

In the text

Fig. 2 Results for the double ring (1, 2). We consider the parameters a = 4 kpc, σ = 24.80 km s-1, ${\mathrm{ℳ}}_{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, κ = 1.2, ρ = 0.12 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 4 kpc, δz = 3 × 10-9 kpc, σv = 0.8σ, and rmax = 4a. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface density profile obtained from the numerical integration of Eq. (6) (full circles). In the region with the highest density the difference is less than 1.7%. Graph b) shows the values of ⟨ v2 ⟩/σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for the case r/a ≈ 2.32, respectively. Graphs e) and f) show the spatial distribution function and the level curves for the case vϕ = 0, respectively.

In the text

Fig. 3 Results for the disc with a ring for the i = 2 model. We consider the parameters a = 4 kpc, σ = 29.74 km s-1, ${\mathrm{ℳ}}_{\mathrm{\left(}\mathit{d}\mathrm{\right)}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}{\mathrm{10}}^{\mathrm{9}}{\mathit{M}}_{\mathrm{\odot }}$, κ = 1.5, ρ = 0.21 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 4 kpc, δz = 2 × 10-9 kpc, σv = σ, and rmax = 4a. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface density profile obtained from the numerical integration of Eq. (6) (full circles). In the region with the highest density, the difference is less than 2%. Graph b) shows the values of ⟨ v2 ⟩/σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for the case r/a ≈ 0.24, respectively. Graphs e) and f) show the spatial distribution function and the level curves for the case vr = 0, respectively.

In the text

Fig. 4 Results for the homogeneous ring model. We consider the following parameters a = 10 kpc, b = 12 kpc, σ = 36.95 km s-1, ℳR = 109M⊙, ρ = 0.24 × 10-13 kg m-3, ma = M⊙, vtyp = 0.8σ, bm = 2 kpc, δz = 2 × 10-8 kpc, and σv = 0.9σ. In graph a) the mass surface density profile of the exact solution (solid line) is compared with the mass surface density profile obtained from the numerical integration of Eq. (6) (full circles). The highest difference between them is less than 2.2%. Graph b) shows the values of ⟨ v2 ⟩ /σ2 within the ring. Graphs c) and d) show the velocity distribution function and the level curves for the case r/a ≈ 1.09, respectively. Graphs e) and f) show the spatial distribution function and the level curves for the case vr = 0, respectively.

In the text