Multiple aperture vector diffraction

Introduction

Electromagnetic diffraction is a phenomenon that occurs when an electromagnetic (EM) wave passes trough an obstacle or aperture. For EM diffraction the electric or magnetic field are described by a complex function. In the case of two or more apertures the spatial distribution of the EM field after the aperture, called the diffraction pattern, is characterized by fringes that vary with the distance between apertures due to the spatial phase difference. This feature is greatly exploited in metrology to measure distance of micro and nano objects.

In isotropic and non-dispersive media, diffraction of EM waves can be modeled by the vector valued source-free Helmholtz partial differential equation [Jackson, 1999]:

Here is the wave number in the medium and is related to the wavelength by . The wavelength in a medium with refractive index is related to free space wavelength by . In some cases only one component of the field is required to describe the full behavior of the field, in that case the field is said to be scalar. The scalar case was studied in the "Single Aperture Scalar Diffraction" tutorial. This tutorial provides a brief introduction to multiple slits using a scalar approach, followed by a focus on the vector properties of EM waves.

Polarization is a property of transverse vector fields that specifies the geometric orientation of the oscillations. These property is exploited in a technique called polarimetry to characterize optically active samples.

As a prerequisite for the content presented in here, it is recommended to work through the material presented in "Single Aperture Scalar Diffraction" as the basics of electromagnetic diffraction and perfectly matched layer are described in there.

Many of the animations of the simulation results shown in this notebook are generated with a call to Rasterize. This is to reduce the disk space this notebook requires. The downside is that the visual quality of the animations will not be as crisp as without it. To obtain high quality graphics remove or comment out the call to Rasterize.

To get high fidelity visualizations comment out the rasterization process.

The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.

Load the finite element package.

Scalar diffraction by multiple apertures

Before modeling the vector Helmholtz equation (1), lets solve the simpler case of scalar diffraction by multiple apertures. This will help to introduce the boundary conditions and element mesh functionalities in a more understandable way.

For this first example, consider a laser radiating an electromagnetic (EM) field propagating in the direction and polarized in the direction. Under those conditions the vector electric field can be consider as a scalar field containing energy only in the component [Born & Wolf, 1999].

This application example considers a Helium-Neon laser with a free space wavelength illuminating an opaque screen with 4 evenly spaced apertures of width of and separation . The apertures are slits assumed to be infinitely large in the direction, so there will be no disturbances to the diffracted wave in the direction, hence a system in 2D can be modeled. The laser and apertures itself will be modelled by a radiation boundary condition on the boundary. The opaque screen is made of a fully absorbing material where light cannot be transmitted or reflected. The boundary indicates the beginning of the perfectly matched layer (PML) region. The PML region is an artificial boundary that absorbs incoming light to mimic an open system where light propagates indefinitely without any reflections [Jin, 1993], therefore the only source of light is the aperture as the field is in the other boundaries. The simulation domain then has dimensions and . The refractive index inside the propagation region is 1.

24.gif

The scheme above is rotated 90 degrees for better visualization, in the following the propagation axis will be the vertical axis. Also note, that the geometry is symmetric along the direction and that the directions points outwards from the screen.

Domain

To make use of PML approach, the propagation region has to be extended by an extra layer of length which will be the PML region. In this region the incoming wave will be damped out.

Create a boundary mesh with the PML region and set up the boundary element markers.
Visualize the PML regions and boundary element markers.

Boundaries 1 and 3 (dotted green) describe the outer and inner boundaries of the PML region. Boundaries with marker 4 (orange) are needed to setup the PML parameters. The boundary 2 (red) describes the opaque screen and apertures. To distinguish between the opaque screen and apertures, a function will be used inside the Neumann boundary condition.

From the "Single Aperture Scalar Diffraction" tutorial, it is observed that the intensity profiles in regions far from the -axis are negligible. A MeshRefinementFunction can be utilized to create a finer mesh at specified coordinates, and "RegionMarker" can be used to divide the PML region into four subregions. In 2D, a MeshRefinementFunction operates based on the element's area. The area of a triangle in the mesh can be approximated using the formula for the area of an equilateral triangle, which is , where represents the edge length of the triangle.

