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

Abstract

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 𝑁-body method and artificial boundary conditions for the spectral methods. With the help of fast Fourier transform, the computational complexity is only 𝑂(𝑀²(log₂𝑀)), where 𝑀 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.

1. Introduction

Let us consider the Poisson equation with singular disk surface density [1–3] ΔΦ(𝐱)=4𝜋𝐺𝜌(𝐱)=4𝜋𝐺𝜎𝑥,𝑦𝛿(𝑧),inℝ³, (1) where𝐱=(𝑥,𝑦,𝑧),𝐺,𝛿, and 𝜌(𝐱)=𝜎𝑥,𝑦𝛿(𝑧), (2) are position, gravitational constant, Dirac symbol, and density, respectively. Hereafter, we set𝐺=1without loss of generality.

One possible method of solving (1) is to discretize the partial differential equation (1) using the finite difference method [4–6]. This discretization is −−Φ_{𝑖+1,𝑗,𝑘}+2Φ_{𝑖,𝑗,𝑘}−Φ_{𝑖−1,𝑗,𝑘}(Δ𝑥)²,−−Φ_{𝑖,𝑗+1,𝑘}+2Φ_{𝑖,𝑗,𝑘}−Φ_{𝑖,𝑗−1,𝑘}Δ𝑦²,−−Φ_{𝑖,𝑗,𝑘+1}+2Φ_{𝑖,𝑗,𝑘}−Φ_{𝑖,𝑗,𝑘−1}(Δ𝑧)²=𝑓_{𝑖,𝑗,𝑘}, (3) whereΦ_{𝑖,𝑗,𝑘}=Φ(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘),𝑓_𝑖,𝑗,0=𝜎_𝑖,𝑗, and𝑓_{𝑖,𝑗,𝑘}=0for𝑘≠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𝑂(𝑀³)unknowns, where𝑀is the number of zones in one dimension. The linear complexity of this approach is at least𝑂(𝑀³), for example, [7]. It seems unfeasible that a thin gaseous disk problem can directly solve the partial differential equation.

The other approach for solving (1) uses the integral form. The solution to (1) can be represented in terms of the fundamental solution,(1/4𝜋)𝒦(𝐱), where 𝒦(𝐱)=1𝑥²+𝑦²+𝑧², (4) as Φ𝑥,𝑦,𝑧=−𝒦𝑥−𝑥,𝑦−𝑦,−𝑧𝜎𝑥,𝑦𝑑𝑥𝑑𝑦. (5) The difficulties encountered in the numerical approach for solving (1) or (5) are related to the extent of the domainℝ³. The spectral method uses trigonometric bases functions and the artificial periodic boundary conditions [8]. However, the potential is not periodic in reality. A direct calculation using the𝑁-body method [1] has the total amount of complexity 𝑂(𝑀⁴)based on the number𝑂(𝑀²)of mesh zones. For the𝑁-body method based on a uniform grid and using the FFT technique, the numerical complexity can be reduced from𝑂(𝑀⁴)to𝑂(𝑀²(log₂(𝑀))²). In other words, this produces a fast algorithm of linear complexity. Unfortunately, this method has a drawback, which is to introduce the softening parameter to overcome the singular integration _𝐷𝒦𝑥,𝑦,0𝑑𝑥𝑑𝑦, (6) where𝐷contains the origin.

The proposed method does not require the artificial periodic boundary condition and softening parameter and further improves the accuracy of the𝑁-body method. This study uses the uniform and logarithmic grid generations in Cartesian and polar coordinates, respectively. The other treatment is linear approximation for the surface density on the cell. This ensures that the proposed method is of second order in Cartesian coordinates. This study uses two examples, a𝐷₂disk [9] and its variation𝐷_2,2, to show numerically that the method achieves second-order accuracy. A comparison study also demonstrates that the proposed method is suitable for potentials on a thin disk.

Before the end of the section, another approach for an improvement of𝑁-body method is the fast multipole methods (FMM) [1, 10] which decomposes the region into local interactions and far-field interactions for every position. The FMM is efficient for fully three-dimensional calculations under the assumption that the far field is smooth. The FFM also relies on the ability to manipulate local multipole expansions for every box in the tree hierarchy [11–13]. The proposed presentation based on the uniform grids is without the assumption of smoothness in the far-field region and does not require extra efforts on the tree hierarchy.

This paper is organized as follows. Sections 2.1 and 2.2 present the numerical methods for Cartesian and spherical coordinates, respectively. Section 4.1 presents the order of accuracy of the proposed methods as verified by a family of finite disks𝐷₂[9]. Section 4.2 presents a comparison with the𝑁-body calculation. Finally, we conclude this study in Section 5.

2. Methods

2.1. Cartesian Coordinates

