Green Scalar Function Method for Analyzing Dielectric Media

1. Introduction

Scattering problems are of great importance in applications such as optical device design [1], medical imaging [2], object detection and materials science. Theories such as the Mie theory [3], the Born approximation [4] and advanced numerical methods are essential tools for tackling these problems [5,6,7]. Some of these methods consists on solving time-dependent Maxwell’s equations [8], obtaining Mie’s coefficients [9], employing method of moments (MoM) or Galerkin Method [10] . However, these theories have certain disadvantages such as: complex and computationally expensive calculations, application to regular surfaces, lack of accuracy with high dielectric constant or small dimensions, hard dicretization of the equations, etc. [11,12]. Therefore, in the present work a scalar formalism with Green’s functions [13,14] is shown. These show the impulsional response of a linear system. Specifically, in the electromagnetic domain, an impulsional response represents a point source in space as if it were a Dirac delta. This makes it possible to describe the electric and magnetic fields and their potentials in terms of these point sources. This allows us to know the distribution of the reflected, refracted and scattered fields. On the other hand, the scalar formulation allows us greater computational efficiency and the possibility of contrasting this methodology with different scattering cases where an analytical solution is presented. To this end, the Helmholtz equation will be derived from Maxwell’s equations and the solution of this equation will be posed in terms of the scalar Green’s functions. In some cases, due to the symmetry of the problem, the dimensionality can be simplified. In the following, different geometries will be presented in order to analyse the resulting electromagnetic scattering, varying the system’s own parameters such as dimensions and refractive index.

2. Theoretical Background

The Green’s function applied to the Helmholtz equation is the solution of the wave equation for a point source. Thus if the solution for this point source is known, then the solution for a general source can be obtained by the superposition principle. This result follows from the linearity of the wave equation itself, since any source can be described by the linear superposition of point sources [15]. For our case we will consider a scalar field $E\left(\mathit{r}\right)$ which has a source $s\left(\mathit{r}\right)$. The Helmholtz equation associated with this system would be

$$({\nabla }^2+k^2)E\left(\mathit{r}\right)=s\left(\mathit{r}\right)$$

(1)

And using the definition of Green’s scalar function, now renamed $g(\mathit{r},{\mathit{r}}^{\prime })$ we have that

$$\left({\nabla }^2+k^2\right)g(\mathit{r},{\mathit{r}}^{\prime })=\delta (\mathit{r}-{\mathit{r}}^{\prime })$$

(2)

Thus, we can establish a relationship between (1) and (2) as

$$E\left(\mathit{r}\right)=-{\int }_{V}g(\mathit{r},{\mathit{r}}^{\prime })s\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }$$

(3)

The solution of the scalar Green’s function in 3D is obtained in Appendix A. This solution will only depend on $R=|\mathit{r}-{\mathit{r}}^{\prime }|$ and in general could be written as

$$g_{3\mathrm{D}}(\mathit{r},{\mathit{r}}^{\prime })=g(\mathit{r}-{\mathit{r}}^{\prime })={\displaystyle \frac{e^{ik|\mathit{r}-{\mathit{r}}^{\prime }|}}{4\pi |\mathit{r}-{\mathit{r}}^{\prime }|}}={\displaystyle \frac{e^{ikR}}{4\pi R}}$$

(4)

Given the equations (3) and (4), we can write the scalar field associated with such a source as

$$E\left(\mathit{r}\right)=-{\int }_V{\displaystyle \frac{e^{ik|\mathit{r}-{\mathit{r}}^{\prime }|}}{4\pi |\mathit{r}-{\mathit{r}}^{\prime }|}}s\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }$$

(5)

2.1. 2D Green Function

As in 3D systems, the scalar Green’s function in 2D geometries will only depend on the relative position between the source and the observer, so it is convenient to define the relative coordinates.

$$\mathit{r}-{\mathit{r}}^{\prime }\equiv \rho -{\rho }^{\prime }=\left(\rho cos\varphi ,\rho sin\varphi ,0\right)=(x-x^{\prime },y-y^{\prime },z-z^{\prime })$$

(6)

Thus, reducing the three-dimensional problem to one in the plane. In the 3D case the scalar Green’s function refers to a point source of electromagnetic fields, located in $(x^{\prime },y^{\prime },z^{\prime })=(0,0,z^{\prime })$. However, the 2D case refers to the field generated by an infinite linear charge source observed in a plane $z=\mathrm{cte}$ whose axial dependence for with the term $exp\left(i{k}_{z}z\right)$. So, one can obtain this 2D Green’s function, $g_{2\mathrm{D}}(\mathit{r},{\mathit{r}}^{\prime })$ as the line integral of this linear source

