Numerical calculation of dipolar quantum droplet stationary states

Several works Saito (2016); Ferrier-Barbut et al. (2016); Wächtler and Santos (2016b); Schmitt et al. (2016); Wächtler and Santos (2016a); Bisset et al. (2016); Baillie et al. (2016); Chomaz et al. (2016); Boudjemâa (2015, 2016, 2017); Ołdziejewski and Jachymski (2016); Macia et al. (2016) have established that the ground states and dynamics of a dipolar condensate in the droplet regime is well-described by the EGPE. In this formalism the time-independent stationary state wavefunction $\Psi$ is a solution of $\mu\Psi=\mathcal{L}_{\text{EGP}}[\Psi]$ where

$\displaystyle\mathcal{L}_{\text{EGP}}[\Psi]\!$

$\displaystyle\equiv\!\biggl{[}-\frac{\hbar^{2}\nabla^{2}}{2M}\!+\!V_{\text{tr}}\!+\!g_{s}|\Psi|^{2}\!+\!g_{dd}\Phi(\mathbf{x})\!+\!\gamma_{\mathrm{QF}}|\Psi|^{3}\biggr{]}\!\Psi.$

(1)

Here $V_{\text{tr}}$ describes any trapping potential, $\mu$ is the chemical potential and $g_{s}=4\pi a_{s}\hbar^{2}/M$ is the $s$-wave coupling constant, with $a_{s}$ being the $s$-wave scattering length. The potential

$\displaystyle\Phi(\mathbf{x})=\int d\mathbf{x}^{\prime}\,I(\mathbf{x}\!-\!\mathbf{x}^{\prime})|\Psi(\mathbf{x}^{\prime})|^{2},$

(2)

describes the effects of the long-ranged DDIs where the kernel is

$\displaystyle I(\mathbf{r})$

$\displaystyle=\frac{3}{4\pi}\frac{1-3\cos^{2}\theta}{r^{3}}.$

(3)

This is written for the case of dipoles polarized along the $z$ axis by an external field, where $\theta$ is the angle between $\mathbf{r}$ and the $z$ axis. The dipole coupling constant is $g_{dd}=4\pi a_{dd}\hbar^{2}/M$, where $a_{dd}=M\mu_{0}\mu_{m}^{2}/12\pi\hbar^{2}$ is the dipole length determined by the magnetic moment $\mu_{m}$ of the particles. The leading-order quantum fluctuation correction to the chemical potential for a uniform system of density $n$ is $\Delta\mu=\gamma_{\mathrm{QF}}n^{3/2}$, where $\gamma_{\mathrm{QF}}\!=\!\frac{32}{3}g_{s}\sqrt{\frac{a_{s}^{3}}{\pi}}(1+\tfrac{3}{2}\epsilon_{dd}^{2})$ and $\epsilon_{dd}\equiv a_{dd}/a_{s}$ Lima and Pelster (2011); Ferrier-Barbut et al. (2016); Bisset et al. (2016). The effects of quantum fluctuations are included in Eq. (1) using the local density approximation $n\rightarrow|\Psi(\mathbf{x})|^{2}$. Stationary EGPE states are also local minima of the energy functional

$\displaystyle E[\Psi]=E_{\text{kin}}+E_{\text{tr}}+E_{\text{int}}+E_{\text{QF}},$

(4)

where

	$\displaystyle E_{\text{kin}}$	$\displaystyle=-\frac{\hbar^{2}}{2M}\int d\mathbf{x}\Psi^{*}\nabla^{2}\Psi,$		(5)
	$\displaystyle E_{\text{tr}}$	$\displaystyle=\int d\mathbf{x}\Psi^{*}V_{\text{tr}}\Psi,$		(6)
	$\displaystyle E_{\text{int}}$	$\displaystyle=\frac{1}{2}\int d\mathbf{x}\Psi^{*}\left(g_{s}|\Psi|^{2}+g_{dd}\Phi\right)\Psi,$		(7)
	$\displaystyle E_{\text{QF}}$	$\displaystyle=\frac{2}{5}\gamma_{\text{QF}}\int d\mathbf{x}|\Psi|^{5},$		(8)

represent the kinetic, trap potential, interaction and quantum fluctuation energies, respectively.