Assume that the numerical computation domain𝐷=[−𝐿,𝐿]×[−𝐿,𝐿]for some number𝐿>0, which contains the support of the surface density, and further discretize the region uniformly as follows. Given a positive integer𝑀, defineΔ𝑥=2𝐿/𝑀,𝑥_𝑖+1/2=−𝐿+𝑖Δ𝑥,Δ𝑦=Δ𝑥, and𝑦_𝑗+1/2=−𝐿+𝑗Δ𝑦, where𝑖,𝑗=0,…,𝑀. Next, define the center of the cell𝐷_𝑖,𝑗=[𝑥_𝑖−1/2,𝑥_𝑖+1/2]×[𝑦_𝑗−1/2,𝑦_𝑗+1/2]as𝑥_𝑖=(𝑥_𝑖−1/2+𝑥_𝑖+1/2)/2and𝑦_𝑗=(𝑦_𝑗−1/2+𝑦_𝑗+1/2)/2, where𝑖,𝑗=1,…,𝑀. Hence, the domain𝐷is discretized into the𝑀²cells.

We recall the𝑁-body method with fast Fourier transform (FFT). Represent the potential functionΦof (1) as Φ𝒦𝑥,𝑦,𝑧=−𝑥−𝑥,𝑦−𝑦,𝑧−𝑧𝜌𝑥,𝑦,𝑧𝑑𝑥𝑑𝑦𝑑𝑧, (7) Where  𝒦(𝑥,𝑦,𝑧)=1/𝑥²+𝑦²+𝑧².According to (2), the potential is Φ𝑥,𝑦,𝑧=𝒦𝑥−𝑥,𝑦−𝑦,−𝑧𝜎𝑥,𝑦𝑑𝑥𝑑𝑦. (8) Calculate (8) using a numerical approach for the plane𝑧=0assuming that the support of the surface density is compact. Based on a restriction to the calculation of potentials on the𝑧=0plane, derive the formulae for all𝑧. Equation (8) is equal to Φ=𝑥,𝑦,𝑧_{𝑀𝑖′𝑀=1𝑗′=1𝐷𝑖′′,𝑗}𝒦𝑥−𝑥,𝑦−𝑦,−𝑧𝜎𝑥,𝑦𝑑𝑥𝑑𝑦. (9) The potential at(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)for the𝑁-body method is approximated by Φ𝑥_𝑖,𝑦_𝑗,𝑧_𝑘≈^𝑀_𝑖′^𝑀_=1𝑗′=1𝒦𝑥_𝑖′−𝑥_𝑖,𝑦_𝑗′−𝑦_𝑗,−𝑧_𝑘𝜎𝑥_𝑖′,𝑦_𝑗′Δ𝑥Δ𝑦. (10) That is, _{𝐷𝑖′′,𝑗}𝒦𝑥−𝑥,𝑦−𝑦,−𝑧𝜎𝑥,𝑦≈𝒦𝑥_𝑖′−𝑥_𝑖,𝑦_𝑗′−𝑦_𝑗,−𝑧_𝑘𝜎𝑥_𝑖′,𝑦_𝑗′Δ𝑥Δ𝑦. (11) The approximated integral is of the first order and the fundamental function  𝒦  is singular whenever(𝑖^′,𝑗^′,0)=(𝑖,𝑗,𝑘). The𝑁-body method applies a softening parameter technique for𝑧=0with𝜖>0, 𝒦≈𝒦_𝜖=1𝑥²+𝑦²+𝜖² (12) to avoid the singularity.

To improve the calculation of (11), approximate the surface density𝜎on𝐷_𝑖,𝑗in (8) linearly by 𝜎𝑥,𝑦≈𝜎_𝑖,𝑗+𝛿^𝑥_𝑖,𝑗𝑥−𝑥_𝑖+𝛿^𝑦_𝑖,𝑗𝑦−𝑦_𝑗, (13) where𝜎_𝑖,𝑗=𝜎(𝑥_𝑖,𝑦_𝑗)and𝛿^𝑥_𝑖,𝑗=(𝜕/𝜕𝑥)𝜎(𝑥_𝑖,𝑦_𝑗)and𝛿^𝑦_𝑖,𝑗=(𝜕/𝜕𝑦)𝜎(𝑥_𝑖,𝑦_𝑗)are constant in the cell𝐷_𝑖,𝑗. The error of the discretization is𝑂((𝑥−𝑥_𝑖)²+(𝑦−𝑦_𝑗)²). If a higher order accuracy is required, additional terms in the Taylor expansion in (13) can be considered.