$$\begin{array}{c}{\displaystyle g_{2\mathrm{D}}(\mathit{r},{\mathit{r}}^{\prime })={\int }_{-\infty }^{\infty }{g}_{3\mathrm{D}}(\mathit{r},{\mathit{r}}^{\prime })e^{i{k}_z{z}^{\prime }}d{z}^{\prime }={\int }_{-\infty }^{\infty }{\displaystyle \frac{e^{ik\sqrt{x^2+y^2+{(z-z^{\prime })}^2}}}{4\pi \sqrt{x^2+y^2+{(z-z^{\prime })}^2}}}{e}^{i{k}_z{z}^{\prime }}d{z}^{\prime }=-{\displaystyle \frac{i}{4}}{H}_0\left(k_{\rho }\rho \right)e^{i{k}_{z}z}}\end{array}$$

(7)

where $H_0\left(k_{\rho }\rho \right)$ denotes the Hankel function of first kind and $k=\left|\mathit{k}\right|=\left|\right({\mathit{k}}_{\rho },k_z\left)\right|$ is the wave vector modulus and $k_{\rho }=\left|{\mathit{k}}_{\rho }\right|$ is the transverse wave vector modulus, respectively.

2.2. Discretization of the Equations

Consider the Helmholtz scalar equation applied to the electric field and applied to a Green scalar function

$${\nabla }^{2}E\left(\mathit{r}\right)+k_0^2{\epsilon }_{B}E\left(\mathit{r}\right)=-k_0^2\Delta \epsilon \left(\mathit{r}\right)E\left(\mathit{r}\right)$$

(8)

$$\left({\nabla }^2+k_0^2{\epsilon }_B\right)g(\mathit{r},{\mathit{r}}^{\prime })=-\delta (\mathit{r},{\mathit{r}}^{\prime })$$

(9)

Where, $\Delta \epsilon \left(\mathit{r}\right)=\epsilon \left(\mathit{r}\right)-{\epsilon }_B$ is the dielectric contrast, indicating the difference in dielectric permittivity between the object $\epsilon \left(\mathit{r}\right)$ and the background ${\epsilon }_B$. The procedure followed to obtain the scattered pattern is analogous to the one described in [16], except that we address the scalar case. On the other hand, if we introduce equation (9) into equation (8), electric field $E\left(\mathit{r}\right)$ can be obtained by

$$E\left(\mathit{r}\right)=E^{\left(1\right)}\left(\mathit{r}\right)+{\int }_{V}g(\mathit{r},{\mathit{r}}^{\prime })·k_0^2\Delta \epsilon \left({\mathit{r}}^{\prime }\right)E\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }$$

(10)

Where the integration is calculated over the whole volume of the chosen simulated region and $E^{\left(1\right)}\left(\mathit{r}\right)$ refers to the incident electric field. The circumstance to be solved that may arise in the study of a reduced volume is the singularity $|\mathit{r}-{\mathit{r}}^{\prime }|\to 0$. Therefore, when $\mathit{r}\approx {\mathit{r}}^{\prime }$, the singularity can be avoided if we exclude that inifinitesimal element of volume $\delta V$ where this singularity is found. We can then treat this singularity separately, rewriting the equation (10) as

$$E\left(\mathit{r}\right)=E^{\left(1\right)}\left(\mathit{r}\right)+\underset{\delta V\to 0}{lim}{\int }_{V-\delta V}g(\mathit{r},{\mathit{r}}^{\prime })·k_0^2\Delta \epsilon \left({\mathit{r}}^{\prime }\right)E\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }-L·{\displaystyle \frac{\Delta \epsilon \left(\mathit{r}\right)}{{\epsilon }_B}}E\left(\mathit{r}\right)$$

(11)

Where L is a scalar term derived from the dyad $\mathit{L}$[16] which depends on the geometry of the excluded volume $\delta V$. Thus, in our work the observation point $\mathit{r}$ is located outside the dispersive medium or “scatterer‘’, whose position is given by ${\mathit{r}}^{\prime }$. To solve equation (11) numerically, we must first define a grating of N nodes which represents the system. Each node i will be centred in a position ${\mathit{r}}_i$ and will present a volume $V_i$ if we are in 3D or an area in 2D and a value of dielectric contrast. Following, electric field $E_i=E\left({\mathit{r}}_i\right)$ must be discretized as well as scalar Green’s function $g_{i,j}=g({\mathit{r}}_i,{\mathit{r}}_j)$.

$$E_i=E_i^{\left(1\right)}+\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{g}_{i,j}·k_0^2\Delta {\epsilon }_j{E}_j{V}_j+M_i·k_0^2\Delta {\epsilon }_i{E}_i-L·{\displaystyle \frac{\Delta {\epsilon }_i}{{\epsilon }_B}}{E}_{i}i=1,\cdots ,N$$