It is worth mentioning that is a recommended initial discretization constraint to properly resolve the EM wave [Jin, 1993]. In some systems a smaller constants has to be used to provide a precise result.

Set up model parameters.

Now, create the mesh with a maximum element area depending on the position. In the region , the recommended edge length is used, as it is the region of interest. Outside that region, the field is practically 0. Therefore, a larger edge length is used to reduce the number of unknowns for the FEM method.

Create an element mesh with variable max cell measure and create 4 different PML subregions.
Visualize the domain and PML regions.

The use of Rasterize above helps to reduce the size of the graphics. For a better visual appearance the Rasterize can be removed.

Parameters setup

The basics of PML are described in the acoustic time and frequency domain tutorials. Roughly speaking, a PML region mimics an absorbing material with an absorbing coefficient . While is set inside the propagation region , an increasing function with maximum value is set inside de PML region [Jin, 1993] and there dampens the wave out.

Set up the Helmholtz operator with an artificial electric susceptibility tensor that correspond to the PML parameters .
Define the absorbing coefficient for all regions.
Visualize the absorbing coefficient on the element mesh.

The use of Rasterize above helps to reduce the size of the graphics. For a better visual appearance the Rasterize can be removed.

In the plots, it is evident that the absorbing coefficient behaves as a wall in each direction outside the region , while it does not alter the electric field inside the region .

Boundary conditions

The incoming wave illuminating the boundary containing the opaque screen and aperture can be modelled with a radiation boundary condition of the form [Jin, 1993]:

with the boundary normal, the incoming EM wave propagation direction and the EM wave amplitude. The boundary condition is a sum of two terms, a light incidence term and an absorbing boundary condition term .

The second term is needed to ensure that the boundary does not reflect EM waves and only the incident light flows outward from the boundary. In this case, both, the boundary normal and the wave propagation direction, are in the z-direction, therefore the inner product is 1. will be used to distinguish between the screen and apertures. From equation (2), it is evident that setting to 0 results in only the absorbing term, which represents an opaque screen. However, setting to represents an aperture with an initial amplitude of .

Set up a function for multiple apertures that returns in the apertures and 0 at the opaque screen.

where is the number of apertures, is the width of each aperture and distance between the centers of any two apertures.

Visualize the apertures of width , spacing and initial amplitude .

Now, set up the radiation boundary condition with the following parameters , , and .

Set up the radiation boundary condition in which will contain both the apertures and screen conditions.

Model Evaluation

Set up the model and replace the PML parameters.
Solve the PDE with NDSolve.

Visualization

Plot the numerical result.

See this note about improving the visual quality of the animation.

Visualize the complete solution as a DensityPlot.

Model verification

Next, the solution is validated against an analytical solution. To verify the model, the integral diffraction formula deduced by Green's method [Ishimarou, 1991] can be used. In 2D, the analytical solution is as follows:

where is the euclidian distance between any point at and a point inside the region . As the aperture is a line, the normal is a vector in the direction. For a detailed deduction follow [Ishimarou, 1991].

Create a helper function to compute the analytical solution.
Visualize the error between the numerical and analytical solution.

A more detailed discussion about the analytical solution is given in "Single Aperture Scalar Diffraction" tutorial.

Vector field diffraction

In the above section, the diffraction of scalar fields through an object was described. Scalar fields have a single value at each point in space. On the other hand, vector field have an amplitude and a direction, therefore the field can be written as with the amplitude and the orientation or polarization angle.

In a medium with a diagonal relative permittivity tensor the components of the electric field are independent and can be solved individually [Born & Wolf, 1999].

In the following, the diffraction of a vector field will be modeled by two apertures with a linear polarizer in each one. A linear polarizer is an optical filter instrument that only allows EM waves with a specific polarization angle to pass through. The polarizer angle is called the transmission axis.

