Self-gravitational Force Calculation of High-Order Accuracy for Infinitesimally Thin Gaseous Disks

TL;DR

This paper introduces a fast, high-order numerical algorithm for calculating self-gravitational forces in infinitesimally thin disks, crucial for galactic and protoplanetary disk studies, using convolution and FFT techniques.

Contribution

The authors develop a linear-complexity, high-order accurate method leveraging B-splines and FFT, avoiding artificial boundary conditions for force calculations in thin disks.

Findings

The method achieves significant accuracy improvements with higher order.

Numerical results agree well with analytic solutions for generalized disks.

The approach is efficient and stable for practical astrophysical simulations.

Abstract

Self-gravitational force calculation for infinitesimally thin disks is important for studies on the evolution of galactic and protoplanetary disks. Although high-order methods have been developed for hydrodynamic and magneto-hydrodynamic equations, high-order improvement is desirable for solving self-gravitational forces for thin disks. In this work, we present a new numerical algorithm that is of linear complexity and of high-order accuracy. This approach is fast since the force calculation is associated with a convolution form, and the fast calculation can be achieved using Fast Fourier Transform. The nice properties, such as the finite supports and smoothness, of B-splines are exploited to stably interpolate a surface density and achieve a high-order accuracy in forces. Moreover, if the mass distribution of interest is exclusively confined within a calculation domain, the method does not require artificial boundary values to be specified before the force calculation. To validate the proposed algorithm, a series of numerical tests, ranging from 1st- to 3rd-order implementations, are performed and the results are compared with analytic expressions derived for 3rd- and 4th-order generalized Maclaurin disks. We conclude that the improvement on the numerical accuracy is significant with the order of the method, with only little increase of the complexity of the method.

Tables

Table 1. Table 1: Each triplet ( α , β , γ ) 𝛼 𝛽 𝛾 (\alpha,\beta,\gamma) in the table represents γ ​ F ~ ​ ( x ℓ , x r , y ℓ , y ​ r , α , β ) 𝛾 ~ 𝐹 subscript 𝑥 ℓ subscript 𝑥 𝑟 subscript 𝑦 ℓ 𝑦 𝑟 𝛼 𝛽 \gamma\tilde{F}(x_{\ell},x_{r},y_{\ell},yr,\alpha,\beta) defined in Eq. ( 28 ), and γ 𝛾 \gamma is the associated coefficient. The associated ranges of integration in x 𝑥 x and y 𝑦 y are indicated in the first row and the first column, respectively.

range $(x_{\ell }/y_{\ell },x_r/y_r)$	$(x_{i^{\prime },0},x_{i^{\prime },1})$
$(y_{j^{\prime },0},y_{j^{\prime },1})$	$(1,0,1)$

Table 2. Table 2: Each triplet ( α , β , γ ) 𝛼 𝛽 𝛾 (\alpha,\beta,\gamma) in the table represents γ ​ F ~ ​ ( x ℓ , x r , y ℓ , y ​ r , α , β ) 𝛾 ~ 𝐹 subscript 𝑥 ℓ subscript 𝑥 𝑟 subscript 𝑦 ℓ 𝑦 𝑟 𝛼 𝛽 \gamma\tilde{F}(x_{\ell},x_{r},y_{\ell},yr,\alpha,\beta) defined in Eq. ( 28 ), and γ 𝛾 \gamma is the associated coefficient. The associated ranges of integration in x 𝑥 x and y 𝑦 y are indicated in the first row and the first column, respectively.

range $(x_{\ell }/y_{\ell },x_r/y_r)$	$(x_{i^{\prime },0},x_{i^{\prime },1})$	$(x_{i^{\prime },1},x_{i^{\prime },2})$
$(y_{j^{\prime },0},y_{j^{\prime },1})$	$(2,1,1)$, $(1,1,-x_{i^{\prime },0})$,	$(2,1,-1)$, $(2,0,y_{j^{\prime },0})$
	$(2,0,-y_{j^{\prime },0})$, $(1,0,x_{i^{\prime },0}{y}_{j^{\prime },0})$	$(1,1,x_{i^{\prime },2})$, $(1,0,-x_{i^{\prime },2}{y}_{j^{\prime },0})$
$(y_{j^{\prime },1},y_{j^{\prime },2})$	$(2,1,-1)$, $(2,0,y_{j^{\prime },2})$	$(2,1,1)$, $(2,0,-y_{j^{\prime },2})$
	$(1,1,x_{i^{\prime },0})$, $(1,0,-x_{i^{\prime },0}{y}_{j^{\prime },2})$	$(1,1,-x_{i^{\prime },2})$, $(1,0,x_{i^{\prime },2}{y}_{j^{\prime },2})$

Table 3. Table 3: Each triplet ( α , β , γ ) 𝛼 𝛽 𝛾 (\alpha,\beta,\gamma) in the table represents γ ​ F ~ ​ ( x ℓ , x r , y ℓ , y ​ r , α , β ) 𝛾 ~ 𝐹 subscript 𝑥 ℓ subscript 𝑥 𝑟 subscript 𝑦 ℓ 𝑦 𝑟 𝛼 𝛽 \gamma\tilde{F}(x_{\ell},x_{r},y_{\ell},yr,\alpha,\beta) defined in Eq. ( 28 ), and γ 𝛾 \gamma is the associated coefficient. The associated ranges of integration in x 𝑥 x and y 𝑦 y are indicated in the first row and the first column, respectively. To avoid the cluttering symbols, we have suppressed the subscripts i ′ superscript 𝑖 ′ i^{\prime} and j ′ superscript 𝑗 ′ j^{\prime} in the triplets. In this table, A x = x 0 + x 1 + x 2 + x 3 subscript 𝐴 𝑥 subscript 𝑥 0 subscript 𝑥 1 subscript 𝑥 2 subscript 𝑥 3 A_{x}=x_{0}+x_{1}+x_{2}+x_{3} , B x = x 0 ​ x 2 + x 1 ​ x 3 subscript 𝐵 𝑥 subscript 𝑥 0 subscript 𝑥 2 subscript 𝑥 1 subscript 𝑥 3 B_{x}=x_{0}x_{2}+x_{1}x_{3} , A y = y 0 + y 1 + y 2 + y 3 subscript 𝐴 𝑦 subscript 𝑦 0 subscript 𝑦 1 subscript 𝑦 2 subscript 𝑦 3 A_{y}=y_{0}+y_{1}+y_{2}+y_{3} , B y = y 0 ​ y 2 + y 1 ​ y 3 subscript 𝐵 𝑦 subscript 𝑦 0 subscript 𝑦 2 subscript 𝑦 1 subscript 𝑦 3 B_{y}=y_{0}y_{2}+y_{1}y_{3} .