Let 𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=_{𝐷𝑖′′,𝑗}1𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘𝑑𝑥𝑑𝒦𝑦,_{𝑥𝑖−𝑖′,𝑗−𝑗′,𝑘}=_{𝐷𝑖′′,𝑗}𝑥−𝑥_𝑖′𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘𝑑𝑥𝑑𝒦𝑦,_{𝑦𝑖−𝑖′,𝑗−𝑗′,𝑘}=_{𝐷𝑖′′,𝑗}𝑦−𝑦_𝑗′𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘𝑑𝑥𝑑𝑦. (14) By substituting (13) into (9), the potential can be approximated by Φ_{𝑖,𝑗,𝑘}≈^𝑀_𝑖′^𝑀_{=1𝑗′=1𝐷𝑖′′,𝑗}𝒦𝑥−𝑥_𝑖,𝑦−𝑦_𝑗,−𝑧_𝑘×𝜎_{𝑖′,𝑗′}+𝛿^𝑥_{𝑖′,𝑗′}𝑥−𝑥_𝑖′+𝛿^𝑦_{𝑖′,𝑗′}𝑦−𝑦_𝑗′𝑑𝑥𝑑𝑦≡Φ⁰_{𝑖,𝑗,𝑘}+Φ^𝑥_{𝑖,𝑗,𝑘}+Φ^𝑦_{𝑖,𝑗,𝑘}, (15) where Φ⁰_{𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=1𝑗′=1𝜎_{𝑖′,𝑗′𝐷𝑖′′,𝑗}1𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘𝑑𝑥𝑑𝑦=^𝑀_𝑖′^𝑀_=1𝑗′=1𝜎_{𝑖′,𝑗′}𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘},Φ^𝑥_{𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=1𝑗′=1𝛿^𝑥_{𝑖′,𝑗′𝐷𝑖′′,𝑗}𝑥−𝑥_𝑖′𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘𝑑𝑥𝑑𝑦=^𝑀_𝑖′^𝑀_=1𝑗′=1𝛿^𝑥_{𝑖′,𝑗′}𝒦^𝑥_{𝑖−𝑖′,𝑗−𝑗′,𝑘},Φ^𝑦_{𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=1𝑗′=1𝛿^𝑦_{𝑖′,𝑗′𝐷𝑖′′,𝑗}𝑦−𝑦_𝑗′𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘𝑑𝑥𝑑𝑦=^𝑀_𝑖′^𝑀_=1𝑗′=1𝛿^𝑦_{𝑖′,𝑗′}𝒦^𝑦_{𝑖−𝑖′,𝑗−𝑗′,𝑘}. (16)

Fortunately, the integrals of (14) can be obtained analytically with the help of the following simple formulae:1𝑥²+𝑦²+𝑧²𝑑𝑥𝑑𝑦=𝑦ln𝑥+𝑥²+𝑦²+𝑧²+𝑥ln𝑦+𝑥²+𝑦²+𝑧²−𝑧arctan𝑥𝑦𝑧𝑥²+𝑦²+𝑧²𝑥+𝐶,𝑥²+𝑦²+𝑧²=1𝑑𝑥𝑑𝑦2𝑦𝑥²+𝑦²+𝑧²+12𝑥²+𝑧²ln𝑦+𝑥²+𝑦²+𝑧²+𝐶. (17) The value𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}is then equal to 𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=𝑦−𝑦_𝑗ln𝑥−𝑥_𝑖+𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘+𝑥−𝑥_𝑖ln𝑦−𝑦_𝑗+𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘−𝑧_𝑘arctan𝑥−𝑥_𝑖𝑦−𝑦_𝑗𝑧_𝑘𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘^{𝑥𝑖′+1/2}_{𝑥𝑖′−1/2}^{𝑦𝑗′+1/2}_{𝑦𝑗′−1/2}, (18) where the notation𝑔(𝑥)|^𝑏_𝑎=𝑔(𝑏)−𝑔(𝑎). The calculations of𝒦^𝑥_{𝑖−𝑖′,𝑗−𝑗′,𝑘}and𝒦^𝑦_{𝑖−𝑖′,𝑗−𝑗′,𝑘}are split into two parts by the identity𝑥−𝑥_𝑖′=(𝑥−𝑥_𝑖)+(𝑥_𝑖−𝑥_𝑖′) and𝑦−𝑦_𝑗′=(𝑦−𝑦_𝑗)+(𝑦_𝑗−𝑦_𝑗′), respectively. This leads to 𝒦^𝑥_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=𝑥_𝑖−𝑥_𝑖′𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}+12𝑦−𝑦_𝑗𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘+12𝑥−𝑥_𝑖²+𝑧²_𝑘×ln𝑦−𝑦_𝑗+𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘^{𝑥𝑖′+1/2}_{𝑥𝑖′−1/2}^{𝑦𝑗′+1/2}_{𝑦𝑗′−1/2},𝒦^𝑦_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=𝑦_𝑗−𝑦_𝑗′𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}+12𝑥−𝑥_𝑖𝑦−𝑦_𝑗²+𝑥−𝑥_𝑖²+𝑧²_𝑘+12𝑦−𝑦_𝑗²+𝑧²_𝑘×ln𝑥−𝑥_𝑖+𝑥−𝑥_𝑖²+𝑦−𝑦_𝑗²+𝑧²_𝑘^{𝑥𝑖′+1/2}_{𝑥𝑖′−1/2}^{𝑦𝑗′+1/2}_{𝑦𝑗′−1/2}. (19)

