lace operator ($-k^2$).

The discrete analog of the Sturm-Liouville theory (e.g., @hil68 [chap. 1.10–1.16]; see also, @atk64) guarantees that the eigenfunctions given in Equation – satisfy discrete orthogonality relations, for example, $$\label{eq:orthogorel}
  \frac{2}{N_x+1} \sum_{i=1}^{N_x}{\cal X}_i^l{\cal X}_i^{l'} = \delta_{ll'}\quad\text{and}\quad\frac{2}{N_x+1} \sum_{l=1}^{N_x} {\cal X}_i^l{\cal X}_{i'}^l = \delta_{ii'},$$ where the symbol $\delta$ denotes the Kronecker delta. These orthogonality relations allow us to expand $\widetilde{\Phi}_{i,j,k}$ and $\rho_{i,j,k}$ as $$\begin{aligned}
  \widetilde{\Phi}_{i,j,k} = \frac{8}{(N_x+1)(N_y+1)(N_z+1)}\sum_{l=1}^{N_x}\sum_{m=1}^{N_y}\sum_{n=1}^{N_z}\widetilde{\Phi}^{lmn}{\cal X}_i^l{\cal Y}_j^m{\cal Z}_k^n,\label{eq:car_Phi_forward}\\
  \rho_{i,j,k} = \frac{8}{(N_x+1)(N_y+1)(N_z+1)}\sum_{l=1}^{N_x}\sum_{m=1}^{N_y}\sum_{n=1}^{N_z}\rho^{lmn}{\cal X}_i^l{\cal Y}_j^m{\cal Z}_k^n,\label{eq:car_rho_forward}\end{aligned}$$ where the expansion coefficients $\widetilde{\Phi}^{lmn}$ and $\rho^{lmn}$ are given by $$\begin{aligned}
  \widetilde{\Phi}^{lmn} &= \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}\sum_{k=1}^{N_z}\widetilde{\Phi}_{i,j,k}{\cal X}_i^l{\cal Y}_j^m{\cal Z}_k^n,\label{eq:car_Phi_backward}\\
  \rho^{lmn} &= \sum_{i=1}^{N_x}\sum_{j=1}^{N_y}\sum_{k=1}^{N_z}\rho_{i,j,k}{\cal X}_i^l{\cal Y}_j^m{\cal Z}_k^n.\label{eq:car_rho_backward}\end{aligned}$$ Plugging Equations – in Equation , we obtain a simple algebraic relation $$\label{eq:car_Poisson_transform}
  \widetilde{\Phi}^{lmn} = \frac{4\pi G \rho^{lmn}}{\lambda_x^l + \lambda_y^m + \lambda_z^n}.$$ Therefore, the Poisson equation in Cartesian coordinates can be solved by the following three steps:

1.  Perform a forward transform $\rho_{i,j,k}\to\rho^{lmn}$ using Equation : ${\cal O}(N_xN_yN_z\log_2[N_xN_yN_z])$.

2.  Convert $\rho^{lmn} \to \widetilde{\Phi}^{lmn}$ using the kernel in Equation : ${\cal O}(N_xN_yN_z)$.

3.  Perform a backward transform $\widetilde{\Phi}^{lmn}\to\widetilde{\Phi}_{i,j,k}$ using Equa