range $(x_{\ell }/y_{\ell },x_r/y_r)$	$(x_{i^{\prime },0},x_{i^{\prime },1})$	$(x_{i^{\prime },1},x_{i^{\prime },2})$	$(x_{i^{\prime },2},x_{i^{\prime },3})$
$(y_{j^{\prime },0},y_{j^{\prime },1})$	$(3,2,1)$, $(3,1,-2y_0)$, $(3,0,y_0^2)$	$(3,2,-2)$, $(3,1,4y_0)$, $(3,0,-2y_0^2)$	$(3,2,1)$, $(3,1,-2y_0)$, $(3,0,y_0^2)$
	$(2,2,-2x_0)$, $(2,1,4x_0{y}_0)$, $(2,0,-2x_0{y}_0^2)$	$(2,2,A_x)$, $(2,1,-2y_0{A}_x)$, $(2,0,A_x{y}_0^2)$	$(2,2,-2x_3)$, $(2,1,4x_3{y}_0)$, $(2,0,-2x_3{y}_0^2)$
	$(1,2,x_0^2)$, $(1,1,-2x_0^2{y}_0)$, $(1,0,x_0^2{y}_0^2)$	$(1,2,-B_x)$,$(1,1,-2B_x{y}_0)$, $(1,0,-B_x{y}_0^2)$	$(1,2,x_3^2)$, $(1,1,-2y_0{x}_3^2)$, $(1,0,x_3^2{y}_0^2)$
$(y_{j^{\prime },1},y_{j^{\prime },2})$	$(3,2,-2)$, $(3,1,A_y)$, $(3,0,-B_y)$	$(3,2,4)$, $(3,1,-2A_y)$, $(3,0,2B_y)$	$(3,2,-2)$, $(3,1,A_y)$, $(3,0,-B_y)$
	$(2,2,4x_0)$, $(2,1,-2x_0{A}_y)$, $(2,0,2x_0{B}_y)$	$(2,2,-2A_x)$, $(2,1,A_x{A}_y)$, $(2,0,-A_x{B}_y)$	$(2,2,4x_3)$, $(2,1,-2A_y{x}_3)$, $(2,0,2x_3{B}_y)$
	$(1,2,-2x_0^2)$, $(1,1,A_y{x}_0^2)$, $(1,0,-B_y{x}_0^2)$	$(1,2,2B_x)$, $(1,1,-A_y{B}_x)$, $(1,0,B_x{B}_y)$	$(1,2,-2x_3^2)$, $(1,1,A_y{x}_3^2)$, $(1,0,-B_y{x}_3^2)$
$(y_{j^{\prime },2},y_{j^{\prime },3})$	$(3,2,1)$, $(3,1,-2y_3)$, $(3,0,y_3^2)$	$(3,2,-2)$, $(3,1,4y_3)$, $(3,0,-2y_3^2)$	$(3,2,1)$, $(3,1,-2y_3)$, $(3,0,y_3^2)$
	$(2,2,-2x_0)$, $(2,1,4x_0{y}_3)$, $(2,0,-2x_0{y}_3^2)$	$(2,2,A_x)$, $(2,1,-2A_x{y}_3)$, $(2,0,A_x{y}_3^2)$	$(2,2,-2x_3)$, $(2,1,4x_3{y}_3)$, $(2,0,-2x_3{y}_3^2)$
	$(1,2,x_0^2)$, $(1,1,-2x_0^2{y}_3)$, $(1,0,x_0^2{y}_3^2)$	$(1,2,-B_x)$, $(1,1,2B_x{y}_3)$, $(1,0,-B_x{y}_3^2)$	$(1,2,x_3^2)$, $(1,1,-2x_3^2{y}_3)$, $(1,0,x_3^2{y}_3^2)$

Table 4. Table 4: This table tabulates the errors and orders of accuracy of the x 𝑥 x and radial forces for the case ( k = 1 , D 3 ) 𝑘 1 subscript 𝐷 3 (k=1,{D_{3}})

$N$	$L_x^1$	$L_x^2$	$L_x^{\infty }$	$L_R^1$	$L_R^2$	$L_R^{\infty }$
32	9.2865e-02	1.3093e-01	2.9581e-01	1.4516e-01	1.8516e-01	3.0346e-01
64	4.4525e-02	6.3186e-02	1.4563e-01	6.9853e-02	8.9359e-02	1.4666e-01
128	2.1765e-02	3.0984e-02	7.1745e-02	3.4178e-02	4.3819e-02	7.1872e-02
256	1.0757e-02	1.5336e-02	3.5542e-02	1.6895e-02	2.1688e-02	3.5558e-02
512	5.3472e-03	7.6281e-03	1.7681e-02	8.3991e-03	1.0788e-02	1.7683e-02
$N_p/N_{p+1}$	$O_x^1$	$O_x^2$	$O_x^{\infty }$	$O_R^1$	$O_R^2$	$O_R^{\infty }$
32/64	1.06	1.05	1.02	1.06	1.05	1.05
64/128	1.03	1.03	1.02	1.03	1.03	1.03
128/256	1.02	1.01	1.01	1.02	1.01	1.02
256/512	1.01	1.01	1.01	1.01	1.01	1.01

Table 5. Table 5: This table tabulates the errors and orders of accuracy of the x 𝑥 x and radial forces for the case ( k = 2 , D 3 ) 𝑘 2 subscript 𝐷 3 (k=2,{D_{3}})

$N$	$L_x^1$	$L_x^2$	$L_x^{\infty }$	$L_R^1$	$L_R^2$	$L_R^{\infty }$
32	1.7712e-02	2.0437e-02	4.5876e-02	2.7583e-02	2.8873e-02	4.6383e-02
64	4.5916e-03	5.3308e-03	1.2064e-02	7.1866e-03	7.5370e-03	1.2085e-02
128	1.1689e-03	1.3609e-03	3.0684e-03	1.8330e-03	1.9244e-03	3.0698e-03
256	2.9472e-04	3.4359e-04	7.7234e-04	4.6251e-04	4.8589e-04	7.7258e-04
512	7.3986e-05	8.6311e-05	1.9378e-04	1.1614e-04	1.2206e-04	1.9379e-04
$N_p/N_{p+1}$	$O_x^1$	$O_x^2$	$O_x^{\infty }$	$O_R^1$	$O_R^2$	$O_R^{\infty }$
32/64	1.95	1.94	1.93	1.94	1.94	1.94
64/128	1.97	1.97	1.98	1.97	1.97	1.98
128/256	1.99	1.99	1.99	1.99	1.99	1.99
256/512	1.99	1.99	1.99	1.99	1.99	2.00

Table 6. Table 6: This table tabulates the errors and orders of accuracy of the x 𝑥 x and radial forces for the case ( k = 3 , D 3 ) 𝑘 3 subscript 𝐷 3 (k=3,{D_{3}})

$N$	$L_x^1$	$L_x^2$	$L_x^{\infty }$	$L_R^1$	$L_R^2$	$L_R^{\infty }$
32	2.0387e-03	3.0324e-03	1.2152e-02	3.3411e-03	4.1530e-03	1.1409e-02
64	1.9115e-04	3.7919e-04	2.4028e-03	3.1242e-04	5.2093e-04	2.3122e-03
128	1.9873e-05	4.8220e-05	4.6035e-04	3.1787e-05	6.6308e-05	4.5209e-04
256	2.1659e-06	6.0452e-06	8.3088e-05	3.4478e-06	8.3172e-06	8.2252e-05
512	2.4457e-07	7.5529e-07	1.4890e-05	3.8775e-07	1.0381e-06	1.4826e-05
$N_p/N_{p+1}$	$O_x^1$	$O_x^2$	$O_x^{\infty }$	$O_R^1$	$O_R^2$	$O_R^{\infty }$
32/64	3.41	3.00	2.34	3.42	3.00	2.30
64/128	3.27	2.98	2.38	3.30	2.97	2.35
128/256	3.20	3.00	2.47	3.20	3.00	2.46
256/512	3.15	3.00	2.48	3.15	3.00	2.47

