Calculation of the volumetric diffracted field with a three-dimensional convolution: the three-dimensional angular spectrum method

Shan Shan Kou; Colin J. R. Sheppard; Jiao Lin

doi:10.1364/OL.38.005296

The angular spectrum (AS) method [1] is a widely used Fourier-transform-based technique for calculating propagation of monochromatic wave fields. It has found numerous applications in areas such as digital holography [2,3], where the calculation is required over a volume. Since the classical AS method is based on a two-dimensional (2D) theory, one has to propagate to obtain a stack of 2D field distributions in order to reconstruct the volumetric wave field. Here we formulate a theory and the associated numerical recipe to show how the entire volume of a wave field can be obtained using a single pair of Fourier transforms (FTs). In fact, the proposed AS method gives the exact relationship of the fundamental first Rayleigh–Sommerfeld diffraction formula (RS1) in terms of a 3D convolution. We therefore term our approach the three-dimensional angular spectrum (3D-AS) method.

It is well known that the RS1 allows the scalar diffracted field at any point to be calculated rigorously given the distribution of the field in the initial plane [4–6]. This is in contrast to the Huygens–Fresnel principle (HFP), which simplifies the RS1 by assuming the point of interest is sufficiently far away as compared to the wavelength [4]. For more generality, in this Letter we will develop the 3D-AS method based on the fundamental RS1, although a similar treatment can be carried out starting from the HFP.

Starting from the RS1 [4], for the field to be calculated at an arbitrary point $Q$ , in terms of the known field $U_{0}$ at point $P$ in the original plane (Fig. 1), we have

U (x, y, z) = \frac{1}{2 π} \iint_{Σ} U_{0} (ξ, η, 0) K (ξ, η, 0) d ξ d η,

where the kernel

K

is given by

K (ξ, η, ζ) = \frac{\partial}{\partial ζ} [\frac{\exp (i k ρ)}{ρ}] .

Here the distance between

P

and

Q

is

ρ = {[{(x - ξ)}^{2} + {(y - η)}^{2} + {(z - ζ)}^{2}]}^{1 / 2}

, and

k = 2 π / λ

, where

λ

is the wavelength of the monochromatic field. Various approximations [1,7,8] for expanding the square root in the expression of the distance can be introduced to simplify the calculation. We have recently demonstrated that under the paraxial approximation, a 3D diffracted field can be calculated by a single 3D FoFT in the spatial domain [9]. This method has its roots in McCutchen’s approach for the Debye approximation [10,11]. In this Letter, we show that RS1 can be formulated exactly as a 3D convolution without introducing any approximation. Therefore, the method is valid even close to the original plane, and no paraxial approximation is made. Equation (1) can be represented as a 3D integral,

U (x, y, z) = \frac{1}{2 π} ∭_{V} V_{0} (ξ, η, ζ) K (ξ, η, ζ) d ξ d η d ζ,

with

V_{0} (ξ, η, ζ) = U_{0} (ξ, η, 0) δ (ζ)

, where

δ (•)

is the Dirac delta function. According to the convolution theorem, Eq. (3) can be evaluated in an exact manner by FTs:

U (x, y, z) = - \frac{1}{2 π} F_{3}^{- 1} {F_{2} {U_{0} (ξ, η, 0)} \times F_{3} {- K (ξ, η, ζ)}} .

Here,

F_{n} {•}

and

F_{n}^{- 1} {•}

denote an

n

-dimensional FT and inverse FT, respectively. Note that the relationship

F_{3} {V_{0} (ξ, η, ζ)} = F_{2} {U_{0} (ξ, η, 0)}

has been used. What differentiates our approach from the traditional 2D-AS method is that the FT of the incident wave function is now multiplied by a 3D transfer function, rather than a 2D counterpart [1,12]. Using the differentiation theorem of FTs, the 3D transfer function can be written as

\hat{K} (f_{x}, f_{y}, f_{z}) = - i 2 π f_{z} \hat{G} (f_{x}, f_{y}, f_{z}),

where

\hat{G} (•)

is the 3D-FT of the free-space Green’s function. Employing the spherical symmetry and evaluating the above using the Fourier sine transform [13],

\hat{G} (f_{r}) = 4 π \int_{0}^{\infty} \frac{\exp (i k r)}{r} \frac{\sin (2 π f_{r} r)}{2 π f_{r} r} r^{2} d r,

where

r = {[ξ^{2} + η^{2} + ζ^{2}]}^{1 / 2}

. In the above,

\hat{G}

can be viewed as a separation into homogeneous and inhomogeneous parts,

\hat{G} = {\hat{G}}_{A} + {\hat{G}}_{B}

[14]. Given the fact that

f_{r} = \sqrt{f_{x}^{2} + f_{y}^{2} + f_{z}^{2}} > 0