We now write the problem in a form utilizing the cylindrical symmetry and introducing natural units. While we do not present any results here that include a trapping potential (we focus on self-bound states in the absence of confinement), for generality we include a cylindrically symmetric trapping potential in our formulation. We take this to be of the form $V_{\mathrm{tr}}=\frac{1}{2}M(\omega_{\rho}^{2}\rho^{2}+\omega_{z}^{2}z^{2})$, where $\rho=\sqrt{x^{2}+y^{2}}$, and $\{\omega_{\rho},\omega_{z}\}$ are the trap frequencies. For this case the system is cylindrically symmetric (since we have also chosen the dipoles to be along $z$) and we can choose to look for stationary solutions of the form

$\displaystyle\Psi(\mathbf{x})=\psi(\rho,z)e^{is\phi},$

(9)

where $\psi$ is real, $\phi=\arctan(y/x)$ and $s$ is an integer specifying the $z$-component angular momentum of the state. By separating variables, and introducing units of length $x_{0}=a_{dd}$ and energy $E_{0}=\hbar^{2}/Ma_{dd}^{2}$, we arrive at the effective cylindrical GPE

	$\displaystyle\mu\psi$	$\displaystyle=\mathcal{L}[\psi],$		(10)
	$\displaystyle\mathcal{L}[\psi]$	$\displaystyle\equiv h_{sp}\psi+4\pi\left({\epsilon_{dd}^{-1}}\psi^{2}+\Phi\right)\psi+\gamma_{\mathrm{QF}}|\psi|^{3}\psi.$		(11)

Here

$\displaystyle h_{sp}=-\frac{1}{2}\left(D_{s}+\frac{\partial^{2}}{\partial z^{2}}\right)+\frac{1}{2}({\bar{\omega}}_{\rho}^{2}\rho^{2}+{\bar{\omega}}_{z}^{2}z^{2}),$

(12)

is the single particle Hamiltonian, which features the Bessel differential operator

$\displaystyle D_{s}\equiv\frac{\partial^{2}}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial}{\partial\rho}-\frac{s^{2}}{\rho^{2}},$

(13)

and $\bar{\omega}_{\nu}=\hbar\omega_{\nu}/E_{0}$, $\nu=\{\rho,z\}$. We also note that for this choice of units the nonlinear coupling constants are $g_{s}\to 4\pi\epsilon_{dd}^{-1}$, $g_{dd}\to 4\pi$, and

$\displaystyle\gamma_{\mathrm{QF}}=\frac{4}{3\pi^{2}}\left(\frac{4\pi}{\epsilon_{dd}}\right)^{5/2}\left(1+\frac{3}{2}\epsilon_{dd}^{2}\right),$

(14)

thus all specified in terms of $\epsilon_{dd}$.

The convolution used to evaluate $\Phi$ is performed in three-dimensions, but the resulting $\Phi$ is a cylindrically symmetric function:

	$\displaystyle\Phi(\rho,z)$	$\displaystyle\equiv\int d\mathbf{x}^{\prime}\,I(\mathbf{x}-\mathbf{x}^{\prime})|\psi(\rho^{\prime},z^{\prime})|^{2},$		(15)
		$\displaystyle=\int\frac{dk_{\rho}\,dk_{z}}{(2\pi)^{2}}e^{ik_{z}z}k_{\rho}J_{0}(k_{\rho}\rho)\tilde{I}(k_{\rho},k_{z})\tilde{n}(k_{\rho},k_{z}),$		(16)

where the Fourier transformed density and DDI kernel are

	$\displaystyle\tilde{n}(k_{\rho},k_{z})$	$\displaystyle=2\pi\int d\rho\,dz\,\rho J_{0}(k_{\rho}\rho)e^{-ik_{z}z}|\psi(\rho,z)|^{2},$		(17)
	$\displaystyle\tilde{I}(k_{\rho},k_{z})$	$\displaystyle=\frac{3k_{z}^{2}}{k_{\rho}^{2}+k_{z}^{2}}-1,$		(18)

respectively. We also note that here we choose the condensate to be normalized to the number of atoms $N$, i.e.

$\displaystyle 2\pi\int_{0}^{\infty}d\rho\int_{-\infty}^{\infty}dz\,\rho|\psi|^{2}=N.$

(19)

In the radial direction we consider $N_{\rho}$ points, that non-uniformly span the interval $(0,R)$, given by

$\displaystyle\rho_{qi}=\frac{\alpha_{qi}}{K_{q}},\qquad i=1,\ldots,N_{\rho},$

(20)