The forms ofΦ⁰_{𝑖,𝑗,𝑘},Φ^𝑥_{𝑖,𝑗,𝑘}, andΦ^𝑦_{𝑖,𝑗,𝑘}in (16) are a type of convolution and the computational complexity can be reduced to𝑂(𝑀log₂𝑀)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𝑂(𝑀²(log₂𝑀)²)and𝑂(𝑀³(log₂𝑀)²)for the plane𝑧=0and the whole 3D space, respectively.

2.2. 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  ℛ=[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𝑀, defineΔ𝜃=2𝜋/𝑀,0<𝛽₀<1,𝛽=𝛽₀(1−Δ𝜃),𝑟_𝑖+1/2=𝛽^𝑀−𝑖𝐿,𝜃_𝑗+1/2=𝑗Δ𝜃,𝑖,𝑗=0,…,𝑀,𝑟_𝑖=(1/2)(𝑟_𝑖−1/2+𝑟_𝑖+1/2), and𝜃_𝑗=(1/2)(𝜃_𝑗−1/2+𝜃_𝑗+1/2)where𝑖,𝑗=1,…,𝑀. The effect of𝛽₀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  ℛ  is discretized into the𝑀²cells,ℛ_𝑖𝑗=[𝑟_𝑖−1/2,𝑟_𝑖+1/2]×[𝜃_𝑗−1/2,𝜃_𝑗+1/2]for𝑖,𝑗=1,…,𝑀. For𝑗=1,…,𝑀, the cellsℛ_1,𝑗do not cover the origin and extra cellsℛ_𝑗=[0,𝑟_1/2]×[𝜃_𝑗−1/2,𝜃_𝑗+1/2]should be included. To simplify notation, denoteℛ_0,𝑗=ℛ_𝑗,𝑗=1,…,𝑀.

The potential functionΦof (1) under the assumption𝐺=1in spherical coordinates can be expressed as Φ𝒦𝑟,𝜃,𝜑=−𝑟,𝑟,𝜃−𝜃,𝜑,𝜑𝜌𝑟,𝜃,𝜑𝑟²sin𝜑𝑑𝑟𝑑𝜃𝑑𝜑, (20) where 𝒦𝑟,𝑟,𝜃−𝜃,𝜑,𝜑=1×𝑟²−2𝑟𝑟cos𝜃−𝜃sin𝜑sin𝜑+𝑟²−2𝑟𝑟cos𝜑cos𝜑^−1/2. (21)

The surface density𝜎onℛ_𝑖,𝑗in (20) can be linearly approximated by 𝜎𝑟,𝜃≈𝜎_𝑖,𝑗+𝛿^𝑟_𝑖,𝑗𝑟−𝑟_𝑖+𝛿^𝜃_𝑖,𝑗𝜃−𝜃_𝑗, (22) where𝜎_𝑖,𝑗=𝜎(𝑟_𝑖,𝜃_𝑗)and𝛿^𝑟_𝑖,𝑗=(𝜕/𝜕𝑟)𝜎(𝑟_𝑖,𝜃_𝑗)and𝛿^𝜃_𝑖,𝑗=(𝜕/𝜕𝜃)𝜎(𝑟_𝑖,𝜃_𝑗)are constants in the cellℛ_𝑖,𝑗. Equation (22) is the Taylor expansion in two dimensions and it follows that the error of the discretization is𝑂((𝑟−𝑟_𝑖)²+(𝜃−𝜃_𝑗)²).

In this section, let 𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=_{ℛ𝑖′′,𝑗}𝑟𝑟²+𝑟²_𝑖−2𝑟𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘𝑑𝑟𝑑𝜃,𝑟_𝑖𝒦^𝑟_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=_{𝑅𝑖′′,𝑗}𝑟𝑟−𝑟_𝑖′𝑟²+𝑟²_𝑖−2𝑟𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘𝑑𝑟𝑑𝜃,𝒦^𝜃_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=_{𝑅𝑖′′,𝑗}𝑟𝜃−𝜃_𝑗′𝑟²+𝑟²_𝑖−2𝑟𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘𝑑𝑟𝑑𝜃. (23)(24)(25) Similarly for𝒦⁰_{𝑖,𝑗−𝑗}^′_,𝑘,𝒦^𝑟_{𝑖,𝑗−𝑗}^′_,𝑘, and𝒦^𝜃_{𝑖,𝑗−𝑗}^′_,𝑘on the central region ℛ_0,𝑗. The term𝑟_𝑖in (24) is for the formulation of a convolution type. According to (20) and (22), Φ_{𝑖,𝑗,𝑘}≈^𝑀_𝑖′^𝑀_{=0𝑗′=1𝑅𝑖′′,𝑗}𝒦𝑟,𝑟_𝑖,𝜃−𝜃_𝑗,𝜑_𝑘,𝜋2×𝜎_{𝑖′,𝑗′}+𝛿^𝑟_{𝑖′,𝑗′}𝑟−𝑟_𝑖′+𝛿^𝜃_{𝑖′,𝑗′}𝜃−𝜃_𝑗′𝑟𝑑𝑟𝑑𝜃∶=Φ^0,𝑝_{𝑖,𝑗,𝑘}+Φ^𝑟_{𝑖,𝑗,𝑘}+Φ^𝜃_{𝑖,𝑗,𝑘}, (26) where Φ^0,𝑝_{𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=0𝑗′=1𝜎_{𝑖′,𝑗′}×_{𝑅𝑖′′,𝑗}𝑟𝑟²+𝑟²_𝑖−2𝑟𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘𝑑𝑟𝑑Φ𝜃,_{𝑟𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=0𝑗′=1𝛿^𝑟_{𝑖′,𝑗′}×_{𝑅𝑖′′,𝑗}𝑟𝑟−𝑟_𝑖′𝑟²+𝑟²_𝑖−2𝑟𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘𝑑𝑟𝑑Φ𝜃,_{𝜃𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=0𝑗′=1𝛿^𝜃_{𝑖′,𝑗′}×_{𝑅𝑖′′,𝑗}𝑟𝜃−𝜃_𝑗′𝑟²+𝑟²_𝑖−2𝑟𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘𝑑𝑟𝑑𝜃. (27) Equation (27) can then be rewritten as Φ^0,𝑝_{𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=1𝑗′=1𝜎_{𝑖′,𝑗′}𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}+^𝑀_𝑗′=1𝜎_0,𝑗′𝒦⁰_{𝑖,𝑗−𝑗}^′_,𝑘,Φ^𝑟_{𝑖,𝑗,𝑘}=𝑟_𝑖^𝑀_𝑖′^𝑀_=1𝑗′=1𝛿^𝑟_{𝑖′,𝑗′}𝒦^𝑟_{𝑖−𝑖,𝑗−𝑗′,𝑘}+𝑟_𝑖^𝑀_𝑗′=1𝛿^𝑟_0,𝑗′𝒦^𝑟_{𝑖,𝑗−𝑗}^′_,𝑘,Φ^𝜃_{𝑖,𝑗,𝑘}=^𝑀_𝑖′^𝑀_=1𝑗′=1𝛿^𝜃_{𝑖′,𝑗′}𝒦^𝜃_{𝑖−𝑖′,𝑗−𝑗′,𝑘}+^𝑀_𝑗′=1𝛿^𝜃_0,𝑗′𝒦^𝜃_{𝑖,𝑗−𝑗}^′_,𝑘. (28)

Define𝐹(𝑟,̃𝜃,𝜑)=1+𝑟̃²−2𝑟̃cos(𝜃)sin(𝜑), where𝑟̃is the dimensionless radius. The evaluation of (23), (24), and (25) can be obtained with the help of the following simple integrals: 𝑟𝑟²+𝑟²−2𝑟(𝜃)𝜑𝑑𝑟cossin𝑟=𝑟𝐹𝑟𝑟(𝜃)𝜑(𝜃)𝜑+,𝜃,𝜑+𝑟cossinln−cossin𝑟𝑟+𝐹𝑟𝑟,𝜃,𝜑+𝐶∶=𝑟𝐻₁𝑟𝑟,𝜃,𝜑+𝐶,𝑟²𝑟²+𝑟²−2𝑟(𝜃)𝜑𝑑𝑟cossin𝑟=𝑟²12𝑟𝑟+32cos(𝜃)sin𝜑𝐹𝑟𝑟,𝜃,𝜑+𝑟²32cos²(𝜃)sin²𝜑−12×ln−cos(𝜃)sin𝜑+𝑟𝑟+𝐹𝑟𝑟,𝜃,𝜑+𝐶∶=𝑟²𝐻₂𝑟𝑟,𝜃,𝜑+𝐶. (29)

The definitions of𝑟_𝑖′+1/2and𝑟_𝑖lead to 𝑟_𝑖′+1/2𝑟_𝑖=2𝛽^{𝑖−𝑖′},𝑟1+𝛽𝑟^𝑖′_𝑖=𝛽^{𝑖−𝑖′}. (30) Calculate the value of the integral 𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}=^{𝜃𝑗′+1/2}_{𝜃𝑗′−1/2}^{𝑟𝑖′+1/2}_{𝑟𝑖′−1/2}𝑟𝑟²+𝑟²_𝑖−2𝑟𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘𝑑𝑟𝑑𝜃=^{𝜃𝑗′+1/2}_{𝜃𝑗′−1/2}𝑟_𝑖𝐹𝑟𝑟_𝑖,𝜃−𝜃_𝑗,𝜑_𝑘+𝑟_𝑖cos𝜃−𝜃_𝑗sin𝜑_𝑘×ln−cos(𝜃)sin𝜑_𝑘+𝑟𝑟_𝑖+𝐹𝑟𝑟_𝑖,𝜃−𝜃_𝑗,𝜑_𝑘^{𝑟𝑖′+1/2}_{𝑟𝑖′−1/2}𝑑𝜃. (31) 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𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}as follows: 𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘}≈12𝑟_𝑖𝐻₁𝑟_𝑖′+1/2𝑟_𝑖,𝜃_𝑗′+1/2−𝜃_𝑗,𝜑_𝑘−𝑟_𝑖𝐻₁𝑟_{𝑖′−1/2}𝑟_𝑖,𝜃_{𝑗′−1/2}−𝜃_𝑗,𝜑_𝑘+𝑟_𝑖𝐻₁𝑟_𝑖′+1/2𝑟_𝑖,𝜃_{𝑗′−1/2}−𝜃_𝑗,𝜑_𝑘−𝑟_𝑖𝐻₁𝑟_{𝑖′−1/2}𝑟_𝑖,𝜃_{𝑗′−1/2}−𝜃_𝑗,𝜑_𝑘×𝜃_𝑗′+1/2−𝜃_{𝑗′−1/2}∶=𝑟_𝑖𝐻₁𝑟𝑟_𝑖,𝜃−𝜃_𝑗,𝜑_𝑘^{𝑟𝑖′+1/2}_{𝑟𝑖′−1/2}^{𝜃𝑗′+1/2}_{𝜃𝑗′−1/2}, (32) where the notation𝑓(⋅)]^𝑏_𝑎=(1/2)(𝑓(𝑏)+𝑓(𝑎))(𝑏−𝑎). Similarly, 𝒦^𝑟_{𝑖−𝑖′,𝑗−𝑗′,𝑘}≈𝑟²_𝑖𝐻₂𝑟𝑟_𝑖,𝜃−𝜃_𝑗,𝜑_𝑘^{𝑟𝑖′+1/2}_{𝑟𝑖′−1/2}^{𝜃𝑗′+1/2}_{𝜃𝑗′−1/2}−𝑟_𝑖′𝑟_𝑖𝒦⁰_{𝑖−𝑖′,𝑗−𝑗′,𝑘},𝒦^𝜃_{𝑖−𝑖′,𝑗−𝑗′,𝑘}≈𝜃−𝜃_𝑗𝑟_𝑖𝐻₁𝑟𝑟_𝑖,𝜃−𝜃_𝑗,𝜑_𝑘^{𝑟𝑖′+1/2}_{𝑟𝑖′−1/2}^{𝜃𝑗′+1/2}_{𝜃𝑗′−1/2}. (33)