and

δ (f_{r} + 1 / λ) = 0

, the two parts can be defined as

{\hat{G}}_{A} = \frac{i}{2 f_{r}} [δ (f_{r} - 1 / λ) - δ (f_{r} + 1 / λ)] = i δ (f_{r}^{2} - 1 / λ^{2}), {\hat{G}}_{B} = \frac{1}{2 π f_{r}} [\frac{1}{(f_{r} - 1 / λ)} + \frac{1}{(f_{r} + 1 / λ)}] = \frac{1}{π (f_{r}^{2} - 1 / λ^{2})} .

We have clearly demonstrated that the RS1 can be regarded as a 3D convolution analytically. If we substitute

(f_{r}^{2} - 1 / λ^{2}) = (f_{z} + M) (f_{z} - M)

, with

M = \sqrt{1 / λ^{2} - f_{x}^{2} - f_{y}^{2}} > 0

, hence neglecting evanescent components, the calculation can be further simplified:

{\hat{G}}_{A} = \frac{i}{2 M} [δ (f_{z} - M) + δ (f_{z} + M)], {\hat{G}}_{B} = \frac{1}{2 π M} [\frac{1}{(f_{z} - M)} - \frac{1}{(f_{z} + M)}] .

Note the fact that

F_{1}^{- 1} (δ (f_{z})) = 1

and

F_{1}^{- 1} (1 / f_{z}) = i π s g n (z)

. For

f_{x}^{2} + f_{y}^{2} \leq 1 / λ^{2}

, it is interesting to note from Eq. (8) that in the forward direction (

z > 0

), the inverse FTs of

{\hat{G}}_{A}

and

{\hat{G}}_{B}

are identical for

f_{z} > 0

, while they cancel each other for

f_{z} < 0

. The opposite is true in the backward direction

z < 0

, where the components for

f_{z} < 0

reinforce. Since we are interested only in the region

z > 0

,

\hat{G}

for the forward direction [neglecting evanescent components for which

(f_{x}^{2} + f_{y}^{2}) > 1 / λ^{2}

] can be reduced to (

z > 0

)