(12)

and

$$M_i=\underset{\delta V\to 0}{lim}{\int }_{V_i-\delta V}g({\mathit{r}}_i,{\mathit{r}}^{\prime })d{\mathit{r}}^{\prime }$$

(13)

Throughout the implementation of the codes and simulations it has become clear that the terms where the scalar Green’s function or the geometric factor L appears, these prevail over $M_i$. However, many authors for more complex geometries than the ones we will deal with in the project show that the $M_i$ term is more relevant than the L term $E_i$ as in [16]. The equations (12) represent a linear system of equations that can be solved with simple numerical methods. Expressing the system of equations in the form

$$\begin{array}{cccc}\hfill E_i^{\left(1\right)}& =E_i·\left[1-M_i{k}_0^2\Delta {\epsilon }_i+L{\displaystyle \frac{\Delta {\epsilon }_i}{{\epsilon }_B}}\right]-\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{g}_{i,j}·k_0^2\Delta {\epsilon }_j{E}_j{V}_j{E}_j\hfill & \hfill E_i& =A_G^{-1}{E}_i^{\left(1\right)}\hfill \end{array}$$

(14)

where $A_G$ is the interaction matrix which contains the scalar Green’s function and the dielectric contrast and geometric factor terms of the source L, defined by:

$$A_G(i,j)=\left\{\begin{array}{ccc}{\displaystyle 1-M_i{k}_0^2\Delta {\epsilon }_i+L\frac{\Delta {\epsilon }_i}{{\epsilon }_B}}\hfill & \mathrm{if}\hfill & i=j\hfill \\ \\ -g_{i,j}{k}_0^2\Delta {\epsilon }_j{V}_j\hfill & \mathrm{if}\hfill & i\ne j\hfill \end{array}\right.$$

(15)

Following our notation; $E_i$ is the simulated electric field matrix, and $E_i^{\left(1\right)}$ is the incident field.

2.3. Scalar Analytical Solution

The scalar approximation simplifies the analysis of electromagnetic wave scattering by reducing the problem to a scalar field, rather than dealing with the full vector nature of electromagnetic fields. This approach is especially useful when the variations in the wave field occur mainly in a single component, allowing for simpler mathematical handling and greater computational efficiency. By using scalar fields, we can take advantage of well-established mathematical techniques, such as Fourier series, to express the fields in terms of their angular harmonics. In order to compare the numerical results, analytical solution present in [17] are used. In order to achieve this, an homogeneous cylinder is considered to obtain the scattered field in cylindrical coordinates in terms of a Fourier series in the angular coordinate. Lets us consider an incident wave front in plane $y=0$ parallel polarized respect z-axis, thus the electric field in cylindrical coordinates expressed as a Fourier series in the angular coordinate is given by:

$$E^{\left(1\right)}(\rho ,\varphi )=E_0{e}^{ikx}=E_0\sum _{m=-\infty }^{\infty }{i}^m{J}_m\left(k\rho \right)e^{im\varphi }$$

(16)

where $E_0$ is the amplitude of the incident wave, $\varphi$ is the angle that the scattered wave forms with respect to the incident wave and J is the Bessel function of first kind. Once we have expressed the incident field in Fourier terms as a function of the angular harmonics, we can give the expression for the scattered field as

$$E^{\left(3\right)}(\rho ,\varphi )=E_0{e}^{ik\rho cos\varphi }=E_0\sum _{m=-\infty }^{\infty }{i}^m{\alpha }_m{H}_m\left(k\rho \right)e^{im\varphi }$$

(17)

This Equation (17) satisfies the Sommerfeld radiation condition at infinity. Where ${\alpha }_m$ is a constant describing the continuity of the dielectric field derived by the boundary conditions as is is obtained in [18] for a perpendicular polarization across the $xz$-plane, is given by

$${\alpha }_m={\displaystyle \frac{n{J}_m^{\prime }\left(ka\right)J_m\left(nka\right)-J_m\left(ka\right)J_m^{\prime }\left(nka\right)}{n{J}_m\left(nka\right)H_m^{\prime }\left(ka\right)-J_m^{\prime }\left(kna\right)H_m\left(ka\right)}}$$

(18)

where $k=2\pi /\lambda$ is the modulus of grating vector $\mathit{k}$, n is the refractive index of the dielectric media and a is the cylinder radius. The scalar approximation provides significant insight into the scattering problem by simplifying the analysis. However, it should be noted that the solutions obtained using this scalar approximation are analogous to those derived from the full vector treatment in the case of TE (Transverse Electric) polarisation. In this context, TE polarisation implies that the electric field is perpendicular to the plane of incidence, making the scalar approximation of the field a valid representation of the vector field under these specific conditions. As described by [16].

3. Numerical Results

First, analytical solution obtained by the Fourier series development of the angular harmonics are compared with the scalar formalism of the Green’s function. In order to achieve this, consider the following coordinate scheme presented in Figure 1. Where a plane wave polarized perpendicular impinges on the cylinder in the $xz$-plane with a tilt angle $\zeta$, and the intensity pattern of the scattered light is measured on a screen. On the other hand, the screen where the intensity pattern is collected can be moved along the axial axis y and the intensity can be gathered at different distances from the cylinder.

In the following simulations we can see the axial evolution of the normalized intensity of the scattered light in Figure 2a, moving towards the screen respect the cylinder, and the pattern of the scattered light on the screen at a distance of $y=10$$\mathrm{mm}$, Figure 2b. We can see that the proposed numerical method with the Green’s function approximates and behaves in an analogous way to the analytical solution proposed by the Equation (17). Moreover, in Figure 2a we observe that for values close to $y=0$ the behaviour of the axial intensity curve presents some abrupt changes in the intensity, which is in agreement with the scalar theory of diffraction being in the Fresnel regime close to the source.

Modifying the dimensions of the dielectric medium, in this case reducing its dimensions by half ($a=20$$\mu m$), we can observe how the pattern of the scattered beam on the screen widens with respect to the intensity pattern for a cylinder twice the size, Figure 3b respect Figure 2b. This result shows that the light scattered by the dielectric medium actually has its origin in the diffraction that the cylinder itself performs when incident on light. On the other hand, the axial behaviour of the intensity, Figure 3a of the scattered beam shows less variation as it approaches the source and decays faster than for a larger scattering medium in comparison with Figure 2a. These two comments are consistent with the scalar theory of light diffraction. On the other hand, the solution obtained by the scalar Green’s method fits better to the analytical solution (17).

Once it has been demonstrated that the numerical method fits the analytical prediction, we decided to try different geometries and to analyze the scattered field in the plane. Thus, the following results represent the light scattering of different dielectric surfaces. For this purpose, 2D scalar Green’s function given by the Equation (9) is used. In Figure 4a we consider a cylinder with an elliptical base with a semi-major axis $a=30$$\mu m$ and semi-minor axis $b=10$$\mu m$ with refractive index $n=1.5$. A planar wave front with a wavelength of $\lambda =633$$\mathrm{nm}$ propagating in air impinges normally on the surface, $\zeta =\pi /2$. We can observe the pattern of the scattered field $E^{\left(3\right)}$ once the electromagnetic wave interacts with the dielectric cylinder. This presents two clear scattered wave fronts which refers to the reflection of the incident wave $E^{\left(1\right)}$. On the other hand, in Figure 4b the intensity of the refracted field inside the dielectric cylinder is represented.

Trying other geometries we opted for a prism with a rectangular base $20\times 30$$\mu m^2$. We have in the following Figure 5a the scattered field by this new dielectric medium for a refractive index $n=1.5$ and the internal field for different refractive indices, $n=1.5$ in Figure 5b and $n=1.05$ in Figure 5c. We can observe how the homogeneity of the transmitted electromagnetic wave has a certain dependence on the difference between refractive indices of the dielectric medium with respect to the external medium, in this case air $n=1$. These results make sense because the intensity of the field transmitted by the dielectric medium depends on the refractive index due to the Fresnel equations, so the greater the relative difference in refractive indices, the less electric field will be transmitted and the less scattered field at the output of the dielectric medium. Furthermore, we can see how the areas of highest intensity are at the interface between refracting surfaces, and this intensity decreases as the electric field propagates through the interior of the prism as it occurs in the dielectric cylinder. Good agreement between scalar Green method and vector form [19,20] is observed.

4. Conclusions

The simplicity of the equations derived from the Green’s functions combined with the proposed numerical method makes this approximate method a good approach to deal with electromagnetic scattering problems. Moreover, this numerical method has been compared with other analytical solutions proposed in the literature and it has been found to fit properly. In addition, phenomena such as light diffraction have been verified by modifying the dimensions of the selected dielectric element, in this work a homogeneous cylinder. Another aspect shown is the transmission of light through a dielectric medium. To sum up, we can say that Green’s method makes it possible to recreate the phenomena of both diffraction and refraction of light. However, as it is a scalar method, information about the degree of polarization of the scattered wave is missed. Thus, this aspect that could be carried out with a vectorial numerical method based on a tensor formalism with Green’s functions.