Table 7. Table 7: This table tabulates the errors and orders of accuracy of the x 𝑥 x and radial forces for the case ( k = 1 , D 4 ) 𝑘 1 subscript 𝐷 4 (k=1,{D_{4}})

$N$	$L_x^1$	$L_x^2$	$L_x^{\infty }$	$L_R^1$	$L_R^2$	$L_R^{\infty }$
32	8.1599e-02	1.2740e-01	3.3661e-01	1.2736e-01	1.8017e-01	3.4076e-01
64	3.8941e-02	6.1282e-02	1.6416e-01	6.1065e-02	8.6666e-02	1.6471e-01
128	1.9026e-02	2.9995e-02	8.0277e-02	2.9873e-02	4.2419e-02	8.0450e-02
256	9.4065e-03	1.4831e-02	3.9722e-02	1.4774e-02	2.0974e-02	3.9740e-02
512	4.6772e-03	7.3733e-03	1.9744e-02	7.3467e-03	1.0427e-02	1.9746e-02
$N_p/N_{p+1}$	$O_x^1$	$O_x^2$	$O_x^{\infty }$	$O_R^1$	$O_R^2$	$O_R^{\infty }$
32/64	1.07	1.06	1.04	1.06	1.06	1.05
64/128	1.03	1.03	1.03	1.03	1.03	1.03
128/256	1.02	1.02	1.02	1.02	1.02	1.02
256/512	1.01	1.01	1.01	1.01	1.01	1.01

Table 8. Table 8: This table tabulates the errors and orders of accuracy of the x 𝑥 x and radial forces for the case ( k = 2 , D 4 ) 𝑘 2 subscript 𝐷 4 (k=2,{D_{4}})

$N$	$L_x^1$	$L_x^2$	$L_x^{\infty }$	$L_R^1$	$L_R^2$	$L_R^{\infty }$
32	1.5794e-02	2.1363e-02	5.8539e-02	2.4639e-02	3.0197e-02	5.9092e-02
64	4.1209e-03	5.5802e-03	1.5555e-02	6.4603e-03	7.8906e-03	1.5600e-02
128	1.0494e-03	1.4225e-03	3.9659e-03	1.6471e-03	2.0116e-03	3.9692e-03
256	2.6472e-04	3.5886e-04	1.0002e-03	4.1563e-04	5.0750e-04	1.0004e-03
512	6.6463e-05	9.0107e-05	2.5103e-04	1.0436e-04	1.2743e-04	2.5105e-04
$N_p/N_{p+1}$	$O_x^1$	$O_x^2$	$O_x^{\infty }$	$O_R^1$	$O_R^2$	$O_R^{\infty }$
32/64	1.94	1.94	1.91	1.93	1.94	1.92
64/128	1.97	1.97	1.97	1.97	1.97	1.97
128/256	1.99	1.99	1.99	1.99	1.99	1.99
256/512	1.99	1.99	1.99	1.99	1.99	1.99

Table 9. Table 9: This table tabulates the errors and orders of accuracy of the x 𝑥 x and radial forces for the case ( k = 3 , D 4 ) 𝑘 3 subscript 𝐷 4 (k=3,{D_{4}})

$N$	$L_x^1$	$L_x^2$	$L_x^{\infty }$	$L_R^1$	$L_R^2$	$L_R^{\infty }$
32	1.2650e-03	1.6606e-03	4.9032e-03	1.9866e-03	2.2916e-03	4.8862e-03
64	1.2213e-04	1.6468e-04	4.4374e-04	1.8915e-04	2.2746e-04	4.3164e-04
128	1.3308e-05	1.8216e-05	4.4007e-05	2.0463e-05	2.5149e-05	4.4150e-05
256	1.5530e-06	2.1442e-06	4.9395e-06	2.3780e-06	2.9586e-06	5.3249e-06
512	1.8790e-07	2.6036e-07	5.9366e-07	2.8684e-07	3.5910e-07	6.5341e-07
$N_p/N_{p+1}$	$O_x^1$	$O_x^2$	$O_x^{\infty }$	$O_R^1$	$O_R^2$	$O_R^{\infty }$
32/64	3.37	3.33	3.47	3.39	3.33	3.50
64/128	3.20	3.18	3.33	3.21	3.18	3.29
128/256	3.10	3.09	3.16	3.11	3.09	3.05
256/512	3.05	3.04	3.06	3.05	3.04	3.03

Table 10. Table 10: This table shows the kernel integrals for α + β = 0 , 1 , 2 , 3 , 4 𝛼 𝛽 0 1 2 3 4 \alpha+\beta=0,1,2,3,4 .

$\alpha$	$\beta$	$F(x,y,\alpha ,\beta )$
0	0	$-\frac{\sqrt{x^2+y^2}}{xy}+C$
1	0	$-\ln \left(y+\sqrt{x^2+y^2}\right)+C$
0	1	$F(y,x,1,0)$
2	0	$y\ln \left(x+\sqrt{x^2+y^2}\right)+C$
1	1	$-\sqrt{x^2+y^2}+C$
0	2	$F(y,x,2,0)$
3	0	$y\sqrt{x^2+y^2}+C$
2	1	$-\frac{1}{2}x\sqrt{x^2+y^2}+\frac{1}{2}{y}^2\ln \left(x+\sqrt{x^2+y^2}\right)+C$
1	2	$F(y,x,2,1)$
0	3	$F(y,x,3,0)$
4	0	$\frac{1}{2}xy\sqrt{x^2+y^2}-\frac{1}{2}{y}^3\ln \left(x+\sqrt{x^2+y^2}\right)+C$
3	1	$-\frac{1}{3}\sqrt{x^2+y^2}\left(x^2-2y^2\right)+C$
2	2	$-\frac{1}{3}xy\sqrt{x^2+y^2}+\frac{1}{3}{y}^3\ln \left(x+\sqrt{x^2+y^2}\right)+\frac{1}{3}{x}^3\ln \left(y+\sqrt{x^2+y^2}\right)+C$
1	3	$F(y,x,3,1)$
0	4	$F(y,x,4,0)$

Equations

$$F^x(x,y,z)=-\frac{\partial \Phi }{\partial x}(x,y,z)=\int \int \int \frac{\partial \mathcal{K}}{\partial x}(x-x^{\prime },y-y^{\prime },z-z^{\prime })\rho (x^{\prime },y^{\prime },z^{\prime })d{x}^{\prime }d{y}^{\prime }d{z}^{\prime },$$

$$\nabla \cdot (\nabla \Phi )(x,y,z)=4\pi G\rho (x,y,z),$$

$$\rho (x,y,z)=\sigma (x,y)\delta (z),$$

$$F^x(x,y,z)=\int \int \int \frac{\partial \mathcal{K}}{\partial x}(x-x^{\prime },y-y^{\prime },z-z^{\prime })\sigma (x^{\prime },y^{\prime })\delta (z^{\prime })d{x}^{\prime }d{y}^{\prime }d{z}^{\prime },$$

$$F^x(x,y,0)=\int \int \frac{\partial \mathcal{K}}{\partial x}(x-x^{\prime },y-y^{\prime },0)\sigma (x^{\prime },y^{\prime })d{x}^{\prime }d{y}^{\prime }.$$