which we refer to as the $q$-order Bessel grid. Here $K_{q}={\alpha_{q\,N_{\rho}+1}}/R$ , and $\{\alpha_{qi}\}$ are the ordered non-zero roots of the Bessel function $J_{q}(x)$ of integer order $q$. In this work we will need several such radial grids, each with the same number of points and range but of different orders. For this reason we need to adopt a notation that explicitly indicates the order. We also introduce the reciprocal space grid of the same order

$\displaystyle k_{qi}=\frac{\alpha_{qi}}{R},\qquad i=1,\ldots,N_{\rho},$

(21)

that spans the interval $(0,K_{q})$.

The radial grid points can be associated with a quadrature-like integration formula Lemoine (1994); Ben Ghanem and Frappier (1999)

$\displaystyle I[g(\rho)]=\int_{0}^{\infty}d\rho\,\rho g(\rho)\approx\sum_{i=1}^{N_{\rho}}w_{qi}\,g_{qi},$

(22)

where $w_{qi}$ are the real-space quadrature weights

$\displaystyle w_{qi}$

$\displaystyle=\frac{2}{K_{q}^{2}|J_{q+1}(\alpha_{qi})|^{2}},$

(23)

and we have used the notation $g_{qi}=g(\rho_{qi})$ to denote the function $g(\rho)$ sampled on the $q$-order grid. This integration requires that the functions of interest have limited spatial range, i.e. $g(\rho>R)=0$.

Similarly, for the reciprocal $k_{\rho}$-space we have a quadrature-like integration formula for a function $\tilde{g}(k_{\rho})$

$\displaystyle I[\tilde{g}(k_{\rho})]=\int_{0}^{\infty}dk_{\rho}\,k_{\rho}\tilde{g}(k_{\rho})\approx\sum_{i=1}^{N_{\rho}}\tilde{w}_{qi}\,\tilde{g}_{qi},$

(24)

where $\tilde{w}_{qi}$ are the $k_{\rho}$-space quadrature weights

$\displaystyle\tilde{w}_{qi}$

$\displaystyle=\frac{2}{R^{2}|J_{q+1}(\alpha_{qi})|^{2}},$

(25)

and $\tilde{g}_{qi}=\tilde{g}(k_{qi})$. Result (24) is only valid on our grid if the function is bandwidth limited, i.e. $\tilde{g}(k_{\rho}>K_{q})=0$.

The Bessel grid is useful because it allows an accurate two-dimensional (2D) Fourier transformation of functions of the form

$\displaystyle F(\bm{\rho})=f(\rho)e^{im\phi},$

(26)

where we have used $\bm{\rho}$ to denote the planar position vector with polar coordinates ($\rho,\phi)$. The 2D Fourier transform of $F$ is

$\displaystyle\tilde{F}(\mathbf{k}_{\rho})$

$\displaystyle=\int d\bm{\rho}\,e^{-i\mathbf{k}_{\rho}\cdot\bm{\rho}}F(\bm{\rho})=2\pi i^{-m}e^{im\phi_{k}}\tilde{f}(k_{\rho}),$

(27)

where

$\displaystyle\tilde{f}(k_{\rho})=$

$\displaystyle\int_{0}^{\infty}\!d\rho\,\rho J_{m}(k_{\rho}\rho)\,f(\rho),$

(28)

is the $m$-th order Hankel transform, arising because the angular integral of $e^{im\phi-i\mathbf{k}_{\rho}\cdot\bm{\rho}}$ yields the $J_{m}$ Bessel function. Here we have introduced $\mathbf{k}_{\rho}$ to represent the 2D $k$-space vector, with polar coordinates $(k_{\rho},\phi_{k})$.

For the case where the angular momentum of the function [i.e., $m$ from Eq. (26)] and the order of the grid used to sample the function are the same (i.e. $q=m$) an accurate discrete Hankel transform can be implemented. For this case the transform yields $\tilde{f}$ sampled on the $k_{qi}$ grid from $f$ sampled on the $\rho_{qi}$ grid (see Guizar-Sicairos and Gutiérrez-Vega (2004)). The explicit form of this discrete transform is obtained by evaluating (28) using the quadrature formula (22)

$\displaystyle\tilde{f}_{qi}=\sum_{j}\mathcal{H}_{ij}{f}_{qj},$

(29)

where

$\displaystyle\mathcal{H}_{ij}=w_{qj}J_{q}(k_{qi}\rho_{qj})=\frac{2J_{q}\left(\frac{\alpha_{qi}\alpha_{qj}}{\alpha_{qN_{\rho}+1}}\right)}{K_{q}^{2}|J_{q+1}(\alpha_{qj})|^{2}},$