\hat{G} = {\begin{matrix} i λ δ (f_{r} - 1 / λ) & f_{z} > 0 \\ 0 & f_{z} < 0 \end{matrix} .

Therefore, the 3D transfer function is given by

\hat{K} = 2 π λ f_{z} δ (f_{r} - 1 / λ), f_{z} \geq 0 & z > 0 .

This equation indicates the interesting finding that the 3D diffracted field distribution under RS1 for propagating components can be calculated by first mapping (subject to a weighting factor of

f_{z}

) the 2D-AS of the incident wave to an infinitesimally thin hemispherical cap (radius equal to

1 / λ

), then performing an inverse 3D FT. To illustrate the proposed 3D-AS method numerically, we first consider the propagation of a Gaussian beam in free space. The parameters used are the wavelength

λ = 633 nm

, the radius of the circular aperture

a = 5 mm

, and the beam waist (located at the aperture)

w_{0} = 1.5 mm

. The propagation of the Gaussian beam is computed using Eqs. (4) and (10), and the results are shown in Fig. 2. The inverse 3D FTs involved in the simulation are implemented with the chirp z-transform algorithm [15]. Good agreement with the well-known Gaussian beam solution [16] to the scalar wave equation is observed.

Fig. 1. Geometry for the diffraction calculation from an aperture illuminated from the left.

Download Full Size | PDF

Fig. 2. 3D field distribution of a propagating Gaussian beam. (a) Transverse intensity and (b) phase when $y = z = 0$ ; (c) axial intensity and (d) phase when $x = y = 0$ . The results from the 3D-AS method (blue circles) are compared with the theoretical prediction (red solid curve).

Download Full Size | PDF

To show that the 3D-AS method is also valid for nonparaxial conditions, the volumetric field distribution at the focus of a highly convergent beam is calculated. We set up wavelength $λ = 633 nm$ , aperture $a = 10 λ = 6.33 um$ , and the subtended half-angle $θ = π / 3$ , where $\sin θ$ is the numerical aperture. Figures 3(a) and 3(c) show the cross-sectional intensities in the resultant 3D focal region from the 3D-AS method. The comparison of the 3D-AS method with some existing models are shown [Figs. 3(b) and 3(d)]. For the lateral case, our numerical result matches well with the Debye model [4], which is considered more accurate than the paraxial model [17] given that the N.A. is large at $\sin (π / 3)$ . For the axial case, the Debye model fails, since it produces only a symmetric axial profile. The numerical 3D-AS result matches well with a closed-form expression [18], denoted as the Smith model, that only exists for axial calculations. From the comparison, it can be concluded that the numerical implementation of the 3D-AS method provides good accuracy and is versatile. It can be utilized across various kinds of focusing and propagating scenarios, since it only assumes that evanescent waves are neglected, making it a less limiting case than many of the currently available methods.

Fig. 3. Intensity distributions along the (a) lateral and (c) axial directions for the 3D-AS method, (b) comparison of transverse line intensities ( $y = z = 0$ ) from the 3D-AS method (blue dots), the Debye model (continuous red line), and a paraxial model (continuous brown line). 3D-AS matches well with the more accurate Debye model. (d) Comparison of axial intensities ( $x = y = 0$ ) around the geometrical focus F from the 3D-AS method (blue dots), the Smith model (red line), and the Debye model (brown line). 3D-AS matches well with the more accurate Smith model.

Download Full Size | PDF

In conclusion, we have two main findings in the proposed 3D-AS method. First we demonstrated analytically that the fundamental RS1 can be viewed as a 3D convolution; hence a single 3D-FT establishes the 3D diffracted electromagnetic field in the reciprocal space without approximations. Next we proposed a numerical recipe based on the proposed 3D convolution relationship to evaluate the RS1, under an approximation of neglecting evanescent waves. In this way, the complete 3D volumetric wave propagation field can be first analyzed and subsequently evaluated using the 3D transfer function of a spherical shell in the K-space. Notably, similar derivations can be obtained for other diffraction theories such as the HFP and Rayleigh–Sommerfeld II (RS2). In fact, for propagating components of the diffracted field in the forward direction of $z$ (i.e., $z > 0$ ), RS1, RS2, and HFP can all be treated to a reduced form of a 3D transfer function with the support of a thin spherical shell. Hence the proposed 3D-AS method is universal, in that this direct 3D approach can be applied to interpret and reduce many known diffraction problems. We expect that the proposed methodology will greatly enhance the capacity for 3D computation of diffracted light, whereby many different types of problems such as 3D imaging and beam profiling can be considered.

S. S. Kou and J. Lin are recipients of the Discovery Early Career Researcher Award funded by the Australian Research Council under projects DE120102352 and DE130100954, respectively. J. Lin acknowledges the financial support from the Defence Science Institute, Australia.

References

1. J. W. Goodman, Introduction to Fourier Optics, 2nd ed. (McGraw-Hill, 1996).

2. E. Cuche, P. Marquet, and C. Depeursinge, Appl. Opt. 38, 6994 (1999). [CrossRef]

3. P. Ferraro, S. Grilli, D. Alfieri, S. D. Nicola, A. Finizio, G. Pierattini, B. Javidi, G. Coppola, and V. Striano, Opt. Express 13, 6738 (2005). [CrossRef]

4. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University, 2005).

5. J. A. C. Veerman, J. J. Rusch, and H. P. Urbach, J. Opt. Soc. Am. A 22, 636 (2005). [CrossRef]

6. R. L. Lucke, Eur. J. Phys. 27, 193 (2006). [CrossRef]

7. A. M. Steane and H. N. Rutt, J. Opt. Soc. Am. A 6, 1809 (1989). [CrossRef]

8. C. J. R. Sheppard and M. Hrynevych, J. Opt. Soc. Am. A 9, 274 (1992). [CrossRef]

9. J. Lin, X. C. Yuan, S. S. Kou, C. J. R. Sheppard, O. G. Rodguez-Herrera, and J. C. Dainty, Opt. Lett. 36, 1341 (2011). [CrossRef]

10. C. W. McCutchen, J. Opt. Soc. Am. 54, 240 (1964). [CrossRef]

11. P. Andrés, M. Martinez-Corral, and Ojeda-Castañeda, Opt. Lett. 18, 1290 (1993). [CrossRef]

12. P. C. Clemmow, The Plane Wave Spectrum Representation of Electromagnetic Fields (Pergamon, 1966).

13. R. N. Bracewell, The Fourier Transform and Its Applications, 3rd ed. (McGraw-Hill, 2000).

14. C. J. R. Sheppard, J. Lin, and S. S. Kou, J. Opt. Soc. Am. A 30, 1180 (2013). [CrossRef]

15. A. V. Oppenheim and R. W. Schafer, Digital Signal Processing (Prentice-Hall, 1975).

16. A. E. Siegman, Lasers (University Science, 1986).

17. Y. Li and E. Wolf, J. Opt. Soc. Am. A 1, 801 (1984). [CrossRef]

18. H. Osterberg and L. W. Smith, J. Opt. Soc. Am. 51, 1050 (1961). [CrossRef]

Calculation of the volumetric diffracted field with a three-dimensional convolution: the three-dimensional angular spectrum method

Abstract

References

Cited By

Figures (3)

Equations (10)

Optics Letters