$$\sigma (x,y)\simeq \sum _{i,j}{C}_{i,j}{\phi }_{i,j}(x,y),$$

$$F^x(x,y,0)\simeq \sum _{i^{\prime },j^{\prime }}{C}_{i^{\prime },j^{\prime }}\int \int \frac{\partial }{\partial x}\mathcal{K}(x-x^{\prime },y-y^{\prime },0){\phi }_{i^{\prime },j^{\prime }}(x^{\prime },y^{\prime })d{x}^{\prime }d{y}^{\prime }.$$

$$F_{i,j}^x=F^x(x_i,y_j,0)\simeq \sum _{i^{\prime },j^{\prime }}{C}_{i^{\prime },j^{\prime }}{w}_{i-i^{\prime },j-j^{\prime }},$$

$$w_{i-i^{\prime },j-j^{\prime }}=\int \int \frac{\partial }{\partial x}\mathcal{K}(x_i-x^{\prime },y_j-y^{\prime },0){\phi }_{i^{\prime },j^{\prime }}(x^{\prime },y^{\prime })d{x}^{\prime }d{y}^{\prime }.$$

$$\displaystyle B^{[1]}_{i}(x)\equiv\left\{\begin{array}[]{cc}1,&\mbox{ if }t_{i}\leq x<t_{i+1}\\ 0,&\mbox{ otherwise}\end{array}\right.,$$

$$B_i^{[k]}(x)\equiv \frac{x-t_i}{t_{i+k-1}-t_i}{B}_i^{[k-1]}(x)+\frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}{B}_{i+1}^{[k-1]}(x).$$