3. Convergence Results

Let us denote|𝐴|as the area of the region𝐴inℝ²and discuss the convergence of the proposed method. The integrals of  𝒦(𝑥,𝑦,0)over a region  𝑅={(𝑥,𝑦)∶𝑥²+𝑦²≤𝑎}and𝐷={(𝑥,𝑦)∶−𝑎≤𝑥≤𝑎,−𝑎≤𝑦≤𝑎}, for some constant𝑎, are _𝑅𝒦1𝑥,𝑦,0𝑑𝑥𝑑𝑦=21𝑎=2|𝑅|𝜋=𝐶|𝑅|,_𝐷𝒦𝑥,𝑦,0𝑑𝑥𝑑𝑦=𝑦ln𝑥+𝑥²+𝑦²+𝑥ln𝑦+𝑥²+𝑦^2𝑎_−𝑎^𝑎_−𝑎=4𝑎ln3+2√2=𝐶|𝐷|, (34)(35) where𝐶is a constant. Furthermore, if𝐷=𝑅=ℝ², then this integral is∞.

Lemma 1. If two compact disk densities 𝜎₁ and 𝜎₂ satisfy |𝜎₁(𝑥,𝑦)−𝜎₂(𝑥,𝑦)|<ℎ then their corresponding potentials Φ₁ and Φ₂, respectively, satisfyΦ₁𝑥,𝑦−Φ₂𝑥,𝑦≤𝐶ℎ, (36) where𝐶is a positive constant and depends on the supports of𝜎₁and𝜎₂.

