Numerical selfconsistent distribution function of flattened ring models
Centro de Ciências Naturais e Humanas, Universidade Federal do
ABC,
09210170
Santo André,
São Paulo,
Brazil
email: cesar.alarcon@ufabc.edu.br; ricardo.calister@ufabc.edu.br; mujevic@ufabc.edu.br
Received:
14
January
2014
Accepted:
31
March
2014
We provide numerical, selfconsistent distribution functions for several flat ring models by simultaneously solving the FokkerPlanck equation and the Poisson equation. In particular, we calculated the distribution function of flat ring systems formed by superposing KuzminToomre disc solutions and an analytic homogeneous ring solution in terms of complete elliptic integrals. We used these geometrical disc solutions, together with physical parameters, to model more realistic physical systems. Moreover, we defined a cutoff radius to handle the infinite KuzminToomre disc families numerically. The FokkerPlanck equation is solved by a direct numerical method using the finite difference method, the left conjugate direction algorithm, and a simple boundary condition that allows us to find good results for a large set of physical parameters and for values of the collision term of the same order of magnitude (or larger) when compared with the other terms in the FokkerPlanck equation. The collision term of the FokkerPlanck equation is explicitly calculated by applying the known Rosenbluth potentials for gravitational encounters. Limitations of the method are also discussed.
Key words: galaxies: kinematics and dynamics / methods: numerical
© ESO, 2014
1. Introduction
The importance of stationary solutions to the Boltzmann or FokkerPlanck 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ínguezTenreiro 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 selfconsistent stellar models for galaxies starting from known wellbehaved gravitational potential models is of broad astrophysical interest. However, the procedure for finding reasonable physical models (either analytical or numerical) of the complete FokkerPlanck equation is not a simple task.
The main goal of our work is to use a direct numerical method to find selfconsistent planar ring models that satisfy the FokkerPlanck 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 FokkerPlanck 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 FokkerPlanck equation has been successfully applied. To our knowledge, there has been no analytical or numerical solution of the FokkerPlanck equation for flat systems with rings.
The article is organized as follows. Section 2 introduces the FokkerPlanck 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 FokkerPlanck 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. FokkerPlanck equation for flat systems with rings
To obtain a physical selfconsistent model for an astrophysical system, its distribution function f must satisfy the FokkerPlanck equation and the Poisson equation simultaneously in a given time t, say, 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 FokkerPlanck collision term. The collision term in the local approximation has the form (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 threedimensional 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 fieldstar 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 FokkerPlanck equation is actually a more complicated integrodifferential 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 FokkerPlanck 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 FokkerPlanck 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 threedimensional 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 FokkerPlanck equation in cylindrical coordinates as (4)with (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 selfconsistent 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 FokkerPlanck 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 FokkerPlanck equation numerically to obtain the distribution function (f). Once we determine the distribution function, we calculate the mass density numerically with the relation (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. , where . 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 FokkerPlanck equation. If both densities coincide, then we conclude that the distribution function is selfconsistent, and if they do not coincide, we modify the values of the physical parameters and repeat the procedure until a selfconsistent 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 (7)where the integrals (6) and (7) are calculated numerically using the Simpson 1/3 method.
We discretized the FokkerPlanck equation by applying a finite differencing method, which uses a centraldifferences secondorder scheme for the first and the secondorder 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 elevenpoint 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 (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 selfconsistent solutions when the collision term is of the same order of magnitude (or bigger) as the larger term on the lefthandside of Eq. (1). For lower values of the collision term, we start having problems finding selfconsistent 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 KuzminToomre 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 KuzminToomre 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)} = 10^{9} 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 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 ⟨ v^{2} ⟩ /σ^{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. 

Open with DEXTER 
4.1.1. Family of onering systems
The family of onering systems is defined by the mass surface density (9)where Σ_{0} 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 (10)and the respective gravitational potential is (in z = 0) (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^{9}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^{2} ⟩/σ^{2} 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 singlering model (1,1) with different values of ℳ^{(1,1)} and a, finding selfconsistent 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 selfconsistency 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 singlering model (1,1) presented here, we also investigated different singlering 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}, , κ = 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 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 ⟨ v^{2} ⟩/σ^{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. 

Open with DEXTER 
4.1.2. Family of doublering systems
In this case, the family of doublering systems is defined by the mass surface density (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 (13)and the respective gravitational potential is (in z = 0) (14)where χ = r^{2} + a^{2}. For this family of rings we consider the same boundary conditions used in the singlering systems.
In Fig. 2 we present the numerical results for the doublering model (1,2) with parameters a = 4 kpc, σ = 24.80 km s^{1}, , κ = 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 singlering model. In Fig. 2b the quantity ⟨ v^{2} ⟩/σ^{2}, 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 doublering model (1, 2) with different values of , a, and κ and find selfconsistent models, given again (as in the singlering 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 selfconsistency occurs when the parameters , m_{a}, a, and δz satisfy the inequality . 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 doublering 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}, , κ = 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 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 ⟨ v^{2} ⟩/σ^{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 v_{r} = 0, respectively. 

Open with DEXTER 
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} = 10^{9}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 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 ⟨ v^{2} ⟩ /σ^{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 v_{r} = 0, respectively. 

Open with DEXTER 
4.1.3. Family of discs with flat rings
The family of discs with flat rings is defined by the mass surface density (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 (16)and the respective gravitational potential is (in z = 0) (17)where again χ = r^{2} + a^{2}.
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}, , κ = 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^{2} ⟩/σ^{2} 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 , a, and κ, again finding selfconsistent 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 selfconsistency occurs when the parameters , m_{a}, a, and δz satisfy the inequality . 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 Σ_{0}, while the mass M_{R} generates the gravitational potential (written in cylindrical coordinates) (18)where η_{a} = 4ar/ (a + r)^{2}, η_{b} = 4br/ (b + r)^{2}, and K and E are complete elliptic integrals of the first and of the second order, respectively. The gravitational acceleration can be written as (19)where we used the relations (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 FokkerPlanck 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^{9}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^{2} ⟩/σ^{2}. 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 selfconsistent 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 selfconsistency 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 selfconsistent data for the FokkerPlanck equation is, in general, not an easy task. In this work we found several solutions for systems of flat rings that satisfy the FokkerPlanck equation and the Poisson equation simultaneously. The collision term of the FokkerPlanck 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 FokkerPlanck 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 twodimensional models and threedimensional structures. Our results provide the basis needed for further evolution studies of the FokkerPlanck 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
 Ashurov, A. E. 2004, ApJ, 127, 2154 [NASA ADS] (In the text)
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton, New Jersey: Princeton University Press) (In the text)
 Catabriga, L., Coutinho, A. L. G. A., & Franca, L. P. 2004, Int. J. Numer. Meth. Eng., 60, 1513 [CrossRef] (In the text)
 Cohn, H. 1979, ApJ, 234, 1036 [NASA ADS] [CrossRef] (In the text)
 Fiestas, J., & Spurzem, R. 2010, MNRAS, 405, 194 [NASA ADS] (In the text)
 Joshi, K. J., Nave, C. P., & Rasio, F. A. 2001, ApJ, 550, 691 [NASA ADS] [CrossRef] (In the text)
 Kalnajs, A. J. 1972, ApJ, 175, 63 [NASA ADS] [CrossRef] (In the text)
 Lass, H., & Blitzer, L. 1983, Celest. Mech. Dyn. Astron., 30, 225 [NASA ADS] [CrossRef] (In the text)
 Mestel, L. 1963, MNRAS, 126, 553 [NASA ADS] [CrossRef] (In the text)
 Nishida, M. T., Watanabe, Y., Fujiwara, T., & Kato, S. 1984, PASJ, 36, 27 [NASA ADS] (In the text)
 Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. 1957, Phys. Rev., 107, 1 [NASA ADS] [CrossRef] [MathSciNet] (In the text)
 Takahashi, K. 1995, PASJ, 47, 561 (In the text)
 Takahashi, K., Lee, H. M., & Inagaki, S. 1997, MNRAS, 292, 331 [NASA ADS] [CrossRef] (In the text)
 Theuns, T. 1996, MNRAS, 279, 827 [NASA ADS] [CrossRef] (In the text)
 Ujevic, M., & Letelier, P. S. 2005, A&A, 442, 785 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Ujevic, M., & Letelier, P. S. 2006, J. Comput. Phys., 215, 485 [NASA ADS] [CrossRef] (In the text)
 Ujevic, M., Letelier, P. S., & Vogt, D. 2011, Int. J. Mod. Phys. D, 20, 2291 [NASA ADS] [CrossRef] (In the text)
 Vogt, D., & Letelier, P. S. 2009, MNRAS, 396, 1487 [NASA ADS] [CrossRef] (In the text)
 Watanabe, Y., Inagaki, S., Nishida, M. T., Tanaka, Y. D., & Kato, S. 1981, PASJ, 33, 541 [NASA ADS] (In the text)
 Yepes, G., & DomínguezTenreiro, R. 1992, ApJ, 387, 27 [NASA ADS] [CrossRef] (In the text)
Appendix A: Explicit determination of the collision term
The collision term for the twodimensional 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 (A.1)where m and m_{a} are the mass of the test particle and the field particle, respectively, , 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. (A.2)where n is the number density of particles, and σ the onedimensional velocity dispersion of the field particles. For this particular case, it is possible to find explicit expressions for the diffusion coefficients (A.3)where , ρ = nm_{a}, ϵ_{i} = v_{i}/σ. We define the new dimensionless diffusion coefficients (A.4)where (A.5)Replacing the expressions (A.3) in the collision term (5), we find(A.6)Evaluating the derivatives we can write the previous equation as (A.7)where (A.8)Finally, the derivatives of the dimensionless diffusion coefficients are determined from the following relationships: (A.9)
where (A.10)In the case where ϵ → 0, we replace in Eq. (A.9) the following limits (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)} = 10^{9} 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 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 ⟨ v^{2} ⟩ /σ^{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. 

Open with DEXTER  
In the text 
Fig. 2 Results for the double ring (1, 2). We consider the parameters a = 4 kpc, σ = 24.80 km s^{1}, , κ = 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 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 ⟨ v^{2} ⟩/σ^{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. 

Open with DEXTER  
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}, , κ = 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 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 ⟨ v^{2} ⟩/σ^{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 v_{r} = 0, respectively. 

Open with DEXTER  
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} = 10^{9}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 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 ⟨ v^{2} ⟩ /σ^{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 v_{r} = 0, respectively. 

Open with DEXTER  
In the text 