$$\displaystyle B^{[2]}_{i}(x)=\left\{\begin{array}[]{ll}\displaystyle\frac{(x-t_{i})}{(t_{i+1}-t_{i})},&\mbox{ if }t_{i}\leq x<t_{i+1}\\ \displaystyle\frac{(t_{i+2}-x)}{(t_{i+2}-t_{i+1})},&\mbox{ if }t_{i+1}\leq x<t_{i+2}\\ \displaystyle 0,&\mbox{otherwise}\end{array}\right.,$$

$$\displaystyle B^{[3]}_{i}(x)=\left\{\begin{array}[]{ll}\displaystyle\frac{(x-t_{i})^{2}}{(t_{i+1}-t_{i})(t_{i+2}-t_{i})},&\mbox{ if }t_{i}\leq x<t_{i+1}\\ \displaystyle\frac{(x-t_{i})(t_{i+2}-x)}{(t_{i+2}-t_{i})(t_{i+2}-t_{i+1})}+\frac{(x-t_{i+1})(t_{i+3}-x)}{(t_{i+2}-t_{i+1})(t_{i+3}-t_{i+1})},&\mbox{ if }t_{i+1}\leq x<t_{i+2}\\ \displaystyle\frac{(t_{i+3}-x)^{2}}{(t_{i+3}-t_{i+1})(t_{i+3}-t_{i+2})},&\mbox{ if }t_{i+2}\leq x<t_{i+3}\\ \displaystyle 0,&\mbox{otherwise}\end{array}\right.,$$

$$w_{i-i^{\prime },j-j^{\prime }}=\int \int \frac{\partial }{\partial x}\mathcal{K}(x-x^{\prime },y-y^{\prime },0)B_{i^{\prime }}^{[k]}(x^{\prime })B_{j^{\prime }}^{[k]}(y^{\prime })d{x}^{\prime }d{y}^{\prime },$$

$$F(x,y,\alpha ,\beta )$$

$$=$$

$$G(x,a,\alpha )=\int \frac{x^{\alpha }}{(x^2+a^2{)}^{3/2}}dx.$$

$$G(x,y,\alpha )=\frac{x^{\alpha -1}}{\alpha -2}\frac{1}{(x^2+y^2{)}^{1/2}}-\frac{\alpha -1}{\alpha -2}{y}^{2}G(x,y,\alpha -2).$$

$$\tilde{F}(x_{\ell },x_r,y_{\ell },y_r,\alpha ,\beta )$$

$$w_{i-i^{\prime },j-j^{\prime }}$$

$$\gamma \tilde{F}(x_{\ell },x_r,y_{\ell },y_r,\alpha ,\beta )$$

$$w_{i-i^{\prime },j-j^{\prime }}$$

$$\sigma (x,y)\simeq \sum _{i,j}{C}_{i,j}{B}_i^{[k]}(x)B_j^{[k]}(y).$$

$$\sigma (x_i,y_j)=\sum _{i^{\prime },j^{\prime }}{C}_{i^{\prime },j^{\prime }}{B}_{i^{\prime }}^{[k]}(x_i)B_{j^{\prime }}^{[k]}(y_j),$$

$$W=[W_{i,j}{]}_{1\le i,j\le 3}=\frac{1}{64}\left[\begin{array}{ccc}1& 6& 1\\ 6& 36& 6\\ 1& 6& 1\end{array}\right],$$

$${\sigma }_{i,j}$$

$$\left[\begin{array}{cccccc}6& 1& 0& 0& \dots & 0\\ 1& 6& 1& 0& \dots & 0\\ 0& 1& 6& 1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1& 6& 1\\ 0& 0& \dots & 0& 1& 6\end{array}\right]\left[\begin{array}{c}{z}_{i,1}\\ z_{i,2}\\ z_{i,3}\\ ⋮\\ z_{i,N-1}\\ z_{i,N}\end{array}\right]=\left[\begin{array}{c}{\sigma }_{i,1}\\ {\sigma }_{i,2}\\ {\sigma }_{i,3}\\ ⋮\\ {\sigma }_{i,N-1}\\ {\sigma }_{i,N}\end{array}\right],$$

$$\frac{1}{64}\left[\begin{array}{cccccc}6& 1& 0& 0& \dots & 0\\ 1& 6& 1& 0& \dots & 0\\ 0& 1& 6& 1& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & 1& 6& 1\\ 0& 0& \dots & 0& 1& 6\end{array}\right]\left[\begin{array}{c}{C}_{1,j}\\ C_{2,j}\\ C_{3,j}\\ ⋮\\ C_{M-1,j}\\ C_{M,j}\end{array}\right]=\left[\begin{array}{c}{z}_{1,j}\\ z_{2,j}\\ z_{3,j}\\ ⋮\\ z_{M-1,j}\\ z_{M,j}\end{array}\right],$$

$${\sigma }_{D_n}(R;\tilde{\alpha })=\{\begin{array}{ll}{\sigma }_0{\left(1-\frac{R^2}{{\tilde{\alpha }}^2}\right)}^{n-1/2},& \text{ for }R<\tilde{\alpha }\\ 0,& \text{ for }R\ge \tilde{\alpha },\end{array},$$

Full text

gravitation - methods: numerical

††journal: ApJ

1 Introduction

Evolution of thin disks is of great interest in the context of galactic dynamics and protoplanetary disks formed around young stars. The formation of massive stars, the substructures, such as spurs and feathers, associated with barred and spirals galaxies suggest that the self-gravity of gaseous disks plays a role in the evolution of galactic disks (Kim et al., 2012; Lee & Shu, 2012; Lin et al., 2013; Lee, 2014; Elmegreen et al., 2014; Seo & Kim, 2014). Furthermore, the existence of starburst rings located at the galactic centers also indicates that the rings themselves are gravitationally unstable (Hsieh et al., 2011; Lin et al., 2013; Seo & Kim, 2014; Krumholz & Kruijssen, 2015; Krumholz et al., 2017; Seo et al., 2019). The self-gravity may be also important to planet formation in protoplanetary disks and was suggested to affect the planet radial migration (Inutsuka et al., 2010; Zhang et al., 2008, 2014).

Although high-order numerical schemes, such as Essentially Non-Oscillatory (ENO; Harten et al. (1987)), Weighted ENO (WENO; Liu et al. (1994); Jiang & Shu (1996) ) and Piecewise Parabolic Method (PPM; Colella & Woodward (1984)), have been developed for hyperbolic systems (for a recent comprehensive review, see Balsara (2017)), a fast and high-order scheme for self-gravitational force calculation for thin disks is relatively less explored. A direct approach of 2nd-order of accuracy for self-gravitational calculation for infinitesimally thin gaseous disks has been proposed by Yen et al. (2012). Within the framework, the approach has been extended to polar coordinates by Wang et al. (2015), and adapted to mesh refinement in Wang et al. (2016). For gravitational forces in $x$-direction, $F^{x}$,:

[TABLE]

where $\displaystyle{\cal K}=\frac{1}{\sqrt{x^{2}+y^{2}+z^{2}}}$ is the three-dimensional kernel and the potential satisfies the Poisson equation:

[TABLE]

where $G$ is the gravitational constant and $\rho$ is the volume density. Relations similar to Eq. (1) for $F^{x}$ can be found for $F^{y}$ and $F^{z}$. The “direct” approach has two meanings: one is to calculate the gravitational forces directly, instead of taking negative gradient from a pre-calculated gravitational potential, $\Phi(x,y,z)$. The other meaning is to compute the integral in Eq. (1) directly and reach high-order of accuracy.

The surface density, $\sigma(x,y)$, of an infinitesimally thin disk relates to the volume density through:

[TABLE]

where $\delta$ represents the Dirac delta function. Substituting Eq. (2) into Eq. (1) formally:

[TABLE]

and integrating over $z^{\prime}$, the expression of $F^{x}$ in the plane $z=0$ reads:

[TABLE]

It is important to realize that the singularity involved in the force evaluation, Eq. (4), is removable and integrable. Moreover, Eq. (4) can be theoretically obtained by applying the Lebesgue’s dominated convergence theorem known in measure theory. That is, an infinitesimally thin disk can be approached by a sequence of models with diminishing disk thickness. In Yen et al. (2012), the 2nd-order accuracy is achieved by approximating the a smooth density distribution using Taylor expansion to the linear term in a local sense. The convolution form in Eq. (4) enables the use of Fast Fourier Transform (FFT) that keeps the numerical complexity linear, i.e., O(N^{2}{\color[rgb]{0,0,0}\ln^{2}}N), where $N$ is the number of zones in one direction.

The objective of this work is to utilize the spline functions to improve the order of accuracy for force calculations in a general way. This approach involves two parts. The first is to approximate the surface density to high-order accuracy using a set of linearly independent splines, {$\varphi_{i,j}$}, where the subscripts $i,j$ are used to distinguish between spline functions. The use of splines for interpolation is well developed in the literature, a rigorous proof and in-depth discussion is beyond the scope of this work. Interested readers are referred to textbooks on numerical analysis (Prenter, 1975; de Boor, 1978; Epperson, 2013). The second part involves the calculation in terms of the discretized convolution form. For the first part, the surface density can be approximated by:

[TABLE]

where $C_{i,j}$ denotes the coefficient associated with the spline $\varphi_{i,j}(x,y)$. Combining Eqs. (4) and (5), we have:

[TABLE]

For fast calculation in the second part, the locations, $(x_{i},y_{j})$, for evaluating forces are chosen such that Eq. (6) can be recast in a convolution form:

[TABLE]

where

[TABLE]

The paper is organized as follows. In § 2, we detail the calculation of Eq. (8) for basis spline functions ($B$-splines, hereafter). Generalization to high-order accuracy is presented. In § 3, we implement the proposed algorithm up to 3rd-order of accuracy and conduct a series of numerical tests. The numerical results are compared with analytic solutions of generalized Maclaurin disks (Schulz, 2009; Wang et al., 2015). Finally, we summarize the works and discuss some theoretical issues in § 4.

2 A direct method using $B$-splines

In this work, we adopt the $B$-splines as the basis functions to interpolate surface density. The first order of $B$-splines in one dimension reads:

[TABLE]

where the upper script represents the order of $B$-splines and the set $\{t_{i}\}$ (a.k.a knots) denotes the locations where different $B$-splines meet. Conventionally, $\{t_{i}\}$ is sorted in a non-decreasing order. The $B$-spline of order $k>1$ is defined recursively:

[TABLE]

Based on the definition of $B$-splines, the support of $B^{[k]}_{i}(x)$ is $[t_{i},t_{i+k})$. For $k=2$ and $k=3$, the explicit expressions are:

[TABLE]

and

[TABLE]

respectively. In this work, the values of $t_{i}$ are all distinct, thus the first $k-1$ derivatives of $B$-splines are continuous across $t_{i}$. It is the finite support and the smoothness properties that make $B$-splines a good choice for approximating the surface density in Eq. (5). The two-dimensional basis set is formed using the tensor product $\varphi_{i,j}(x,y)=B^{[k]}_{i}(x)B^{[k]}_{j}(y)$. Substituting this expression of $\varphi_{i,j}(x,y)$ into Eq. (8), it follows:

[TABLE]

for $(x,y)=(x_{i},y_{j})$. Suppose that the surface density is confined within a finite rectangular computational domain ${\cal D}=[L,R]\times[B,T]$, where $L<R$ and $B<T$, and the domain is uniformly subdivided by rectangular cells, with $\Delta x$ and $\Delta y$ being the cell sizes in $x$- and $y$-direction, respectively. The grid points are $x_{i+1/2}=L+i\Delta x$ for $i=0,1,...,M$ and $y_{j+1/2}=B+j\Delta y$ for $j=0,1,...,N$. The cell center at $(i,j)$ is $\displaystyle(x_{i},y_{j})=(\frac{x_{i-1/2}+x_{i+1/2}}{2},\frac{y_{j-1/2}+y_{j+1/2}}{2})$. Figure 1a shows a spline $B^{[3]}(x)$ peaked at the cell center $x_{i}$, with $x_{i-1/2}$ and $x_{i+1/2}$ being the cell boundaries. For this particular $B$-spline of $k=3$, the size of its support is $3\Delta x$ and its first and second derivatives are continuous over the whole computational domain. Figure 1b shows the tensor product $B^{[3]}(x)B^{[3]}(y)$ peaked at $(x_{i},y_{j})$, with white dashed-lines indicating the cell boundaries located at half-integer subscripts. The map of this tensor are color-coded and the values shown in the figure are evaluated for the cell center $(x_{i},y_{j})$ (black- font number) and for those centers surrounding it (white- font numbers).

2.1 Calculation of Equation (22)

As will be shown in § 2.2, the calculation of Eq. (22) involves an integration of the form:

[TABLE]

where $\alpha$ and $\beta$ are non-negative integers (see Appendix B for the derivation of Eq. (24)) and

[TABLE]

Equation (26) can be evaluated through the following recursion relation:

[TABLE]

The derivation of Eq. (27) is given in Appendix B. Although Eq. (24) appears to apply only when $\alpha+\beta>1$, this is not a restrictive limitation since the close forms for $\alpha+\beta=0,1,2,3,4$ exist and are listed in Appendix A222Symbolic mathematical software, Maple, fails to find the close form for $(\alpha,\beta)=(4,0),(2,2)$ and $(0,4)$. For the convenience of discussion below, we define $\tilde{F}(x_{\ell},x_{r},y_{\ell},y_{r},\alpha,\beta)$ to be the definite integral of $F(x,y,\alpha,\beta)$ over the domain $[x_{\ell},x_{r}]\times[y_{\ell},y_{r}]$:

[TABLE]

2.2 Examples

Case 1. ($k=1$)

For $(x,y)=(x_{i},y_{j})$ and from Eq. (11), Eq. (22) reads:

[TABLE]

where in Eq. (29) we have applied the change of variable $(\tilde{x},\tilde{y})=(x^{\prime}-x,y^{\prime}-y)$ and introduced auxiliary variables $x_{i^{\prime},r}=x_{i^{\prime}-1/2+r}-x$ and $y_{j^{\prime},r}=y_{j^{\prime}-1/2+r}-y$ for $r=0,..,k$ to denote the ranges of integration. Eq. (29) is simply $\tilde{F}(x_{i^{\prime},0},x_{i^{\prime},1},y_{j^{\prime},0},y_{j^{\prime},1},1,0)$. In Table 1 the integration Eq. (29) is presented in a more concise way. Each entry in the table denotes a triplet $(\alpha,\beta,\gamma)$, where $\gamma$ is the associated coefficient.

In general, as will be shown and become more clear in the following cases, if one attempts to expand the integral $w_{i-i^{\prime},j-j^{\prime}}$ term by term, each term would end up with a form:

[TABLE]

Each term has a different range of integration (see the definition of $B$-splines but generalize it for a tensor in 2D), different $(\alpha,\beta)$ and a constant coefficient $\gamma$ associated with it. Take $k=3$ for example, one would expect nine different integration ranges (three in $x$ and three in $y$, respectively; see also the nine blocks shown in Fig. 1b), which corresponds the size of Table 3. Within one integration range, one would expect nine terms associated with different $(\alpha,\beta)$ and a coefficient $\gamma$ attached to each of them. In the following cases, to avoid the clutter but straightforward formulas, we summarize that information in a table for each case.

Case 2. ($k=2$)

For $(x,y)=(x_{i},y_{j})$ and from Eq. (16), Eq. (22) reads:

[TABLE]

The associated integrations are listed in Table 2. We note that the final result would need to sum over all the integration results listed in Table 2 and multiply by a factor $1/(\Delta x\Delta y)$.

Case 3. ($k=3$) With the recursion relation Eq. (24), it is straightforward to generalize to a third-order ($k=3$) of accuracy or higher. For $(x,y)=(x_{i},y_{j})$ and from Eq. (21), Eq. (22) reads:

[TABLE]

The calculation for $k=3$ involves 81 integrations as listed in Table 3. We emphasize that, the complicated evaluation of Eq. (22) need only to be done once at the beginning of simulations, and the results can be reused afterwards when needed. We note that the final result would need to sum over all the integration results listed in Table 3 and multiply by a factor $1/(4(\Delta x)^{2}(\Delta y)^{2})$.

2.3 Solving a linear system of $C_{i,j}$

Using the $B$-splines as the basis functions, Eq. (5) now reads:

[TABLE]

That is, we are using $B$-splines to interpolate the surface density. By construction, a $B$-spline has a finite support and a symmetric point. Except for $k=1$ (piecewise constant), the symmetric center corresponds to the location of maximum. We found that it is advantageous to have the symmetric centers of $B$-splines coincide with cell centers, $(x_{i},y_{j})$. Applying the condition:

[TABLE]

we are solving an unknown matrix $C_{i,j}$ with a size the same as $\sigma_{i,j}\equiv\sigma(x_{i},y_{j})$. Due to the finite support of $B$-splines, unless $k$ is large, usually only those $B$-splines in the vicinity of $(x_{i},y_{j})$ are involved. For example, when $k=1,2$, one may simply take $C_{i,j}=\sigma(x_{i},y_{j})$ given that $B^{[k]}_{i}(x)B^{[k]}_{j}(y)$ is normalized. For higher order $B$-splines, we adopt $k=3$ to detail the procedure as follows. We define a weight matrix:

[TABLE]

where the factor in front of the matrix is used to normalize the matrix such that the summation over all the entries of W is unity. The entries in $W$ describe the weight distribution of a $B$-spline tensor $B^{[3]}_{i}(x)B^{[3]}_{j}(y)$ for the cell centers located within its support. A visual interpretation of the weight matrix is shown in Fig. 1b, which describes how a $B^{[3]}(x)B^{[3]}(y)$ tensor peaked at $(x_{i},y_{j})$ contributes to itself and those cell centers adjacent to it. If a set of $B$-spline tensors is used to interpolate a density located at $(x_{i},y_{j})$, it is obvious that only those nine tensors surrounding $(x_{i},y_{j})$ will contribute. The finite support of $B$-spline therefore greatly simplifies the linear system for density interpolation. For a linear system of $k=3$, the surface density $\sigma_{i,j}$ has the following relation with $C_{i,j}$:

[TABLE]

where $z_{i,j}\equiv\frac{1}{64}C_{i-1,j}+\frac{6}{64}C_{i,j}+\frac{1}{64}C_{i+1,j}$.

This suggests to solve a symmetric tridiagonal linear system for $z_{i,j}$:

[TABLE]

where $i=1,...,M$. One can efficiently solve this matrix using Gaussian elimination. Once the matrix $z_{i,j}$ is solved, we may proceed to solve the a matrix of the same type for $C_{i,j}$:

[TABLE]

where $j=1,...,N$. From Eq. (36) to Eqs. (37)(38), we have applied the isolated boundary condition $C_{i,j}=0$ when $i<1$ or $j<1$ or $i>M$ or $j>N$. We note that although the matrix presentation in Eq. (37) appear to solve for single column $i$, in practice, we may solve for all columns of $z_{i,j}$ independently and simultaneously using Message Passing Interface (MPI) if the density matrix is distributed along one-dimension as the data arrangement needed for using the package FFTW. After the matrix $z_{i,j}$ is solved, a transposing operation is needed before solving the system Eq. (38). This can be done efficiently using the transposing functions provided by FFTW. We emphasize that the step of Gaussian elimination is fast, and the associated numerical complexity is linear. This procedure for solving $C_{i,j}$ is systematic and can be easily extended to a higher order method with $k>3$.

3 simulations and Results

In this section, we implement the algorithm detailed in § 2 for $k=1,2,3$ and compare the results with analytic solutions of the generalized Maclaurin disks, $D_{3}$ and $D_{4}$ (Wang et al., 2015). We will show that the methods implemented for the $B$-spline of order $k=1,2,3$ correspond to the 1st-, 2nd- and 3rd-order of accuracies in terms of the one, two, and infinite norms.

Now, we have all the ingredients needed for the force calculations, i.e., Eq. (7). To validate the methods, we need analytic density-force pairs that are sufficiently smooth for testing the order of accuracy. In addition, for exploring the order of accuracy, the mass distribution of disk models has to be finite to be entirely embedded in a computational domain. To our knowledge, the family of generalized Maclaurin disks satisfies both requirements.

The surface density of a generalized Maclaurin disk has the form:

[TABLE]

where $R$ is the polar radius, $\sigma_{0}$ is the disk central density, $\tilde{\alpha}$ represents the disk size and $n=0,1,2,...$ the order of disk. The analytic expressions of gravitational forces for $n=0,1,2$ disks are reported by Schulz (2009) and generalized to arbitarly positive integers by Wang et al. (2015, 2016). The analytic forms of radial forces for $D_{3}$ and $D_{4}$ disks in the $z=0$ plane are:

[TABLE]

and

[TABLE]

repectively, and $\eta\equiv\tilde{\alpha}/R$ is a dimensionless variable. These analytic expressions are derived using the recursion formulas Eqs.(53)-(55) in Wang et al. (2015). It will be shown below that the smoothness of a generalized Maclaurin disk is important for tests. In general, the $n$th derivative of $\sigma_{D_{n}}$ will encounter singularities at the edge of the disk. Thus a $D_{n}$ disk is not suitable for the test of an algorithm of $n$th-order of accuracy. That is, to test an algorithm of $n$th-order accuracy, one needs at least a disk of $(n+1)$th order.

Without loss of generality, we set the computational domain $\Omega=[-1,1]\times[-1,1]$, the parameters for the disks $\sigma_{0}=G=1$ and the disk size $\tilde{\alpha}=0.5$. The $\Omega$ is uniformly subdivided by $N\times N$ square cells, with $N=32,64,128,256,512$ being the numbers of zones in both $x$ and $y$-directions. The errors compared to the analytic solutions are measured in terms of the one, two and infinite norms:

[TABLE]

where $L^{1}$, $L^{2}$, $L^{\infty}$ are the one norm, two norm and infinite norm of error, and $F^{{\rm num}}_{i,j}$, $F^{{\rm ana}}_{i,j}$ are numerical and exact forces at locations indexed by $(i,j)$, respectively. When using $L^{1}$ and $L^{2}$, we evaluate the total variation and energy in a global sense, while using $L^{\infty}$ we monitor the convergence of maximum errors in a pointwise sense.

3.1 Results

Six cases are considered to study the order of accuracy of the algorithm. Tables 4 to 9 show the numerical results of $(k=1,{D_{3}})$, $(k=2,{D_{3}})$, $(k=3,{D_{3}})$, $(k=1,{D_{4}})$, $(k=2,{D_{4}})$ and $(k=3,{D_{4}})$, respectively. In each case, errors of $x$-forces and radial forces are measured using the one, two and infinity norms, as denoted by $L_{x}^{1}$, $L_{x}^{2}$, $L_{x}^{\infty}$, $L_{R}^{1}$, $L_{R}^{2}$, $L_{R}^{\infty}$, respectively. The upper tables tabulate the numerical values of error measurements with increasing $N$, while the lower tables tabulate the orders of accuracy between different spatial resolutions.

The results of the case $(k=1,D_{3})$ are shown in Table 4. The one norm errors of $x$-forces ($L^{1}_{x}$) are listed in the second column in the order of increasing cell numbers, $N$. It is obvious that the numerical accuracy is improved with increasing $N$, i.e., increasing spatial resolution. The degree of improvement for $x$-forces when doubling the spatial resolution is calculated using $O^{1}_{x}=\log_{2}(L^{1}_{x,N_{p}}/L^{1}_{x,N_{p+1}})$, where the subscript $p=1,2,3,4$ is used to indicate the comparisons between different cell numbers. For instance, when $(N_{p}/N_{p+1})=(32/64)$, the $O^{1}_{x}=\log_{2}(\mbox{9.2865e-02/4.4525e-02})=1.06$ measures the degree of improvement, i.e., the order of accuracy of a method using one norm. The order of accuracy depends on how we measure the errors. The most rigorous measurement is the accuracy improvement for the maximum errors, as denoted by $O^{\infty}_{x}$ in the fourth column. From the second to the fourth columns, we measure the numerical errors using different error norms and list the corresponding order of accuracy for $x$-forces, while from the fifth to the seventh columns, the corresponding measurements are listed for the radial forces. Similar interpretations apply to Tables 5 to 9.

From Table 4 $(k=1,D_{3})$, we observe that the order of accuracy converges steadily toward one in terms of one, two and infinite norms. As expected, the result shows that the method using $B$-splines of $k=1$ is of 1st-order accuracy. Similarly, as shown in Table 5$(k=2,D_{3})$, the method using $B$-splines of $k=2$ is indeed of 2nd-order accuracy for all error norms. However, in Tabel 6$(k=3,D_{3})$, although the method using $B$-splines of $k=3$ converges to 3rd-order accucracy with the one and two norms, $O^{\infty}_{x}$ and $O^{\infty}_{R}$ degrade to be less than 2.5. The reason is that the ${D_{3}}$ disk is not sufficiently smooth at the edge of the disk. This is clearly seen in the Fig. 2a, where the maximum errors of the radial forces coincide with the edge of the $D_{3}$ disk. This also indicates that we need to test the method using $B$-splines of $k=3$ with a smoother disk. We, therefore, adopt the $D_{4}$ disk for the tests and the results of $k=1,2,3$ are shown in Tables 7 to 9, respectively.

From Tables 7 $(k=1,D_{4})$ and 8$(k=2,D_{4})$, the orders of accuracy converge to one and two, respectively, for all error norms. For $k=1,2$ the $D_{4}$ disk behaves the same as the $D_{3}$ disk. However, as shown in Table 9 for the case $(k=3,D_{4})$, all the error norms are of 3rd-order. If we observe the error distribution for the case $(k=3,D_{4})$ shown in Fig. 2b, the errors converge more uniformly over the disk. Similar conclusion was also drawn for the $D_{2}$ disk in Wang et al. (2015), where they studied an algorithm of 2nd-order accuracy in polar coordinates.

For a given disk model, one may compare the results between different orders of $B$-splines. We observe that the improvement in numerical accuracy is significant with increasing order of $B$-spline. For example, for the cases with $N=512$, the improvement is more than one order of magnitude between $k=1$ and $k=2$ for both ${D_{3}}$ and ${D_{4}}$ disks. The improvement can be comfortably more than two orders of magnitude between $k=2$ and $k=3$ for the ${D_{4}}$ disk for all error norms. Comparing the gain from accuracy improvement with increasing complexity of the algorithm, it is much more effective to go for an algorithm of higher order.

4 Summary and Discussions

4.1 Summary

In this work, we developed a new algorithm that is fast and of high-order accuracy for calculating self-gravity of infinitesimally thin disks. The method is fast since the force calculations can be represented by a convolution form, i.e., Eq. (7). Based on the convolution theorem, a convolution associated with two functions can be evaluated as a pointwise product of their Fourier transforms. Thus, the calculation of convolution can be accelerated by the use of FFT, which is of linear complexity. The high-order accuracy is achieved by approximating the surface density using high-order $B$-splines. The use of $B$-splines has several advantages. First, $B$-splines of high-order are defined recursively. Thus, the mathematical form of basis functions for a high-order method is well-defined. Second, $B$-splines have finite supports. This property directly leads to a linear system with a sparse matrix described in (§ 2.3), which can be solved efficiently with linear numerical complexity. Third, the first $k-1$ derivatives of a $B$-spline of order $k$ are continuous. The increasing smoothness with $k$ makes $B$-splines an ideal basis for expanding a smooth function.

Although the concepts are straightforward, generalizing to high-order of accuracy for practical uses is not obvious. For the calculation of the kernel integral Eq. (22), we derived the recursion relations Eqs. (24) and (27). In § 2.2, detailed calculations are presented for $B$-spline tensors of order $k=1,2,3$. By properly arranging the locations of $B$-spline tensors, we have shown that for $k=1,2$, the calculation of $C_{i,j}=\rho(x_{i},y_{j})$ is easy, while for $k=3$, a linear system needs to be solved. The linear system is symmetric and sparse with respect to the diagonal band. We solve the system using Gaussian elimination, which is also of linear complexity. All the numerical steps involved in the force calculations are of linear complexity. It concludes that the overall complexity of the algorithm is therefore linear. The algorithm can be deployed on a distributed memory system and benefit the parallel computation using MPI.

To validate the algorithm, we implemented the methods and calculated the self-gravitational forces for the generalized Maclaurin disks $D_{3}$ and $D_{4}$. The numerical results are compared with the analytic solutions discussed in § 3. We adopted the one, two and infinite norms to quantify the order of accuracy. We identified that a $D_{3}$ disk is not sufficiently smooth for validating a method of 3rd-order. By adopting a $D_{4}$ disk, we have proved that the method using the $B$-spline tensors of $k=3$ is indeed of 3rd-order accuracy.

4.2 Discussions

The use of $B$-splines for surface density interpolation has more advantages compared to the Taylor expansion used by Yen et al. (2012). First, it reduces the number of performing FFTs. In Yen et al. (2012), three times of forward FFTs are needed for both $x-$ and $y$-direction, respectively, and two inverse FFTs are required to get back the $x$ and $y$ forces. That is, for the 2nd-order accuracy, five FFTs are required before we obtain the self-gravitational forces. For the 3rd-order accuracy using Taylor expansion, eight FFTs need to be performed for every force calculation. However, in this work, two Gaussian eliminations (not required for $k=1,2$), one forward FFT for $C_{i,j}$ and two inverse FFTs for getting the forces are required for all orders of accuracy. Most of the heavy workload is shifted to the kernel calculations, which needs to be done only once at the beginning of simulations and the results can be reused afterward when needed.

Second, if one observes the behavior of the order of accuracy for all tables, and compares it with those tables presented by Yen et al. (2012), it is obvious that the algorithm using $B$-spline for interpolation is more stable than the algorithm adopting Taylor expansion. Instead of oscillating in the values of order of accuracy as in the work of Yen et al. (2012), the method using $B$-splines converges to the expected order smoothly and steadily. In reality, if a disk has discontinuity in the density, the interpolation would introduce the so-called Gibbs phenomenon in the vicinity of discontinuities. However, since the force calculations are related to the surface density through integrations, this side effect can be ignored.

As shown in the Fig. 2, adopting an appropriate model that is sufficiently smooth is important for validating a method of a given order of accuracy. For example, for the $B$-splines of $k=2$, a $D_{3}$ disk is sufficiently smooth for validating a method of 2nd-order accuracy in terms of the infinite norm. However, for $(k=3,D_{3})$, the order of accuracy degrades in terms of infinite norm since the edge of the disk is not sufficiently smooth to the third-order. This might not be clearly observed with $L^{1}$ and $L^{2}$ since the edge is an one-dimensional structure while these norms evaluate the average values of errors over the whole calculation domain. The squareness of errors seen in Fig. 2b is due to the squareness of cells and the squareness of the $B$-spline tensors, which can be also observed in Fig 1b. Since our testing cases are axisymmetric and our coordinates are Cartesian, the conclusions drawn are generally applicable to more complicated density distribution.

The authors acknowledge the support of the Theoretical Institute for Advanced Research in Astrophysics (TIARA) based in Academia Sinica Institute of Astronomy and Astrophysics (ASIAA). C.C.Y. thanks the Institute of Astronomy and Astrophysics, Academia Sinica, Taiwan for their constant support. C.C.Y. is supported by Ministry of Science and Technology of Taiwan, under grant MOST-107-2115-M- 030-005-MY2. HHW thanks the support by the grant from the Research Grants Council of Hong Kong: General Research Fund 14308217, and the two supports by the Research Committee Direct Grant for Research from CUHK: 4053229, 4053309. M.-C. Shiue was supported in part by the Grant MOST-106-2115-M-009-011-MY2.

[FIGURE:]

Appendix A Formula

For the integral $\displaystyle F(x,y,\alpha,\beta)=\int\!\!\!\int\frac{x^{\alpha}y^{\beta}}{(x^{2}+y^{2})^{3/2}}\,dxdy$ with $\alpha$ and $\beta$, some of its antiderivatives are listed below.

Appendix B Derivation of Equations (15)(18)

Before deriving Eq. (15), we shall derive Eq. (18) first.

[TABLE]

where Eq. (B1) is obtained using the integration by part and the integral in the third term of Eq. (B2) is the same as the original integral. From Eq. (B2), we obtain a recursion relation, Eq. (18):

[TABLE]

Now, we are in a position to derive Eq. (15). To begin with,

[TABLE]

where Eq. (B4) is obtained by multiplying Eq. (B1) with $y^{{\beta-2}}$ and integrating the resulting in the third term in Eq. (B5) is the same as the original one and the integral in the fourth term in Eq. (B5) is symmetric to the original one, it can be inferred from Eq. (B5) that

[TABLE]

By the same token, we obtain a relation for the symmetric expression:

[TABLE]

Substituting Eq. (B6) into Eq. (B7), we have:

[TABLE]

Substituting Eq. (B7) into Eq. (B6), we have another expression:

[TABLE]

An expression that reflects the symmetry in $x$ and $y$ is obtained by $F(x,y,\alpha,\beta)=(\mbox{B8+B9})/2$:

[TABLE]

Since the close form of the integral Eq. (B3) exists for $\alpha=0,1,2$, the evaluation of Eq. (B10) is straightforward. As a reference, the analytic expressions for $\alpha=0,1,2$ are listed below:

[TABLE]

where $C$ is an integration constant. For example, to evaluate $F(x,y,2,1)$ from Eq. (B10):

[TABLE]

From Eq. (B14) to Eq. (B15), we have used Eqs. (B12) and (B13). To illustrate the use of Eq. (B3), we evaluate $F(x,y,4,0)$ from (B10) as follows:

[TABLE]

From Eq. (B17) to Eq. (B18), we have applied Eq. (B3) to evaluate the first integral and Eq. (B11) to the second. From Eq. (B18) to (B19), we have applied (B13) to the integral.

For a general $F(x,y,\alpha,\beta)$, the evaluation procedure is the following:

Apply Eq. (B10). 2. 2.

If the exponent of the integrand numerator is larger than 2, apply Eq. (B3) to reduce the exponent by 2. 3. 3.

If the exponent of the integrand numerator is equal or less than 2, apply Eqs. (B11) to (B13) directly; go back to step 2, otherwise.