 grids (Section \[sec:cart\_int\]) and cylindrical grids (Section \[s:interior\_solver\_cylindrical\]). We note that in principle, alternative fast and efficient solvers (such as multigrid) could be employed for computing the interior potential given Dirichlet boundary conditions for the potential.

The interior solver is employed in three different instances in our method. The first instance is our use of the interior solver to obtain the numerical DGF, as described in Appendix \[s:calc\_dgf\]; this is done once at the beginning of any simulation. The other two instances are in the first and fourth step of the James’s method; each time the Poisson solution is required, two calls to the interior solver are made.

In the first step of James’s method, the interior potential from the original density distribution is computed subject to the zero boundary condition. In the fourth step of James’s method, the interior potential must be computed subject to a known boundary potential $\Phi^{\rm B}$; a method to obtain $\Phi^{\rm B}$ will be presented in Section \[s:bc\]. Formally, we define $\Phi^{\rm B}$ as having nonzero value only in a single layer of ghost zones immediately outside the active domain. Then we can write the desired potential as $\Phi = \widetilde{\Phi} + \Phi^{\rm B}$ where the required boundary condition is $\widetilde{\Phi}=0$. From the definition of $\Phi^{\rm B}$, $\Delta^2
\Phi^{\rm B}$ will be nonzero only in the single layer of active zones adjoining the domain boundaries. We can thus define a modified density distribution $\rho \rightarrow \rho - \Delta^2
\Phi^{\rm B}/(4\pi G)$ which is the same as the original density distribution everywhere except in the layer just inside the domain boundaries. We then employ this modified density distribution following the procedure of Section \[sec:cart\_int\] or Section \[s:interior\_solver\_cylindrical\] to compute $\widetilde{\Phi}$. Within the interior, where $\Phi^{\rm B}=0$, this solution is then the desired final solution $\Phi$.

Note that in S