(30)

is the $q$-order Hankel transformation matrix and ${f}_{qj}={f}(\rho_{qj})$ etc.

Similarly, the inverse 2D transform

$\displaystyle F(\bm{\rho})=\int d\mathbf{k}_{\rho}\frac{e^{i\mathbf{k}_{\rho}\cdot\bm{\rho}}}{(2\pi)^{2}}\tilde{F}(\mathbf{k}_{\rho})=\frac{i^{-q}e^{iq\phi}}{2\pi}{f}(\rho),$

(31)

is accomplished by the inverse Hankel transform

$\displaystyle f(\rho)=\int_{0}^{\infty}dk_{\rho}\,k_{\rho}J_{q}(k_{\rho}\rho)\tilde{f}(k_{\rho}),$

(32)

with discrete form ${f}_{qi}=\sum_{j}\mathcal{H}_{ij}^{-1}\tilde{f}_{qj}$, with $\mathcal{H}_{ij}^{-1}=\frac{K_{q}^{2}}{R^{2}}\mathcal{H}_{ij}$. The Discrete Hankel transform matrices are not exact inverses of each other, but for typical grid sizes ($N_{\rho}\sim 10^{2}$) we have that $\sum_{j}\mathcal{H}_{ij}\mathcal{H}_{jk}^{-1}\approx\delta_{ik}+O(10^{-9})$, which is adequate for our purposes.

Hankel transforms are particularly useful for accurately evaluating the kinetic energy operator in radially symmetric cases (e.g. see Bisseling and Kosloff (1985)). To see this we note that by separating variable the 2D Laplacian $\nabla^{2}_{\bm{\rho}}$ acting on the function $f(\rho)e^{is\phi}$ is equivalent to the Bessel differential operator $D_{s}$ (13) acting on $f(\rho)$ [cf. Eq. (12)]. Using that $J_{s}$ are eigenfunctions of $D_{s}$, i.e. $D_{s}J_{s}(k_{\rho}\rho)=-k_{\rho}^{2}J_{s}(k_{\rho}\rho)$, we can utilize the Hankel transform to act on a radial function with $D_{s}$:

$\displaystyle[D_{s}f]_{i}\approx-\sum_{jk}\mathcal{H}^{-1}_{ij}{k_{sj}^{2}}\mathcal{H}_{jk}f_{k}.$

(33)

Allowing us to implement a spectrally accurate radial kinetic energy matrix [i.e. discretized version of $T_{\rho}=-\frac{1}{2}D_{s}$, c.f. Eq. (12)]:

$\displaystyle T_{\rho,ij}=\frac{1}{2}\sum_{jk}\mathcal{H}^{-1}_{il}k_{sl}^{2}\mathcal{H}_{lj}.$

(34)

Our system has reflection symmetry along $z$ so functions of interest will be of definite parity and here we will concern ourselves with the case of even parity. Utilizing this symmetry we use a half-grid on the interval $(0,Z)$ spanned by $N_{z}$ equally spaced points

$\displaystyle z_{j}=(j-\tfrac{1}{2})\Delta z,\qquad j=1,2,\ldots N_{z},$

(35)

where $\Delta z=Z/N_{z}$. The corresponding reciprocal grid

$\displaystyle k_{j}=(j-\tfrac{1}{2})\Delta k,\qquad j=1,2,\ldots N_{z},$

(36)

spans the interval $(0,K_{z})$ where $\Delta k=\pi/Z$ and $K_{z}=\pi/\Delta z$. The appropriate quadrature for this grid is the rectangular rule

$\displaystyle I[g(z)]=\int_{-\infty}^{\infty}dz\,g(z)\approx\sum_{i=1}^{N_{z}}g_{i}\,2\Delta z,$

(37)

where $g_{i}=g(z_{i})$.

The 1D Fourier transform for an even parity function is given by the cosine transform

	$\displaystyle\tilde{f}(k_{z})$	$\displaystyle=\int_{-\infty}^{\infty}dz\,e^{-ik_{z}z}{f}(z),$		(38)
		$\displaystyle=2\int_{0}^{\infty}dz\cos(k_{z}z){f}(z).$		(39)

We can discretize this transform onto our chosen grid as

$\displaystyle\tilde{f}_{i}$

$\displaystyle=\Lambda_{ij}f_{j},$

(40)

where

$\displaystyle{\Lambda}_{ij}=2\cos(k_{i}z_{j})\Delta z,$

(41)

