A Numerical Convolution Representation of Potential for a Disk Surface Density in 3 D

This study presents a numerical convolution representation of a potential induced from a disk of surface density, which has often been investigated in the center region of galaxies. The advantage of this representation is to release the softening length for the N-body method and artificial boundary conditions for the spectral methods. With the help of fast Fourier transform, the computational complexity is only O(M(log 2 M)), where M is the number of zones in one dimension. Numerical results show an almost second order of accuracy on a Cartesian coordinate system. A comparison study also demonstrates that this method can calculate the potential for disk surface density based on the uniform grids.

One possible method of solving (1) is to discretize the partial differential equation (1) using the finite difference method [4][5][6].This discretization is where Φ ,, = Φ(  ,   ,   ),  ,,0 =  , , and  ,, = 0 for  ̸ = 0, based on the uniform mesh grids (  ,   ,   ) with the mesh sizes Δ, Δ, and Δ in , , and  directions, respectively.In this approach, the values of the potential on imposed boundaries should be calculated and a fully 3D calculation must be performed.Any numerical solutions of this partial differential problem involve ( 3 ) unknowns, where  is the number of zones in one dimension.The linear complexity of this approach is at least ( 3 ), for example, [7].It seems unfeasible that a thin gaseous disk problem can directly solve the partial differential equation.

Methods
Based on a restriction to the calculation of potentials on the  = 0 plane, derive the formulae for all .Equation ( 8) is equal to The potential at (  ,   ,   ) for the -body method is approximated by That is, The approximated integral is of the first order and the fundamental function K is singular whenever (  ,   , 0) = (, , ).The -body method applies a softening parameter technique for  = 0 with  > 0, to avoid the singularity.

Let
By substituting ( 13) into ( 9), the potential can be approximated by where Fortunately, the integrals of ( 14) can be obtained analytically with the help of the following simple formulae: The value K 0 −  ,−  , is then equal to The forms of Φ 0 ,, , Φ  ,, , and Φ  ,, in ( 16) are a type of convolution and the computational complexity can be reduced to (log 2 ) with the help of FFT in the  and  directions.In the  direction, the computational complexity is ().This implies that the total complexity of this method is ( 2 (log 2 ) 2 ) and ( 3 (log 2 ) 2 ) for the plane  = 0 and the whole 3D space, respectively.

Spherical and Polar Coordinates.
While the method appears to be second-order-accurate for discretization on a Cartesian grid, it loses second-order accuracy using spherical polar coordinates.An algorithm for solving the Poisson equation in spherical polar coordinates has been used extensively in protostellar collapse and fragmentation calculations, where often a thin disk formed at the central core of the collapsing cloud [14].The algorithm implemented in [15] is subsequently modified, see [16], to deal with filamentary density distributions.This method in [16] employed finite differences for the radial dependence of the potential and spherical harmonics for the angular dependence, resulting in a second-order method through convergence testing.
We are going to develop our approach for spherical coordinates (, , ), which can be reduced to polar coordinates ( = /2).The singular integral disappears, but the high order of accuracy is not attained because there is no explicit form for the integral of an elliptic function.Although the proposed calculation of potential has some drawbacks, it is still better than -body calculations and the method in [17].
The support of the surface density is compact, contained in a region R = [0, ] × [0, 2] for some number  > 0. Confine the discretization on the radial region in logarithmic form and the azimuthal region uniformly as follows to achieve a fast calculation.Given a positive integer , where ,  = 1, . . ., .The effect of  0 is to refine the mesh.Here, the proposed method for polar coordinates is of first order because a singular integration occurs (see below).Furthermore, the region R is discretized into the for ,  = 1, . . ., .For  = 1, . . ., , the cells R 1, do not cover the origin and extra ] should be included.To simplify notation, denote R 0, = R ,  = 1, . . ., .
The potential function Φ of (1) under the assumption  = 1 in spherical coordinates can be expressed as where The surface density  on R , in (20) can be linearly approximated by where  , = (  ,   ) and   , = (/)(  ,   ) and   , = (/)(  ,   ) are constants in the cell R , .Equation ( 22) is the Taylor expansion in two dimensions and it follows that the error of the discretization is (( −   ) 2 + ( −   ) 2 ).
In this section, let where Equation ( 27) can then be rewritten as Define (r, , ) = √1 + r2 − 2r cos() sin(), where r is the dimensionless radius.The evaluation of ( 23), (24), and (25) can be obtained with the help of the following simple integrals: The definitions of    +1/2 and   lead to Calculate the value of the integral . ( The last integral in this equation is an elliptic integral, and this study uses a trapezoidal rule for its evaluation.This numerical approach approximates the value K 0 −  ,−  , as follows: where the notation (⋅)]   = (1/2)(() + ())( − ).Similarly,
Since || is finite and the lemma is obtained.
The order of accuracy between analytic and numerical solutions can be observed from the truncated errors of the disk surface (, ) and truncated density   (, ),       (, ) −  (, )     ≤ (Δ)  .
The potential calculated numerically is denoted by Φ  and the analytic solution is Φ.We estimate We deduce the following theorem by aforementioned argument and Lemma 1.

