Open Access
Add to BibSonomy
Abstract
We present a surface integral algorithm, utilizing Fourier integrals to solve optical fields within a volume bounded by a complicated polygonal surface. The method enables the full electric field to be solved from electric field values on the bounding surface at any point within the volume. As opposed to FDTD and FEM methods, volume discretization and the need to iteratively solve the E-field at every discrete volume element is not needed with this method. Our new surface integral algorithm circumvents the limitations that exist in current surface methods. Namely, in present methods, the need to determine a Green’s function only allows for simple bounding surfaces, and these methods generally use integrals that cannot utilize computationally fast Fourier integrals. Here, we prove the algorithm mathematically, show it with a numerical example, and outline important cases where the algorithm can be used. These cases include the design of free-form reflectors and near field optical scanning microscopy (SNOM). We then briefly analyze the algorithm's computational scaling.
© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
An efficient algorithm to find the complex electric field distribution of propagating optical modes within a volume bounded by a complicated and, in general, asymmetric surface is of high interest in applied optics. For example, obtaining this efficient algorithm would enable the design of freeform light sources to create a desired lighting effect within the bounded volume [1–5]. As well, the algorithm could efficiently compute optical solutions for holography using spatial-light modulators (SLMs), microelectromechanical mirrors (MEMMs), or freeform lenses on curved surfaces, which has immense applications from experimental optics to industry [6–9].
We present an efficient algorithm to solve the volume Helmholtz equation bounded by arbitrary convex polygonal surfaces. This algorithm avoids the usual bottlenecks for computational efficiency in numerical solvers, such as volume meshes or iterative solvers. All that is needed is the electric field distribution of light along the surface, entering the volume.
Currently, volumetric methods do exist to solve the optical field equation – i.e., the Helmholtz equation, such as finite difference time domain (FDTD) or finite element methods (FEM) [10,11]. However, these methods require that both bounding surfaces and volumes are meshed. The full 3-D meshing must be used, and the fields iteratively calculated across all voxels (3-D pixels) to find the electric field at any point in the volume, which is an exhaustive approach. Moreover, reflections along the bounding surfaces, or the property that light crosses the boundary in both directions of the normal to the surface, yields that additional assumptions that may not reflect the system's physical nature must be placed. An example of an additional assumption is the need to use perfect matched layers (PMLs) to mimic perfect absorbers in the domain outside of the bounded volume. Since the PML occupies space outside the region of interest, the undesired effect of expanding the simulation domain and, therefore, computational resources and time must be present.
Another class of solvers that are more efficient than FDTDs or FEMs are surface integral methods (SIMs). Here, only the electric fields on the bounding surface are used in an integral expression always over the surface mesh elements to calculate the electric field at any point in the enclosed volume. In contrast to FDTDs or FEMs intermediate optical fields at other regions in the volume do not need to be found. Then surface integral methods reduce the computational time/complexity by firstly, only needing the two-dimensional surface pixels, and secondly, only using prior given information in the calculation of the optical field value at a point in the volume. An additional advantage that lies with surface integral methods is that since the iterative need of FDTDs or FEMs is not present, parallelization across many computers can simultaneously evaluate all volume points of interest.
However, state-of-the-art surface integral methods are limited in scope to simple bounding surfaces where a Green’s function must be found [12–16] or to systems where surface currents along the bounding surfaces must be present [12]. The need for surface currents in existing SIMs negates these methods’ applicability for homogenous light propagation (e.g., in air). Additionally, the integral expressions derived with current SIM methods are numerically complex to evaluate, as they do not rely on simple FFT integrals.
The Green’s function's necessary use for current SIMs limits the validity and accuracy of the methods for more general bounding surfaces as many sources of error emerge with the Green function kernel. In general, the Green’s function integral is not over a complete Hilbert basis, but over a sub-space that results in issues with the Green’s approach's generality to solve any surface geometry and field.
The loss of generality manifests itself, for example, in one of the chief practical problems with the kernel, i.e., that components of electromagnetic radiation traveling parallel to convex polygons on the bounding surface cannot be incorporated well within the Green’s function approach. The lack of these wave components emerge because the Green’s function, usually corresponding to a point source, has discontinuities that emerge within its momentum representation, or the lack of evanescent components [17].
The difficulty with incorporating parallel wave components in the Green’s kernel leads to a nonzero electric field outside of the polygonal surface element being considered on the aperture plane. These are redundant artifact fields that lead to inaccurate results. In past approaches [18], for singular planar apertures, a mirror sink term was introduced in the Green’s function to render the fields zero outside the aperture window; however, this ad-hoc addition fails in the general case where a sink term may not be present.
The problem of artifact electric fields outside of the aperture windows is symptomatic of a more general problem of the Green’s function approach; that is, it’s most useful as a far-field diffraction method. In the near-field, the artifact fields are not only present, but the kernel becomes divergent as there is a distance-to-the-aperture-point-source term in the denominator that goes to zero. Thus, for small enclosed volumes or locations in the inner volume that lie close to the bounding surface, the approach can dramatically fail.
To circumvent the current difficulties with the Green’s function approach, new methods rely on obtaining function kernels (i.e., not based on point sources) that incorporate more geometric information of the bounding surfaces [16]. However, these methods rely on the specific problem and are hard to derive, and moreover, since they cannot span the entire Hilbert space of possibilities, they are limited to specific optical field distributions.
Our method presented in this paper is a three dimensional SIM based on the Fourier integral distributions at the bounding surface and can be used without far-field preconditions. Moreover, it can be used for any convex polygonal surface, regardless of the complexity, and is compatible with Fast Fourier transform (FFT) eigenmode solvers. A Fourier based SIM has been explored in two-dimensions for smooth surfaces (i.e., not discrete polygonal surfaces) in [16]. However, this method is still based on the usual Green’s function, albeit solved with a Fourier series approach, so it possesses the full constraints of the Green’s approach. In contrast, our Fourier surface method is not based on surface currents or Green’s functions, and its theoretical formulation spans the full Fourier Hilbert space, without limit to a sub-space of functions (it is only limited in its numerical implementation by the constraints imposed by discretizing to a Fourier series). We will mathematically prove and rigorously describe the scalar version of the algorithm here, leaving the vectorial formulation and numerical verification for future work.
2. Notation and important concepts in Fourier optics used in the method
Our algorithm seeks to solve the linear optical propagation of light within a bounded volume containing no sinks or source terms but can contain perfect scatters. The algorithm solves the well-known homogeneous linear Helmholtz equation given as [11],
with boundary conditions defined on a bounding polygonal surface. The bounding surface must be convex, such that the outer surface normal of surface polygons do not intersect other surface elements.The a priori information needed for the algorithm is the optical electric field distribution at every surface polygon of the optical modes that propagate across the surface aperture towards the inner volume (propagate in the negative ${z_j}$ direction). We note here the advantage of our method that the full electric field distribution then is not necessary, as the fields propagating in the outward ${z_j}$ direction are not needed.
To start the description of our surface integral algorithm, we start by defining relevant terms. The term “aperture” will be used to describe a region of a plane surface that can be oriented in any arbitrary way in space. The polygonal surface enclosing the volume of interest will be composed of such apertures. These apertures may physically correspond to planar polygonal light sources or to transparent surfaces that transmit light. There are two types of coordinate systems used: local and global. The local right-handed orthogonal coordinate system, ${x_j},{y_j},{z_j}$ is defined explicitly for an aperture labeled with index j. ${z_j}$ is the dimension corresponding to the normal vector of the aperture. The outward normal direction (away from the enclosed volume) on surface apertures corresponds to negative ${z_j}$. A visual description of the bounding surface composed of apertures is given in Fig. 1.
Fig. 1. a) An example schematic of the bounding convex surface. The local coordinate system is shown for an aperture labeled as j. Also shown is the general reference coordinate system. The electric field penetrating the aperture from outside to the inner volume is denoted in the above example. b) A reflector arrangement that can be solved by the algorithm. The reflector is composed of a transparent entrance aperture for the incoming optical field. All fields are reflected along the inner reflector surface denoted in gray and exit through the entrance aperture. The algorithm provides a solution to the reflected optical field coming out from the reflector (i.e., $u_1^{0 - }$) as well as the internal reflected field in the volume.
The aperture transmission function is generally defined as:
The optical field distribution on a given aperture is a function of ${x_j},{y_j}$ and is labeled as ${u_j}({{x_j},{y_j},{z_j} = 0} )$ or $u_j^0$ in short form. In terms of an arbitrary incident optical field on the entire aperture plane, ${u_{inc}}({{x_j},{y_j},{z_j} = 0} )$, represented in the local coordinate system,
For simplicity, all-optical fields are monochromatic at the same angular frequency, $\omega $. The time-dependent phase, ${e^{ - i\omega t}}$, is not explicitly shown in the electric field representation. The results in this paper can easily be extended to the polychromatic case. There, this methodology would solve for the spatial distribution per frequency.Following Fourier optics [19], the monochromatic optical field with boundary condition ${u_j}$ (i.e., the value of the aperture optical at ${z_j} \ne 0$) is given and fully described at any point in three-dimensional space as (Fourier integral constants in front of all Fourier integrals are omitted for clarity.):
Where, $k_x^j,k_y^j$ are the angular wavenumbers (spatial frequencies) corresponding to the local coordinate variables ${x_j},{y_j}$. $k_z^j$ is the component of the wavenumber along the normal of the aperture and can be real or imaginary. Using the dispersion relationship of monochromatic light in isotropic homogeneous medium,
From Eq. (5), it can be seen that $k_z^j$ can be either positive or negative, thus, $u_j^0$ can be composed of the summation of two counter-propagating fields labeled as $u_j^{0 + }$ for total field that propagates toward the inside normal direction (+z so $k_z^j > 0$) and $u_j^{0 - }$ for the outside normal direction (-z so $k_z^j < 0$).
Then,
As will be seen in the algorithm’s proof, only the inward propagating fields are needed, so we omit $u_j^{0 - }$ dependent terms from the summation integrals. For evanescent waves when $k_z^j$ is imaginary, the wave is propagating only in the two-dimensional plane spanned by ${x_j},{y_j}$ and exponentially decays as a function of ${z_j}$ (these imaginary values can occur both in $u_j^{0 + },\; u_j^{0 - }$) [19].
To describe our method, one aperture field will be evaluated at another aperture (i.e., it will act as part of the incident field on the other aperture). Through Eqs. (4)–(6), the field values of the field emanating from the first aperture, labeled as ${u_p}({{x_m},{y_m},{z_m} = 0} )$ in Eq. (5), will be found on the second aperture and represented as a function of the second aperture's local coordinate basis. Notation-wise, this will be represented as:
$u_{p@m}^0$ is a summation: $u_{p@m}^0 = u_{p@m}^{0 + } + u_{p@m}^{0 - }$ where the superscript $+ $ denotes the incident field propagating in the positive ${z_m}$ direction and the ${-}$ denotes the incident field propagating in the negative ${z_m}$ direction. These fields are obtained by breaking the integral in Eq. (4) for ${u_p}$, e.g., for $u_{p@m}^{0 + }$ as follows:When $k_x^p,k_y^p\epsilon R|({k_x^p\widehat {{{\boldsymbol x}_{\boldsymbol p}}} + k_y^p\widehat {{{\boldsymbol y}_{\boldsymbol p}}} + k_z^p\widehat {{{\boldsymbol z}_{\boldsymbol p}}}} )\cdot \widehat {{{\boldsymbol z}_{\boldsymbol m}}} = 0$, these fields are grouped with the inward to outward propagating fields, such that they are omitted from consideration for the algorithm (described in the next section). Once the positive and negative incident fields on aperture m are found, $u_{p@m}^0$ can be propagated to any arbitrary point (in the coordinate system of aperture $m$) by applying Eq. (4).
3. Scalar version of the new bounding surface method
Having defined all relevant notation and needed concepts in Fourier optics, we now turn our attention to the derivation of the general method, presented as the following algorithmic steps:
- 1. Define the bounding surface by a convex polygonal meshing. Each polygon is a surface aperture.
- 2. Arbitrarily number each surface polygon starting at $j = 1$.
- 3. Define a series of fields, incrementing with j and terminating at the last surface aperture ($j = F)$, as such:
These $\rho $’s are called the reduced fields of the corresponding apertures.
- 4. Then, the series of the above fields give u within the volume of the bounding surface, $V$, i.e.,$$u({{x_g},{y_g},{z_g}} )= \mathop \sum \nolimits_{j = 1}^{j = F} {\rho _j}({{x_g},{y_g},{z_g}} )- \mathop \sum \nolimits_{j = 1}^{j = F} \mathop \sum \nolimits_{m = 1}^{m = j - 1} {\rho _{j@m}}({{x_g},{y_g},{z_g}} ), $$where, ${x_g},{y_g},{z_g}$ are the coordinates in a predefined global coordinate system. In general, the double summation can be omitted as its summands have negligible absolute values compared to the summands of the single summation for coordinate points on the bounding surface the E-field equates only to the single summation term.
The procedure outlined in steps 1-4 are reminiscent of the Gram-Schmidt procedure to obtain an orthogonal set of basis vectors. Since the order of how the Gram-Schmidt procedure is carried out influences the numerical round-off error [20], it is then expected that the order of how the summations of Eq. (11) are evaluated also influences this error. Thus, a modified procedure following how subtraction is ordered in the ‘modified’ Gram-Schmidt algorithm [20], which produces less numerical round-off error, would be to evaluate the summation by subtracting the reduced field of a surface aperture from all other apertures “at once” when it is found, and not in an iterative fashion.
The exact numerical discretization needed to propagate surface aperture fields to accomplish steps 3, 4, and maintain a manageable error after necessary coordinate grid transformations will not be discussed in detail here. In this regard, general methodologies must be developed for this method and will utilize already known techniques for maintaining sampling conditions for the Fourier propagation integral at different imaging planes [21,22].
In brief, the numerical discretization must, at the least, satisfy the Nyquist conditions for the two-dimensional FFT Ref. [22]. Given a monochromatic field of frequency,f, and associated magnitude of wavenumber, ${k_o} = \frac{f}{c}$, the Nyquist conditions become:
Where, $\Delta x,y$ are the sampling step-distance of the aperture coordinates, $\Delta {k_{x,y}}$ is the momentum step-interval, ${|{x,y} |_{max}}$ is the maximal distance in the coordinate in the bounding surface aperture set. If these conditions are followed at every aperture, undersampling that could result from propagating one aperture field to another is avoided since one aperture's coordinate interval will be mapped to a more considerable or equivalent distance on another aperture plane.
3.1 Proof of scalar method
The proof of our surface integral method is best transmitted through the following thought-experiment. In the thought experiment, we place a perfect absorber directly underneath an aperture, i.e., on the side facing the inner volume—So that the light passing through the aperture is still described by Eq. (3). The absorber placement is simulated by mimicking what the absorber’s effects would be, i.e., by subtracting the corresponding aperture mode propagating towards the inner volume, from the internal volume and other surface elements. Absorbers are placed iteratively until all apertures are covered. When all apertures are covered, no internal fields are present within the volume, i.e., the electric field is zero everywhere, since there are no optical source or sink terms within the volume. Thus, all aperture mode fields subtracted due to the placement of the perfect absorbers are equal to the original electric field within the volume. It will be shown that this is equivalent to Eq. (11). Thus, we prove that Eq. (11) is equal to the electric field within the volume given the bounding surface boundary conditions.
To proceed with the outlined proof approach, we start by placing a perfect absorber on aperture $J = 1$. Then, the electric field within the volume becomes,
- 1. $u - {\rho _1}$
while now placing the next perfect absorber on aperture $J = 2$, yields that
- 2. $u - {\rho _1} - {\rho _2} + {\rho _{2@1}}.$
The term ${\rho _{2@1}}$, i.e., the inward propagating optical mode from ${\rho _2}$ crossing aperture 1 emerges in the summation because the perfect absorber at $J = 1$ also blocks any field impingement from aperture $J = 2$ on $J = 1$. However, in general, this term can be omitted as this contribution is negligible in comparison to ${\rho _2}$.
Continuing by placing a perfect absorber on aperture $J = 3$, yields,
- 3. $u - {\rho _1} - {\rho _2} + {\rho _{2@1}} - {\rho _3} + {\rho _{3@2}} + {\rho _{3@1}}$
and placing perfect absorbers until aperture F then yields:
- 4. $\begin{aligned} u - \mathop \sum \nolimits_{j = 1}^{j = F} {\rho _j} &+ \mathop \sum \nolimits_{j = 1}^{j = F} \mathop \sum \nolimits_{m = 1}^{m = j - 1} {\rho _{j@m}} = 0\; \forall \; ({{x_g},{y_g},{z_g}} )\in V\; \\ &\therefore u = \mathop \sum \nolimits_{j = 1}^{j = F} {\rho _j} - \mathop \sum \nolimits_{j = 1}^{j = F} \mathop \sum \nolimits_{m = 1}^{m = j - 1} {\rho _{j@m}}\forall \; ({{x_g},{y_g},{z_g}} )\in V\; \; \end{aligned}$
We complete the proof by noting that the convex bounding surface condition is necessary to reduce the algorithm's complexity. The complexity is reduced because the convex property that only enables electric fields propagating towards the inner volume at each surface aperture is considered Other directions of propagation can then be omitted from the analysis.
The proof reveals that electric fields’ directionality can reduce the amount of surface elements of the algorithm at each surface aperture. Exactly one surface aperture with no electric field propagating across it towards the inside volume of the bounding surface (as boundary condition) can be omitted from the set of surface apertures considered for the algorithm. Omitting this surface aperture is allowed because, if perfect absorbers cover all other apertures, there would still be no fields in the inner volume; thus, step 4 in the algorithm is still valid. Moreover, since the field across the omitted surface aperture would also be zero everywhere, Eq. (11) would also be valid to describe the field across this omitted surface aperture.
4. Numerical example
While the algorithm was proved in Section 3, it is prudent to show the practical implementation of the algorithm in a numerical example on the MATLAB platform, using the built-in FFT optimization. To illustrate the algorithm we choose to simulate a plane wave propagating, such that it passes through a square box surface, of side length a (for our simulation $a = 17{k_o}$), aligned to have the input face be perpendicular to the phase fronts, as seen in Fig. 3. This example, over a plane-wave eigenmode of the corresponding system, is representative of any solution field here, which would be a summation over plane-wave eigenmodes.
We now solve for the aperture field for aperture 6 to verify that it matches the true solution. We can choose this aperture as the “unknown” aperture of E-field distribution to be solved by our algorithm since no electric field goes from the outside to the inside of aperture 6. The true solution for aperture 6 would be a constant intensity across the full aperture, as it is aligned to be perpendicular to the plane wave impinging on it.
To start the numerical solution, we note that the original boundary fields, i.e., the electric field of wave components going from the outside to the inside of apertures 2-6 are zero, while the electric field going from the outside to the inside of aperture 1 is a constant value (we choose as 1) across the entire aperture. Since all apertures are square, including aperture 1,
Using ${\rho _1}$, we then numerically solve for all other fields in step 1-4, whereby the sampling conditions are given by Eqs. (12)–(13). Here we note that, since all apertures have the same square aperture function, the grid spacing remains the same. Thus, it is sufficient only to set the “z” component coordinate spacing of a given aperture to be the same as that of the aperture grid. We also use the rectangular symmetry to reduce the numerical load of the algorithm by setting aperture 4 to the mirror inversion of aperture 3 and the field from aperture 1 to aperture 5 as the mirror inversion of that of aperture 2. Finally the field from aperture 1 to aperture 3 is the field from aperture 2 rotated to match the coordinate system of aperture 3.
The E-field result on aperture 6 from the algorithm is shown in momentum space and in real space in Figs. 2(b) and 2(c) (across the $k_x^6$ coordinate), respectively. The dotted line curve in Fig. 2(c), shows the ideal $k_x^6$ momentum distribution of the true solution. The field profile remains close to constant across the entire aperture, deviating slightly, especially at the edges (maximal deviation 18% of average value), due to aliasing and sampling errors. The simulation is done on a standard desktop laptop, to show the versatility of the algorithm, so sampling constraints, resulting in the error, are restricted by limitations imposed by the computer.
Fig. 2. a) Representative numerical test system of the bounding surface method. b) Normalized to peak momentum spectrum along the x coordinate. Dotted line represents the true solution which matches the simulation with an almost perfect match. c) real-space electric-field amplitude distribution over aperture 6.
5. Application examples of the algorithm
The omission of a surface aperture with only outward propagating optical fields expands the scope of problems our algorithm can be applied to. For example, if the givenness of the problem under consideration is that one polygonal surface element of the bounding surface has an unknown outward propagating electric field, then the algorithm can solve for the field across this unknown surface element as well. When designing intricate lighting patterns, this is of particular use to illuminate a “ground” or sample surface, which finds application from room lighting to designing microscope illumination on sample planes. These applications are of great importance to industry and optical engineering. Another application is to obtain an incident field on a sample object in diffractive imaging such as in x-ray free-electron lasers (XFELs) [23] or on near field scanning optical microscope (NSOM) fiber tips. Solving for optical fields across a sample object in diffractive microscopy is essential to the field, and obtaining the field distributions at the output of an NSOM tip is at the heart of near field microscopy.
We illustrate the algorithm's high potential use in near field microscopy through a conceptual example. In this example, the unknown aperture is the base of the NSOM tip-, and light enters the NSOM tip through the remaining boundaries from the sample [24–26]. Figure 3 illustrates this concept, where a pyramidal NSOM tip is used to image light emitted from red fluorophores. The light propagating out of the tip from the aperture base surface, labeled as $u_1^{0 - }$, is used to reconstruct the sample image and can be found by our algorithm. The base field information, such as its spatial mode content, could effectively design downstream optics in the NSOM apparatus, such as zone plates, transporting the field in fiber optic cabling, and validating experimental findings.
Fig. 3. A conceptual example of an NSOM tip imaging a fluorophore sample. Our algorithm can find the exiting field on the base aperture of the pyramidal tip, given known fields across other aperture elements. The base field information, such as its spatial mode content, could be used to effectively design downstream optics in the NSOM apparatus, such as zone plates, the transport of the field in fiber optic cabling, and to validate experimental findings.
Another application centers on reflector surfaces, where the boundary conditions represent light that is completely reflected from the surfaces and then outputted through one transparent surface aperture, shown in Fig. 1(b). All fields are reflected along the inner reflector surface. The reflected field can only exit the input aperture (labeled as aperture 1). Thus, for fields u, aperture 1 only has inner to outer propagating optical fields and can be omitted from the algorithm's requirements. The algorithm provides as a solution the reflected optical field coming out from the reflector (i.e., $u_1^{0 - }$) as well as the internal reflected field in the volume. The output aperture also acts as the input aperture for the incident field on the reflector; however, this input field can be omitted from being considered without impact on the results, provided the reflected optical fields are given at the reflective aperture surfaces. The algorithm’s use in reflector design can be of particular importance in the design of non-planar MEMMS setups, where light undergoes reflections on a complicated convex surface, or in the design of directional retroreflectors for night-time driving [27].
6. Amount of computational operations needed
Our surface integral algorithm's theoretical scaling reveals that it scales with fewer operations than the corresponding volume discretized FDTD and FEM methods. To start counting the required operations, we start first by the subtraction that must take place from the propagated solutions of previous surface boundary apertures to a given considered surface boundary aperture (i.e., step 1 in the algorithm). A Fourier transform to the momentum domain will have to be performed to propagate the reduced solution from a previous surface aperture. Let ${N_j}$ represents the number of grid points at a surface aperture j. Then the Fourier transform to the momentum domain would be: $\sigma \left[ {\sqrt {{N_j}} \textrm{log}({{N_j}} )} \right]$. The momentum spectrum would then have to be multiplied by the propagator, and the two-dimensional Fourier inverse transform would have to be performed at each grid point on the next surface. These points are at a different coordinate system and thus, the double integral may (in the worst case) have to be evaluated without using the fast Fourier transform technique (FFT) requiring $\sigma [{{N_j}} ]$.
Using these facts, the total amount of operations ($S$) would be:
If the number of grid points are the same, the expression for the total amount of operations becomes:
In general, ${F^2} \ll \sqrt {{N_j}} $, since the amount of surface apertures do not need to be varied for the convergence of the algorithm. However, the discretization of each surface element, represented by ${N_j}$ must be a large number for the numerical FFTs to converge in the algorithm. Then the algorithm scales less than ${N_j}^{\frac{3}{2}}$ which corresponds to the scaling of volume element methods. Furthermore, the iterative need of FDTDs or FEMs on volume discretization is not present in our method, so parallelization across many computers can simultaneously evaluate all volume points of interest, significantly decreasing the required simulation time.
7. Conclusions
A new semi-analytic numerical algorithm and solution to the inverse Helmholtz problem has been developed in this paper. From the electric field's values on a bounding surface, the field values at any point within the bounded volume can be found by our surface integral method. Since the discretization only takes place over surfaces, our algorithm scales less than volume element methods, such as FDTD and FEM, and requires less time to implement. Moreover, once found, the result of the algorithm, a semi-analytic expression describing the electric field, can be used to find values without the need to iteratively solve for all volume elements (as opposed to FDTD and FEM methods). By treating the electric field on the bounding surfaces as integrals over Fourier eigenmodes in the fashion described here, the computation of how light acts on reflective and transmissive boundaries in complex geometries such as NSOM tips can be found efficiently using this method. We have also demonstrated this algorithm with a pertinent numerical example using a cubic convex bounding surface. The algorithm can be expanded to a vectorial version, which will be explored in future work.
Disclosures
The author declares no conflicts of interest.
References
1. 1. J. B. Jiang, C. F. Cheung, S. To, K. W. Cheng, H. Wang, and W. B. Lee, “Design and fabrication of freeform reflector for automotive headlamp,” 2006 2nd International Conference on Power Electronics Systems and Applications, Hong Kong, 2006, pp. 220–224.
2. K. P. Thompson, P. Benítez, and J. P. Rolland, “Freeform Optical Surfaces: Report from OSA’s First Incubator Meeting,” Opt. Photonics News 23(9), 32–37 (2012). [CrossRef]
3. . Top of Form
4. B. G. Assefa, T. Saastamoinen, M. Pekkarinen, V. Nissinen, J. Biskop, M. Kuittinen, J. Turunen, and J. Saarinen, “Realizing freeform lenses using an optics 3D-printer for industrial based tailored irradiance distribution,” OSA Continuum 2(3), 690–702 (2019). [CrossRef]
5. J. C. Miñano, P. Benítez, and A. Santamaría, “Free-form optics for illumination,” Opt. Rev. 16(2), 99–102 (2009). [CrossRef]
6. M. Khorasaninejad and F. Capasso, “Metalenses: Versatile multifunctional photonic components,” Science 358(6367), eaam8100 (2017). [CrossRef]
7. F. Chen, S. Liu, K. Wang, Z. Liu, and X. Luo, “Free-form lenses for high illumination quality light-emitting diode MR16 lamps,” Opt. Eng. 48(12), 123002 (2009). [CrossRef]
8. H. Joonku, H. Kim, Y. Lim, G. Park, and B. Lee, “Wide viewing angle dynamic holographic stereogram with a curved array of spatial light modulators,” Opt. Express 16(16), 12372–12386 (2008). [CrossRef]
9. B. Parkyn and D. Pelka, “Free-form illumination lenses designed by a pseudo-rectangular lawnmower algorithm,” Nonimaging Optics and Efficient Illumination Systems III. Vol. 6338. International Society for Optics and Photonics, 2006.
10. S. D. Gedney, “Introduction to the finite-difference time-domain (FDTD) method for electromagnetics,” Synth. Lect. Comp. Electromag. 6(1), 1–250 (2011). [CrossRef]
11. I. Babuška, F. Ihlenburg, E. T. Paik, and S. A. Sauter, “A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution,” Comp. Meth. Appl. Mech. Eng. 128(3-4), 325–359 (1995). [CrossRef]
12. S. R. Arridge and M. Schweiger, “Inverse methods for Optical Tomography,” In: H. H. Barrett and A. F. Gmitro, eds. Information Processing in Medical Imaging. IPMI 1993. Lecture Notes in Computer Science, vol 687. Springer, Berlin, Heidelberg, (1993).
13. S. Kronberg, M. J. H. Anthonissen, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Modelling Light Scattering in the Context of Freeform Optical Design,” EPJ Web of Conferences. Vol. 238. EDP Sciences, 2020.
14. J. H. M. ten Thije Boonkkamp and W. L. IJzerman, “Illumination freeform design using Monge-Ampère equations,” in Optical Design and Fabrication 2019 (Freeform, OFT), OSA Technical Digest (Optical Society of America, 2019), paper FT2B.1.
15. N. Guan, S. Habu, K. Takenaga, K. Himeno, and A. Wada, “Boundary element method for analysis of holey optical fibers,” J. Lightwave Technol. 21(8), 1787–1792 (2003). [CrossRef]
16. R. M. James, “On the use of Fourier series/FFT's as global basis functions in the solution of boundary integral equations for EM scattering,” IEEE Trans. Antennas Propag. 42(9), 1309–1316 (1994). [CrossRef]
17. C. J. Sheppard, S. S. Kou, and J. Lin, “The Green-function transform and wave propagation,” Front. Phys. 2, 67 (2014). [CrossRef]
18. J. Saatsi and P. Vickers, “Miraculous success? Inconsistency and untruth in Kirchhoff’s diffraction theory,” British J. Phil. Sci. 62(1), 29–46 (2011). [CrossRef]
19. E. G. Steward, Fourier optics: an introduction. Courier Corporation, 2004.
20. W. Hoffmann, “Iterative algorithms for Gram-Schmidt orthogonalization,” Computing 41(4), 335–348 (1989). [CrossRef]
21. H. Zia, “Simulation of white light generation and near light bullets using a novel numerical technique,” CNSNS 54, 356–376 (2018). [CrossRef]
22. S. Mehrabkhani and T. Schneider, “Is the Rayleigh-Sommerfeld diffraction always an exact reference for high speed diffraction algorithms?” Opt. Express 25(24), 30229–30240 (2017). [CrossRef]
23. A. P. Mancuso, O. M. Yefanov, and I. A. Vartanyants, “Coherent diffractive imaging of biological samples at synchrotron and free electron laser facilities,” J. Biotechnol. 149(4), 229–237 (2010). [CrossRef]
24. P. Bazylewski, S. Ezugwu, and G. Fanchini, “A review of three-dimensional scanning near-field optical microscopy (3D-SNOM) and its applications in nanoscale light management,” Appl. Sci. 7(10), 973 (2017). [CrossRef]
25. A. El Eter, N. M. Hameed, F. I. Baida, R. Salut, C. Filiatre, D. Nedeljkovic, E. Atie, S. Bole, and T. Grosjean, “Fiber-integrated optical nano-tweezer based on a bowtie-aperture nano-antenna at the apex of a SNOM tip,” Opt. Express 22(8), 10072–10080 (2014). [CrossRef]
26. A. Casandruc, G. Kassier, H. Zia, R. Bücker, and R. D. Miller, “Fiber tip-based electron source,” J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom. 33(3), 03C101 (2015). [CrossRef]
27. A. Rajib, A. A. Rifat, M. U. Hassan, A. K. Yetisen, and H. Butt, “Phase-conjugated directional diffraction from a retroreflector array hologram,” RSC Adv. 7(41), 25657–25664 (2017). [CrossRef]
