\rm out\to inn}) \sigma^m_{k'}({\rm out}){\cal V}_{N_R+1},\end{aligned}$$ $$\begin{aligned}
\label{eq:Theta_out}
  \Theta^m_k({\rm out}) &= \sum_{i'=1}^{N_R} {\cal G}^m_{k,i'}({\rm top\to out}) \sigma^m_{i'}({\rm top}){\cal V}_{i'} + \sum_{i'=1}^{N_R} {\cal G}^m_{k,i'}({\rm bot\to out}) \sigma^m_{i'}({\rm bot}){\cal V}_{i'}\nonumber\\
                        &+ \sum_{k'=1}^{N_z} {\cal G}^m_{k-k'}({\rm inn\to out}) \sigma^m_{k'}({\rm inn}){\cal V}_{0} + \sum_{k'=1}^{N_z} {\cal G}^m_{k-k'}({\rm out\to out}) \sigma^m_{k'}({\rm out}){\cal V}_{N_R+1},\end{aligned}$$ where we use symbolic notations such that ${\cal G}^m_{i,k'}({\rm inn\to top}) = {\cal G}^m_{i,0,N_z+1-k'}$, $\sigma^m_{i'}({\rm top}) = \sigma^m_{i',N_z+1}$, etc. Finally, we apply an inverse Fourier transform to $\Theta^m_i({\rm top}), \Theta^m_i({\rm bot}), \Theta^m_k({\rm inn})$, and $\Theta^m_k({\rm out})$ to obtain the boundary potential $\Theta_{i,j,k}^{\rm B}$ due to the surface charges. Then, the desired boundary potential $\Phi^{\rm B}$ due to the original charge $\rho$ is given by $$\label{eq:boundary_condition}
  \Phi^{\rm B}_{i,j,k} = \Psi^{\rm B}_{i,j,k} -\Theta_{i,j,k}^{\rm B} = -\Theta_{i,j,k}^{\rm B},$$ which gives the required Dirichlet boundary condition for the interior solver. Note that the boundary potential calculations explained above involve FFTs on 2D arrays (e.g., $\sigma^m_{i'}({\rm top})$) together with the summations of the Green’s function amounting to ${\cal O}(N^3)$ operations. Therefore, the overall computational cost of the boundary potential calculation is of order ${\cal O}(N^3 + N^2\log N)$, similarly to the case with a Cartesian grid.

We note that the above formulation for the boundary potential is valid for a mass distribution under $P$-fold symmetry in $\phi$. In Appendix \[s:P-fold\_symm\], we directly demonstrate that this is really the case as long as ${\cal G}_{i,i',j-j',k-k'}$ in Equation properly accounts for the contributions to the boundary potential from all periodic images of the mass density.

TEST RES