is the transformation matrix, $f_{j}=f(z_{j})$, and $\tilde{f}_{i}=\tilde{f}(k_{i})$. This can be identified as the type-IV discrete cosine transformation (e.g. see Strang (1999)).

The inverse Fourier transformation ${f}(z)=\frac{1}{2\pi}\int_{-\infty}^{\infty}dk_{z}e^{ik_{z}z}\tilde{f}(k_{z})$ similarly can be mapped to the inverse discrete cosine transform $f_{j}=\Lambda^{-1}_{ij}\tilde{f}_{j}$, where ${\Lambda}_{ij}^{-1}=\frac{K_{z}}{2\pi Z}{\Lambda}_{ij}$, with $\sum_{j}{\Lambda}_{ij}^{-1}{\Lambda}_{jk}=\delta_{ik}$.

Utilizing that the derivative operator is diagonal in $k_{z}$-space, we can implement a discrete second derivative operator as

$\displaystyle\left[\frac{d^{2}f}{dz^{2}}\right]_{i}$

$\displaystyle\approx-\sum_{jk}\Lambda_{ij}^{-1}k_{j}^{2}{\Lambda}_{jk}f_{k}$

(42)

This allows us to define a spectrally accurate axial kinetic energy matrix [i.e. discretized version of $T_{z}=-\frac{1}{2}\frac{d^{2}}{dz^{2}}$, c.f. Eq. (12)]:

$\displaystyle T_{z,ij}=\frac{1}{2}\sum_{jk}\Lambda_{ij}^{-1}k_{l}^{2}{\Lambda}_{lj}.$

(43)

The single particle operator (12) acting on the discretized field $\psi_{ij}=\psi(\bm{\rho}_{ij})$ can be evaluated as

$\displaystyle h_{ij}[\psi]=\sum_{k}T_{\rho,ik}\psi_{kj}+\sum_{k}T_{z,jk}\psi_{ik}+V_{\mathrm{tr},ij}\psi_{ij},$

(44)

where $V_{\mathrm{tr},ij}=\frac{1}{2}(\bar{\omega}_{\rho}^{2}\rho_{si}^{2}+\bar{\omega}_{z}^{2}z_{j}^{2})$, and we have used the Bessel (34) and cosine (43) derivative operators.

The contact interaction term in Eq. (11) is a nonlinear term that is local in position space and is evaluated as

$\displaystyle C_{ij}\left[\psi\right]\equiv 4\pi\epsilon_{dd}^{-1}\psi_{ij}^{3},$

(45)

and similarly for the quantum fluctuation term

$\displaystyle Q_{ij}\left[\psi\right]\equiv\gamma_{\mathrm{QF}}\left|\psi_{ij}\right|^{3}\psi_{ij}.$

(46)

Note: we have assumed that $\psi_{ij}$ is real, but modulus sign is needed on the quantum fluctuation term to properly deal with any case where $\psi_{ij}$ is negative.

The DDI potential $\Phi$ can be evaluated using the convolution theorem in cylindrical coordinates, as given in Eq. (16). To evaluate this expression we first need to Fourier transform the density (17), which we obtain by applying the cosine transform along $z$ and the quadrature rule to evaluate the radial transform

$\displaystyle\tilde{n}_{ij}^{\prime}=2\pi\sum_{kl}\Lambda_{jl}w_{sk}J_{0}(k_{0i}\rho_{sk})\psi_{kl}^{2}.$

(47)

Here we use the prime to indicate that the Fourier transformed density is evaluated on the cylindrical $k$-space mesh of grid points $\mathbf{k}_{ij}^{\prime}=(k_{0i},k_{j})$, noting that the radial $k$-space points are defined for the 0-order Bessel grid¹¹1The radial transform in Eq. (47) is identical to (29) if $s=0$.. We then evaluate the potential $\Phi$ from Eq. (16) as

$\displaystyle\Phi_{ij}$

$\displaystyle=\frac{1}{2\pi}\sum_{kl}\Lambda^{-1}_{jl}\tilde{w}_{0k}J_{0}(k_{0k}\rho_{si})\tilde{I}(\mathbf{k}_{kl}^{\prime})\tilde{n}^{\prime}_{kl},$

(48)

and thus the full DDI term

$\displaystyle{D}_{ij}\left[\psi\right]\equiv 4\pi\Phi_{ij}\psi_{ij}.$

(49)

