Content deleted Content added
dimensions |
Swasthyapro (talk | contribs) Link suggestions feature: 3 links added. |
||
| (8 intermediate revisions by one other user not shown) | |||
Line 64:
</math>
where <math>k</math> is the [[wavenumber]], <math>\mathbf{I}</math> is the [[identity matrix]], and <math>\mathbf{r}</math> is the vector from the source dipole to the observation point. Evaluating the derivatives leads to the explicit form:
:<math>
Line 73:
</math>
where <math>\hat{\mathbf{r}} = \mathbf{r} / |\mathbf{r}|</math> is the [[unit vector]] pointing from the source to the observation point.
This Green’s tensor describes the electric field generated by a dipole in a homogeneous medium. It is used to compute the off-diagonal blocks of the interaction matrix in DDA, that is, the interaction between distinct dipoles <math>j \ne k</math>. The singular self-term <math>\mathbf{G}(\mathbf{r} = 0)</math> is excluded and replaced by a prescribed local term involving the inverse polarizability tensor <math>\boldsymbol{\alpha}_j^{-1}</math>.
Line 246:
Each dipole has three vector components (<math>x</math>, <math>y</math>, <math>z</math>), so we can rearrange the unknown vector <math>\mathbf{P}</math> by grouping all x-components together, then y-components, then z-components:
<math>▼
:<math>
\mathbf{P} =
\begin{bmatrix}
Line 252 ⟶ 253:
\mathbf{P}_y \\
\mathbf{P}_z
\end{bmatrix}
\in \mathbb{C}^{3N} \quad \text{where} \quad \mathbf{P}_x =
\begin{bmatrix}
Line 277 ⟶ 280:
Similarly, the incident field can be grouped as:
:<math>
<math>▼
▲\mathbf{E}^{\mathrm{inc}} =
\begin{bmatrix}
\mathbf{E}_x^{\mathrm{inc}} \\
Line 285 ⟶ 287:
\mathbf{E}_z^{\mathrm{inc}}
\end{bmatrix}
\in \mathbb{C}^{3N}
</math>
Because the system is linear, we can equivalently rewrite it in block matrix form, that describe how the <math>\beta</math>-component of polarization affects the <math>\alpha</math>-component of the resulting field:
Line 307 ⟶ 311:
\end{bmatrix}
</math>
Each matrix-vector multiplication <math>\mathbf{A}_{\alpha\beta} \mathbf{P}_\beta</math> can be computed as a convolution when the dipoles are arranged on a regular grid, allowing the use of Fast Fourier Transforms (FFTs) to accelerate the solution. ▼
The expanded form of the equations is:
Line 314 ⟶ 316:
<math>
\begin{aligned}
-\mathbf{
(\boldsymbol{\alpha}^{-1} \mathbf{P})_x &= \mathbf{E}_x^{\mathrm{inc}} \\ -\mathbf{
(\boldsymbol{\alpha}^{-1} \mathbf{P})_y &= \mathbf{E}_y^{\mathrm{inc}} \\ -\mathbf{
(\boldsymbol{\alpha}^{-1} \mathbf{P})_z &= \mathbf{E}_z^{\mathrm{inc}} \end{aligned}
</math>
Each block <math> \mathbf{G}_{ij} \in \mathbb{C}^{N \times N} </math> and the total system size is <math> 3N \times 3N </math>. The interaction matrix <math> \mathbf{G} \in \mathbb{C}^{3N \times 3N} </math> is composed of 9 blocks:
<math> \mathbf{G}_{xx}, \mathbf{G}_{xy}, \mathbf{G}_{xz}, \mathbf{G}_{yx}, \mathbf{G}_{yy}, \mathbf{G}_{yz}, \mathbf{G}_{zx}, \mathbf{G}_{zy}, \mathbf{G}_{zz} </math> (only 6 of them need to be evaluated due to symmetry).
▲Each matrix-vector multiplication <math>\mathbf{
Let <math>\boldsymbol{\beta}_j = \boldsymbol{\alpha}_j^{-1}</math> denote the inverse polarizability tensor for dipole <math>j</math>. Each <math>\boldsymbol{\beta}_j</math> is a complex-valued <math>3 \times 3</math> matrix. This gives:
▲<math>
(\boldsymbol{\beta} \mathbf{P})_x =
\mathrm{diag}(\beta_{1,xx}, \dots, \beta_{N,xx}) \, \mathbf{P}_x +
\mathrm{diag}(\beta_{1,xy}, \dots, \beta_{N,xy}) \, \mathbf{P}_y +
\mathrm{diag}(\beta_{1,xz}, \dots, \beta_{N,xz}) \, \mathbf{P}_z
</math>
▲<math>
(\boldsymbol{\beta} \mathbf{P})_y =
\mathrm{diag}(\beta_{1,yx}, \dots, \beta_{N,yx}) \, \mathbf{P}_x +
\mathrm{diag}(\beta_{1,yy}, \dots, \beta_{N,yy}) \, \mathbf{P}_y +
\mathrm{diag}(\beta_{1,yz}, \dots, \beta_{N,yz}) \, \mathbf{P}_z
</math>
<math>
(\boldsymbol{\beta} \mathbf{P})_z =
\mathrm{diag}(\beta_{1,zx}, \dots, \beta_{N,zx}) \, \mathbf{P}_x +
\mathrm{diag}(\beta_{1,zy}, \dots, \beta_{N,zy}) \, \mathbf{P}_y +
\mathrm{diag}(\beta_{1,zz}, \dots, \beta_{N,zz}) \, \mathbf{P}_z
</math>
In the special case of an isotropic and homogeneous particle, the polarizabilities <math>\boldsymbol{\alpha}_j</math> are identical for all dipoles and proportional to the identity matrix: <math>\boldsymbol{\alpha}_j = \alpha \, \mathbf{I}</math>. Then, the inverse becomes <math>\boldsymbol{\beta}_j = \alpha^{-1} \, \mathbf{I}</math>, all off-diagonal elements vanish, and the expressions reduce to a simple element-wise division:
<math>
(\boldsymbol{\beta} \mathbf{P}) = \frac{1}{\alpha} \, \mathbf{P}
</math>
Note on practical implementation. In Fortran and MATLAB, arrays such as <math>\mathbf{P}(n_x,n_y,n_z,3)</math> or <math>\mathbf{E}^{\mathrm{inc}}(n_x,n_y,n_z,3)</math> are stored in column-major order, where the first index varies fastest in memory (anti-lexicographic). This means that all x-components <math>\mathbf{P}_x = \mathbf{P}(:,:,:,1)</math> are contiguous in memory, followed by all y-components <math>\mathbf{P}_y = \mathbf{P}(:,:,:,2)</math>, and then all z-components <math>\mathbf{P}_z = \mathbf{P}(:,:,:,3)</math>. In contrast, Python (NumPy) uses row-major order by default (lexicographic, last index varies fastest). To achieve the same contiguous layout of <math>\mathbf{P}_x</math>, <math>\mathbf{P}_y</math>, <math>\mathbf{P}_z</math> in memory, the array should be defined in Python as <math>\mathbf{P}(3, n_x, n_y, n_z)</math>, with the vector component index (x, y, z) first. This ensures that <math>\mathbf{P}_x = \mathbf{P}[0, :, :, :]</math> is stored contiguously in memory, followed by <math>\mathbf{P}_y</math> and <math>\mathbf{P}_z</math>.
== Green's function in Fourier space==
In the Fourier ___domain, the free-space dyadic Green's function for the vector Helmholtz equation can be written as
<math>
\widehat{\mathbf{G}}(\mathbf{k}) =
S(\mathbf{k})\left( \mathbf{I} - \frac{\mathbf{k}\,\mathbf{k}^{\mathrm{T}}}{|\mathbf{k}|^{2}} \right),
\qquad
S(\mathbf{k}) = \frac{1}{|\mathbf{k}|^{2} - k_{0}^{2} + i0}\,.
</math>
Here <math>\widehat{\mathbf{G}}(\mathbf{k})</math> is the Fourier transform of the dyadic Green's function <math>\mathbf{G}(\mathbf{r})</math>, <math>\mathbf{k} = (k_{x},k_{y},k_{z})</math> is the spectral wavevector, and <math>\mathbf{I}</math> is the identity tensor. The term <math>\mathbf{k}\,\mathbf{k}^{\mathrm{T}}/|\mathbf{k}|^{2}</math> is the longitudinal projector onto the direction of <math>\mathbf{k}</math>, so <math>\mathbf{I} - \mathbf{k}\,\mathbf{k}^{\mathrm{T}}/|\mathbf{k}|^{2}</math> is the transverse projector. The scalar factor <math>S(\mathbf{k})</math> contains the [[Helmholtz equation|Helmholtz denominator]] <math>|\mathbf{k}|^{2} - k_{0}^{2}</math>, where <math>k_{0} = \omega / c</math> is the wavenumber in the embedding medium. The infinitesimal imaginary term <math>+i0</math> implements the [[Sommerfeld radiation condition]], ensuring that the inverse Fourier transform produces an outgoing-wave solution.
Line 328 ⟶ 381:
In practice, the dominant computational cost in DDA arises from the repeated evaluation of matrix-vector products during the iteration process. When the vector <math>\mathbf{P}</math> is stored in component-block form (as <math>\mathbf{P}_x</math>, <math>\mathbf{P}_y</math>, <math>\mathbf{P}_z</math>), the action of <math>\mathbf{A}</math> reduces to evaluating nine sub-products of the form <math>\mathbf{A}_{\alpha\beta} \mathbf{P}_\beta</math>, where <math>\alpha,\beta \in \{x,y,z\}</math>. These operations can be computed efficiently using convolution and FFT-based techniques when the dipole geometry is grid-based.
==Fast Fourier Transform for fast convolution calculations==
The use of the fast Fourier transform (FFT) to accelerate convolution operations in the discrete dipole approximation (DDA) was introduced by Goodman, Draine, and Flatau in 1991{{r|Goodman1991}}. Their method employed the three-dimensional FFT algorithm (GPFA) developed by Clive Temperton{{r|Temperton1983}}and required extending the interaction matrix from its original size <math>(n_x,n_y,n_z)</math> to <math>(2n_x,2n_y,2n_z)</math>. This extension was achieved by reversing and mirroring the Green’s function tensor blocks so that all positive and negative spatial offsets (lags) were represented in a single array, with an inserted zero plane between the positive and negative sides along each axis. This arrangement ensured that the discrete convolution of the Green’s function with the polarization vector could be performed as a cyclic convolution using FFTs, avoiding aliasing from wraparound effects. The sign-flipping of the Green’s function in the frequency ___domain and the block extension procedure became standard steps in efficient DDA implementations. Several alternative formulations have since been proposed.
A similar variant to that of Goodman, Draine, and Flatau was adopted in the 2021 MATLAB implementation by Shabaninezhad and Ramakrishna{{r|matlab2021}}. In this approach, the computational ___domain for the polarization vector is zero-padded to <math>(2n_x-1)\times (2n_y-1)\times (2n_z-1)</math> instead of the <math>2n_x\times 2n_y\times 2n_z</math> sizes used in ''DDSCAT''. The stored interaction matrix <math>G</math> differs from ''DDSCAT'' in that there is no zero plane inserted between the positive and negative offsets along each axis. The FFT is performed as a sequence of one-dimensional transforms along the <math>x</math>, <math>y</math>, and <math>z</math> axes, which is mathematically equivalent to performing a full 3-D FFT on the padded ___domain.
Sequence of 1D FFTs was used by MacDonald in his Ph. D. thesis <ref name=mcdonald2009/>.
Barrowes method is a numerical technique for multiplying an <math>n</math>-dimensional block Toeplitz matrix by a vector using the fast Fourier transform (FFT). In three dimensions with grid sizes <math>n_x</math>, <math>n_y</math>, and <math>n_z</math>, the method embeds the block Toeplitz array into a larger block circulant array of size <math>(2n_x-1)\times (2n_y-1)\times (2n_z-1)</math>, which ensures that the corresponding convolution is free of cyclic wraparound. The kernel—containing all positive and negative offsets with the self term set to zero—is reversed along the offset axes, flattened into a one-dimensional array, and transformed by a single long FFT. The input vector is likewise placed in a zero-padded ___domain of the same size, flattened, and transformed. An element-wise product in the frequency ___domain corresponds to the spatial-___domain convolution; an inverse FFT is then reshaped and cropped back to the physical ___domain to obtain the result. The method applies to arbitrary dimension <math>n</math> and block size, and was used originally in the discrete dipole approximation.<ref name="Barrowes2001" />
Line 396 ⟶ 457:
<ref name=chaumet2022discrete>{{cite journal |last=Chaumet |first=Patrick C. |title=The discrete dipole approximation: A review |journal=Mathematics |volume=10 |issue=17 |pages=3049 |year=2022 |doi=10.3390/math10173049 |doi-access=free }}</ref>
<ref name=Barrowes2001>{{cite journal |last1=Barrowes |first1=B. E. |last2=Teixeira |first2=F. L. |last3=Kong |first3=J. A. |title=Fast algorithm for matrix–vector multiply of asymmetric multilevel block-Toeplitz matrices in 3-D scattering |journal=Microw. Opt. Technol. Lett. |volume=31 |issue=1 |pages=28–32 |year=2001 |doi=10.1002/mop.1348}}</ref>
Line 407 ⟶ 466:
<ref name="Moncada-Villa2022">{{Cite journal | last1 = Moncada-Villa | first1 = E. | last2 = Cuevas | first2 = J. C. | year = 2022 | title = Thermal discrete dipole approximation for near-field radiative heat transfer in many-body systems with arbitrary nonreciprocal bodies | journal = Physical Review B | volume = 106 | issue = 23 | pages = 235430 | doi = 10.1103/PhysRevB.106.235430 | arxiv = 2206.14921 | bibcode = 2022PhRvB.106w5430M }}</ref>
<ref name="Temperton1983">C. Temperton. "Self-sorting mixed-radix fast Fourier transforms." Journal of Computational Physics, 52.1 (1983): 1–23.</ref>
}}
| |||