Theorem 2.
If the truncated errors of the Taylor series of the disk surface density (, ) are (ℎ  ), where ℎ is the mesh size in one dimension and  > 0, then the errors between the numerical and analytic potential also are (ℎ  ).

Numerical Results
This study verifies that the proposed method achieved is of second-order accuracy through the following two examples: a  2 disk [9] and a nonaxisymmetric disk  2,2 consisting of two superposed  2 disks.The  2 disk has the surface density where  = √ 2 +  2 and  is a given constant.The corresponding potential is where (42)  Σ  2 ( 1 ; ) + Σ  2 ( 2 ; ) with  1 = √ ( − 0.2) 2 +  2 and  2 = √ ( + 0.2) 2 +  2 .Figure 1 presents the profiles of the surface density, the potential on  = 0 plane and the residual between analytic and numerical solutions for  = 512.In  (c) Figure 2: The numerical solutions of the  2 disk for  = 512 to investigate the calculation of potentials in polar coordinates.From left to right, the profiles are surface density, the potential on the plane  = 0, and the difference between analytic and numerical solutions in natural common logarithm, respectively.Table 1, column   represents the error between the analytic and numerical solutions in -norm,  = 1,2, and ∞.Column   in Table 1 represents the order of accuracy as measured by log 2 (  (2 −1 )/  (2  )) for  = 6 to 9 and similarly for   .These errors show that the proposed method achieves nearly second-order accuracy.More precisely, the order of convergence is approximately 1.9 or 2.0.
Continue to use the  2 disk as an example and a unit disk (0, 1) = Ω = [0, 1] × [0, 2] as the computational domain to investigate the potential in polar coordinates.Set the value  0 = 0.99. Figure 2 presents the profiles of the surface density, potential, and the difference between analytic and numerical solutions for  = 512.According to Table 2, the order of accuracy is only approximately 1.Although the surface density at the origin is smooth, the singular elliptic integral introduces significant error there.