For example, when an EM wave with a polarization angle pass through a linear polarizer with transmission axis , the EM wave polarization angle after the polarizer will be , but with a reduced intensity. The intensity will be and the remaining intensity will be absorbed by the polarizer. This is called the Malus' law. One important note, is that if and are orthogonal () the final intensity will be zero.

A linear polarizer with transmission axis is characterized by a transmission matrix [Born & Wolf, 1999].

Assume that the linear polarizer located at the upper aperture is fixed at , i.e. only the x-component of the incident field passes through the aperture. The other polarizer located at the lower aperture has a variable transmission axis .

96.gif

The scheme above is rotated 90 degrees for better visualization, in the following the propagation axis will be the vertical axis. Also note, that the directions points outwards from the screen.

Parameter setup

Set up the vector Helmholtz operator for the field and the diagonal anisotropic (orthotropic) permittivity for the PML region. In this case an additional rank 3 tensor is needed to set up both PDEs simultaneously, see the DiffusionPDETerm reference page for more information.

Set up the material parameters.
Set up the dependant variables and permittivity .
Set up the vector Helmholtz equation with the PML parameters.

Domain definition

The same domain, ElementMesh and PML parameters are used as for the 2D PML scalar example presented above.

Reuse the previously defined mesh.
Reuse the absorbing coefficient for all regions.

Boundary conditions

Before establishing the boundary conditions for the apertures with the linear polarizers, an auxiliary function for the polarizer transmission matrix needs to be created.

Create an auxiliary function for the transmission matrix of the polarizers.

Now, the same approach as in the scalar example can be used. Set up a function for the amplitude for each aperture. In this case the amplitude is not a scalar, it is a 2D vector. As linear polarizer is placed in each aperture, one with transmission axis at 0° and the other has an arbitrary transmission axis , the transmission matrix needs to be multiplied to the incoming vector field amplitude.

Set up the function for the two apertures with the polarizers considering a incident field with amplitude .

is the width and is the spacing between the two apertures center.

Visualize the behavior of the incident field after the polarizers.

See this note about improving the visual quality of the animation.

Set up the radiation boundary conditions for the two components of the field at an arbitrary polarization angle .

Model Evaluation

For this example, the same mesh element and PML regions employed in the 2D PML example will be used.

Set up the PDE and PML parameters.
Use ParametricNDSolveValue to evaluate the solution for different transmission axis of the variable polarizer.

Now, do a sweep for different values of between and .

Visualize the change of each field component for different angles.

See this note about improving the visual quality of the animation.

When , both apertures have the same polarization and a result is obtained equal to the scalar case studied previously. When starts increasing, the component of the EM wave going out of the upper aperture will interfere with the component of the EM wave going out of the lower aperture. The amount of energy in the component varies with . In the other hand, The EM wave going out the upper aperture does not contain energy in the component, therefore the component of the EM wave going out of the lower aperture does not interfere with anything, producing a single slit pattern which magnitude varies with . When the fields diffracted at the upper and lower aperture are orthogonal, therefore showing two single aperture patterns.

The intensity of the field will be proportional to the sum of the squares of all the components i.e Ex2+Ey2.

Nomenclature

SymbolDescriptionUnit
EElectric field[V/m]
λnWavelength[m]
knWavenumber[1/m]
μrRelativepermeabilityN/A
ϵrRelativePermittivityN/A
LPropagationregionlength[μm]
HPropagationregionwidth[μm]
dAperturesize[m]
ϕAperturesize[m]
ΓscreenOpaquescreenboundaryN/A
ΓapertureApertureboundaryN/A
ΓABCAbsorbingBoundaryN/A
ΓinScreenandapertureboundaryN/A
ΩComputationaldomainN/A

References

1.  Born, M. and Wolf, E.(1999). Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (7th ed.). Cambridge: Cambridge University Press.

2.  Jackson, John D. Classical Electrodynamics. New York :Wiley, 1999.

3.  Jin, Jian-Ming. The finite element method in electromagnetics. John Wiley & Sons, 2015.

4.  A. Ishimarou, Electromagnetic Wave Propagation, Radiation, and Scattering, Prentice Hall, 1991.