Proof. Let the supports of𝜎₁and𝜎₂be contained in a finite region𝐷. Following (5) and (35), it yields that Φ₁𝑥,𝑦−Φ₂=𝑥,𝑦_𝐷𝒦𝑥−𝑥,𝑦−𝑦,0𝜎₁𝑥,𝑦−𝜎₂𝑥,𝑦𝑑𝑥𝑑𝑦≤_𝐷𝒦𝑥−𝑥,𝑦−𝑦,0𝜎₁𝑥,𝑦−𝜎₂𝑥,𝑦𝑑𝑥𝑑𝑦≤𝐶ℎ. (37) 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𝜎_𝑛(𝑥,𝑦), 𝜎_𝑛𝑥,𝑦−𝜎𝑥,𝑦≤𝐶(Δ𝑥)^𝑠. (38)

The potential calculated numerically is denoted byΦ_𝑛and the analytic solution isΦ. We estimate Φ_𝑛≤𝑥,𝑦−Φ𝑥,𝑦_{𝑀𝑀𝑖=1𝑗=1𝐷𝑖,𝑗}𝒦𝑥−𝑥,×𝜎𝑦−𝑦,0₁𝑥,𝑦−𝜎₂𝑥,𝑦𝑑𝑥𝑑𝑦≤𝐶(Δ𝑥)^2𝑀𝑀_{𝑖=1𝑗=1𝐷𝑖,𝑗}𝒦𝑥−𝑥,𝑦−𝑦,0𝑑𝑥𝑑𝑦≤𝐶(Δ𝑥)^𝑠_𝐷𝒦𝑥−𝑥,𝑦−𝑦,0𝑑𝑥𝑑𝑦≤𝐶(Δ𝑥)^𝑠. (39)

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 𝑂(ℎ^𝑠).