A Comparison
For polar coordinates [1], the value ofKcan be approximated by where    = ln(   ) and   = ln(  ).Note that when (  ,   ) = (, ), K is undefined.Previous research has presented a way to avoid the singularity problem [1].However, the proposed method avoids the singularity problem by directly evaluating the forces, thereby increasing the order of accuracy.This study compares the proposed method with the body method using the simulations of the  2 and  2,2 disks in the previous section.For Cartesian coordinates, choose the softening parameters as the mesh size  = Δ.Table 3 shows the errors for the disks  2 and  2,2 that the accuracy when using the softening parameter approach for the  2 and  2,2 disks is of first order.This comparison confirms that the proposed method is more accurate and has a higher order of accuracy.
This study calculates potentials on a disk of surface density with as few restrictions as possible.Alternatively, it is possible to solve the reduced equation given by (48) The approach in [17] transforms the polar coordinate (, ) into the coordinate (, ) by setting  =   or  = ln().The potential-density pair in terms of the reduced surface density and reduced potential is given in [17], and it is where  is real and positive and is defined as The bounded unit disk (0, 1) = [0, 1] × [0, 2] can be transformed to the unbounded domain  = (−∞, 0] × [0, 2] and apply this method to the  2 disk using the polar coordinates.In this special case, compute  = 0 and truncate Figure 3: The variation of the potential with respect to radius using the Kalnajs method (a) and the proposed method (b).The residuals appear in the small window in each panel and show that the Kalnajs method produces significant errors near the origin, which are eliminated by the proposed method.
Table 4: The errors of the potential on the -axis and order accuracy of the proposed method on Cartesian coordinates for various zones,  = 2  from  = 5 to 9 for the  2 disk.This table shows that the accuracy is nearly second order.where the value  min is to approximately −∞.The truncation process produces a hole in the unit disk and can introduce significant errors at the origin.Given a positive integer  and based on the discretization for the radial region in the previous subsection, calculate (52) and (50) using the trapzoidal rule.Figure 3 shows the variation of the potential with respect to radius.The profile on the left panel shows that the numerical and analytic solutions for the Kalnajs method agree well except for those close to the origin of  = 512.The small window embedded within the panel zooms in on the residuals between the numerical and analytic solutions to the interval [0, 0.3].The truncated portion contributes to significant errors near the origin.In contrast, the application of the proposed method to the calculation of potentials leads to the results shown in the right panel of Figure 3.Although the singular integration remains because of the unbounded domain, the proposed method is preferable for both Cartesian and polar coordinates because a hole near the origin is not introduced.

Potentials on 𝑧-Axis.
Although this study concentrates on the calculation of potentials for a disk, the proposed method can also calculate the forces on the whole space compared with [18].The simulation of Example 1 in Section 4.1 shows that the order of accuracy for the -axis is of nearly second order (see Table 4), and Figure 4 shows the profiles, the potential on -axis and the residual between analytic and numerical solutions.

Conclusion
This study presents an improved method of the -body calculation for solving Poisson equation induced by a disk of surface density.The proposed method does not require artificial boundary conditions and it also does not require a softening length parameter and the computational complexity is ( 2 (log 2 () 2 )), which is the same order as the body method.The proposed method also achieves secondorder accuracy.

Figure 1 :
Figure 1: A simulation of the  2 (up row) and  2,2 (down row) disks for  = 512.The profiles are surface density (left), the potential on the plane  = 0 (middle), and the residual between analytic and numerical solutions in natural common logarithm, respectively.

Lemma 1 .
If two compact disk densities  1 and  2 satisfy | 1 (, ) −  2 (, )| < ℎ then their corresponding potentials Φ 1 and Φ 2 , respectively, satisfy This study investigates the order of accuracy of the proposed method in norms  1 ,  2 , and  ∞ to demonstrate the convergence in total variation, energy, and pointwise senses, respectively.These norms for a function  on a domain Ω are defined as 4.1.Order of Accuracy.

Table 1 :
The errors and order accuracy of the proposed method on Cartesian coordinates for various number of zones,  = 2  from  = 5 to 9 for the  2 (up table) and  2,2 (down table) disks.This table shows that the accuracy is nearly second order.
Surface density (x, y)

Table 2 :
The errors and order accuracy of the proposed method on polar coordinates for various number of zones,  = 2  from  = 5 to 9 for the  2 disk.This table shows that the accuracy is only first order.

Table 3 :
The top (bottom) table demonstrates the errors and order accuracy of the -body simulations in Cartesian coordinates for various number of zones,  = 2  from  = 5 to 9 for the  2 ( 2,2 ) disk.This table shows that the accuracy is only first order.