The function $\tilde{I}(\mathbf{k})$ is of critical importance for accurate numerical calculations of the DDIs. For clarity we hereon refer to the analytic form of this function introduced in Eq. (18) as the bare $k$-space potential $\tilde{I}_{\mathrm{bare}}$. This result is singular since the $k\to 0$ limit does not exist, reflecting the long-ranged anisotropic character of the interaction. As such the quadrature in (48) will converge slowly as the density of $k$-space points increases, or equivalently as the spatial ranges $R,Z$ are made larger (also see discussion in Refs. Ronen et al. (2006); Jiang et al. (2014)). This slow convergence can be understood as the $\Phi$ having contributions from periodic copies (arising from using Fourier transforms on a finite interval) of the density distribution. These periodic copies, separated by twice the grid range, interact with each other causing a shift in the interaction energy. The unphysical influence of these copies decays with the cube of the grid range, making it impractical to calculate accurate results using $\tilde{I}_{\mathrm{bare}}$.

One approach to deal with this problem is to use spherical coordinates in which the Jacobian introduced by the coordinate transform removes the singularity of $\tilde{I}_{\mathrm{bare}}$, but this requires the use of a non-uniform Fourier transform Jiang et al. (2014); Bao et al. (2016) in the algorithm. An alternative approach, first proposed for dipolar BECs in Ref. Ronen et al. (2006) (also see Exl et al. (2016); Mennemann et al. (2019)), is to introduce a truncation of the DDI, i.e. restrict the range of the DDI to the physical extent of the grids used, so that interactions between periodic copies are formally zero (also see Antoine et al. (2018), and related treatments for the Coulomb case Jarvis et al. (1997); Rozzi et al. (2006)). We choose to follow this approach here since it allows us to work with the cylindrical grid and associated transforms that we have introduced. The issue now becomes how to obtain the appropriate truncated DDI potential in $k$-space.

First let us consider a spherical truncation derived and applied to dipolar BECs in Ref. Ronen et al. (2006) (also see Wilson et al. (2009)). Here the truncated interaction in position space is

$\displaystyle{I}_{\mathrm{sph}}({\mathbf{r}})=\begin{cases}\frac{3}{4\pi r^{3}}(1-3\cos^{2}\theta),&r\leq r_{c}\\ 0,&\mathrm{otherwise},\end{cases}$

(50)

i.e. a cutoff radius $r_{c}$ such that the DDIs do not occur between any two points separated by a distance of more than $2r_{c}$. In practice choosing $r_{c}=\min\{R,Z\}$ ensures that no interactions can occur between periodic copies of the density. Fortunately the Fourier transform of ${I}_{\mathrm{sph}}$ can be calculated analytically

$\displaystyle\tilde{I}_{\mathrm{sph}}(\mathbf{k})=$

$\displaystyle\left[1+3\frac{\cos(kr_{c})}{(kr_{c})^{2}}-3\frac{\sin(kr_{c})}{(kr_{c})^{3}}\right]\tilde{I}_{\mathrm{bare}}(\mathbf{k}),$

(51)

and is seen to regularize the behaviour of $\tilde{I}_{\mathrm{bare}}(\mathbf{k})$ near $k=0$. The spherically truncated potential is most useful for situations where the condensate is nearly spherical so it is natural to choose $R\approx Z$.

For highly anisotropic cases the natural grid choice is $R\ll Z$ or $R\gg Z$ for cigar or pancake shaped condensates, respectively. Here the spherical truncation is impractical as the range of the narrow dimension would have to be extended to match the range of the other dimension, introducing many redundant grid points. In these situations it is natural to introduce a cylindrical truncation

$\displaystyle I_{\mathrm{cyl}}({\mathbf{r}})=\begin{cases}\frac{3}{4\pi r^{3}}(1-3\cos^{2}\theta),&\mathbf{r}\in\{\rho<R,|z|<Z\}\\ 0,&\mathrm{otherwise}.\end{cases}$

(52)

The Fourier transform of this truncated kernel, $\tilde{I}_{\mathrm{cyl}}$, does not have an analytic result and needs to be calculated numerically. We refer the interested reader to Ref. Lu et al. (2010) for details about the numerical calculation.

In the testing we perform we evaluate the DDI term (49) using $\Phi_{ij}$ evaluated according to (48) with $\tilde{I}$ set to one of $\{\tilde{I}_{\mathrm{bare}},\,\tilde{I}_{\mathrm{sph}},\,\tilde{I}_{\mathrm{cyl}}\}$. If not specified otherwise, we use $\tilde{I}_{\mathrm{cyl}}$