4. Numerical Results

This study verifies that the proposed method achieved is of second-order accuracy through the following two examples: a𝐷₂disk [9] and a nonaxisymmetric disk𝐷_2,2consisting of two superposed𝐷₂disks. The𝐷₂disk has the surface density Σ_𝐷2(𝑅;𝛼)=𝜎₀𝑅1−𝛼^23/2,for𝑅<𝛼,0,for𝑅>𝛼, (40) where𝑅=𝑥²+𝑦²and𝛼is a given constant. The corresponding potential is Φ_𝐷2(𝑅,𝑧;𝛼)=−𝜋𝜎₀𝐺64𝛼³68𝛼⁴−8𝛼²𝑅²+16𝛼²𝑧²+3𝑅⁴−24𝑅²𝑧²+8𝑧⁴sin⁻¹𝑓₁−𝑓₂+√2𝑅2𝛼18𝛼²−9𝑅²+26𝑧²𝑓₁𝑓₂−𝑓₃−√2|𝑧|58𝛼²−55𝑅²+50𝑧²𝑓₁𝑓₂+𝑓₃, (41) where 𝑓₁=𝑧²+(𝑅+𝛼)²,𝑓₂=𝑧²+(𝑅−𝛼)²,𝑓₃=𝛼²−𝑅²−𝑧². (42)

4.1. Order of Accuracy

This study investigates the order of accuracy of the proposed method in norms𝐿¹,𝐿², and𝐿^∞to demonstrate the convergence in total variation, energy, and pointwise senses, respectively. These norms for a function𝑓on a domainΩare defined as 𝑓_𝑝=_Ω𝑓(𝐱)^𝑝𝑑𝐱^1/𝑝,for𝑝=1,2,𝑓_∞=ess_Ωsup𝑓(𝐱). (43) Without loss of generality for studying the order of accuracy, set the computational domainΩ=[−1,1]×[−1,1], and 𝜎₀=𝐺=1. Set𝛼=0.25for a𝐷₂disk simulation and𝛼=0.4for a𝐷_2,2disk simulation, whereΣ_𝐷2,2=Σ_𝐷2(𝑅₁;𝛼)+Σ_𝐷2(𝑅₂;𝛼)with𝑅₁=(𝑥−0.2)²+𝑦²and𝑅₂=(𝑥+0.2)²+𝑦². Figure 1 presents the profiles of the surface density, the potential on𝑧=0plane and the residual between analytic and numerical solutions for𝑁=512. In 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 bylog₂(𝐸^𝑝(2^𝑘−1)/𝐸^𝑝(2^𝑘))for𝑘=6to9and similarly for𝑂^𝑝. These errors show that the proposed method achieves nearly second-order accuracy. More precisely, the order of convergence is approximately1.9or2.0.

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 𝐷₂ (up table) and 𝐷_2,2 (down table) disks. This table shows that the accuracy is nearly second order.

(a)

(b)

(c)

(d)

(e)

(f)

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.

Continue to use the𝐷₂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.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.

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 𝐷₂ disk. This table shows that the accuracy is only first order.

(a)

(b)

(c)

Figure 2 

The numerical solutions of the𝐷2disk for𝑁=512to 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.

4.2. A Comparison Study

Start with Φ𝒦𝑥,𝑦,0=−𝐺𝑥−𝑥,𝑦−𝑦,0𝜎𝑥,𝑦𝑑𝑥𝑑𝑦=−𝐺^𝑁𝑁_{𝑖=1𝑗=1𝐷𝑖,𝑗}𝒦𝑥−𝑥,𝑦−𝑦,0𝜎𝑥,𝑦𝑑𝑥𝑑𝑦. (44) And introduce a softening parameter𝜖to approximate _{𝐷𝑖,𝑗}𝒦𝑥−𝑥,𝑦−𝑦𝜎𝑥,𝑦𝐺≈−𝜖²+𝑥_𝑖′−𝑥_𝑖²+𝑦_𝑗′−𝑦_𝑗^{2𝐷𝑖,𝑗}𝜎𝑥,𝑦𝑑𝑥𝑑𝑦. (45) For polar coordinates [1], the value of𝒦can be approximated by 𝒦_{𝑖′−𝑖,𝑗′−𝑗}𝐺∶=−2cosh𝑢_𝑖′−𝑢_𝑖−cos𝜃_𝑗′−𝜃_𝑗, (46) where𝑢_𝑖′=ln(𝑥_𝑖′)and𝑢_𝑖=ln(𝑥_𝑖). Note that when(𝑖^′,𝑗^′)=(𝑖,𝑗),   𝒦  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𝐷₂and𝐷_2,2disks in the previous section. For Cartesian coordinates, choose the softening parameters as the mesh size𝜖=Δ𝑥. Table 3 shows the errors for the disks𝐷₂ and 𝐷_2,2 that the accuracy when using the softening parameter approach for the𝐷₂and𝐷_2,2disks is of first order. This comparison confirms that the proposed method is more accurate and has a higher order of accuracy.

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) disk. This table shows that the accuracy is only first order.

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 Φ𝜎𝑥,𝑦,0=−𝐺𝑥,𝑦(𝑥−𝑥)²+𝑦−𝑦²𝑑𝑥𝑑𝑦 (47) or Φ()𝜎𝑟,𝜃,0=−𝐺𝑟,𝜃𝑟²+𝑟²−2𝑟𝑟cos𝜃−𝜃𝑟𝑑𝑟𝑑𝜃. (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 𝑒^3𝑢/2𝜎𝑒^𝑢=1,𝜃4𝜋²_𝑚^∞_−∞𝐴_𝑚(𝛼)𝑒^{𝑖(𝑚𝜃+𝛼𝑢)}𝑑𝛼,𝑒^𝑢/2Φ𝑒^𝑢1,𝜃=−2𝜋𝐺_𝑚^∞_−∞𝐾(𝛼,𝑚)𝐴_𝑚(𝛼)exp[𝑖(𝑚𝜃+𝛼𝑢)]𝑑𝛼, (49)(50) where𝐾is real and positive and is defined as 𝐾()≡1𝛼,𝑚2Γ[(𝑚+1/2+𝑖𝛼)/2]Γ[(𝑚+1/2−𝑖𝛼)/2]Γ)]Γ)].[(𝑚+3/2+𝑖𝛼/2[(𝑚+3/2−𝑖𝛼/2 (51) 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𝐷₂disk using the polar coordinates. In this special case, compute𝑚=0and truncate 𝐴₀(𝛼)=⁰_−∞𝑒^3𝑢/2𝜎𝑒^𝑢𝑒^{−𝑖𝛼𝑢}≈𝑑𝑢_0𝑢min𝑒^3𝑢/2𝜎𝑒^𝑢𝑒^{−𝑖𝛼𝑢}𝑑𝑢, (52) where the value𝑢_minis 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.

(a)

(b)

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.

4.3. 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 Example1in 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.

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 𝐷₂ disk. This table shows that the accuracy is nearly second order.

(a)

(b)

Figure 4 

The simulation of a𝐷2disk for𝑁=512. The profiles are the potential on𝑧-axis (a) and the residuals between analytic and numerical solutions (b) in natural common logarithm.

5. 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 𝑂(𝑁²(log₂(𝑁)²)), which is the same order as the𝑁-body method. The proposed method also achieves second-order accuracy.

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The author is grateful to anonymous referees for their valuable comments. This work was supported in part by the Ministry of Science and Technology, Taiwan, under the Grant NSC 101-2115-M-030-004.