COMSOL 6.4 - Hexagonal Grating

A hexagonal grating is an infinite structure that is periodic with hexagonal (or rhomboid unit) cells. Figure 1 shows the hexagonal domain used for this model. The reflecting perfectly conducting surface consists of regularly spaced protruding hemispheres.

Figure 1: The hexagonal domain, used for computing the diffraction from the hexagonal grating.

As shown in Figure 2, for a hexagonal cell of side length a, the corresponding unit cell is a rhomboid with side length

. In Figure 2, the side vectors for the hexagonal cell starts from the point P and are denoted a1 and a2. The angle between a1 and a2 is 120 degrees. Similarly, for the rhomboid unit cell, the primitive vectors are denoted u1 and u2 and starts from the hexagon center point Q. The angle between the two primitive vectors is also 120 degrees.

Figure 2: Schematic showing the hexagonal cells with side length a and side vectors a1 and a2. The primitive cells are defined by the primitive vectors u1 and u2.

If the incident plane wave have a wave vector defined by

(1)

where k|| is the wave vector component parallel to the periodic boundary and k⊥ is the component orthogonal to the periodic boundary, the in-plane wave vector component for diffraction order mn is given by

(2)

where the reciprocal lattice vectors G1 and G2 are defined from the primitive vectors u1 and u2 as

(3)

and

(4)

where n is the normal vector (length 1) to the periodic boundary.

Since the out-of-plane wave vector component for mode mn is defined by

(5)

it is clear that for propagating modes, where the out-of-plane wave vector component above must be real, the length of the in-plane wave vector component must be smaller than the material wave number k. Figure 3 shows that it is only the modes inside the circle with radius k that will be propagating. In the example shown in Figure 3, there are five modes that will be propagating, in this case the modes m = n = 0 (the reflected wave), m = -1, n = 0, m = 0, n = -1, m = -1, n = -1, and m = -2, n = -1. All other modes will be evanescent and damped out.

Figure 3: The reciprocal lattice, showing the reciprocal lattice vectors G1 and G2, the in-plane wave vector component k||, and the circle with radius k (the material wave number) enclosing the propagating mode points (larger dots. The dotted hexagon indicates that also the reciprocal lattice is a hexagonal point lattice. The dashed rhomboid indicates the unit cell spanned by the reciprocal lattice vectors.

Model Definition

In this model, the unit cell is small compared to the wavelength, so there will only be two modes that are propagating, the modes m = 0, n = -1 and m = -1, n = -1. For wavelengths longer than approximately 1.01 μm (the critical wavelength), the mode m = 0, n = -1 will be evanescent.

First a wavelength sweep will be made for an incident field having the polarization perpendicular to the plane of incidence (spanned by the wave vector for the incident wave and the normal to the periodic boundary) (so called s-polarization). Thereafter another wavelength sweep is made, but now with the polarization in the plane of incidence (p-polarization).

Results and Discussion

Figure 4 shows the electric field norm and the propagation directions for the incident, the reflected and the diffracted waves. Notice that the diffracted waves come in pairs (both have the same mode numbers), one wave having the polarization in the plane-of-diffraction and the other wave have orthogonal polarization to the plane-of-diffraction. The plane-of-diffraction is spanned by the wave vector for the diffracted wave and the normal to the periodic boundary. The wavelength is close the critical wavelength for the m = 0, n = -1 mode. This is evident from the plot, as the wave vector for that mode (the yellow arrows) is almost parallel to the periodic boundary.

Figure 4: The electric field norm and the propagation directions for the incident wave (red arrows), the reflected wave (blue arrows) and the two diffraction orders (green and yellow arrows). The wavelength is 1.01 μm, which is close to the critical wavelength for the mode m = 0, n = -1, and the polarization of the incident wave is perpendicular to the plane of incidence.

Figure 5 shows the reflectance (for mode m = n = 0) and the diffraction efficiencies for the diffracted waves. Notice that both the reflectance and the diffraction efficiency for the in-plane-polarized m = -1, n = -1 mode show resonances (peaks or dips) close to the critical wavelength for the m = 0, n = -1 modes.

Figure 5: Diffraction efficiencies for the reflected wave and the diffracted waves. The polarization of the incident wave is perpendicular to the plane of incidence.

Figure 6 shows polarization plots for the same case as Figure 4 and Figure 5. It is clear that even though the polarization for the incident wave (the dotted horizontal line) is linear, the polarization for the reflected diffraction orders are elliptical. The size of the polarization ellipses indicates the diffraction efficiency. So, the zeroth diffraction order is the strongest and the m = -1, n = -1 order is the weakest one. The blue solid line indicates the area in the m-n-plane where there can be diffraction orders for propagating waves. As there does not exist any other mode numbers in this area, this is a verification that the automatically added diffraction orders are the correct ones to include in the simulation.

Figure 6: Polarization plot showing the polarization ellipses for the three diffraction orders that are not evanescent at the first wavelength, 0.9 µm. The linearly polarized input field is indicated by the dotted line. The solid blue line encircles the area in which there can be diffraction orders representing propagating plane waves.

Figure 7 shows a similar plot as Figure 4, but here the polarization of the incident wave is parallel to the plane of incidence.

Figure 7: Similar plot as in Figure 4, but here the polarization of the incident wave is parallel to the plane of incidence.

Figure 8 shows that for p-polarization both the reflected wave and the two m = -1, n = -1 modes show resonances close to the critical wavelength for the m = 0, n = -1 mode.

Figure 8: Similar plot as in Figure 5, but here the polarization of the incident wave is parallel to the plane of incidence.

Notes About the COMSOL Implementation

To define the periodic unit cell, a Periodic Structure node is added. This node automatically adds and configures its subnodes — the periodic port and the three Floquet periodic conditions.

RF_Module/Tutorials/hexagonal_grating

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Global Definitions

First add some parameters that defines the geometry and the incident electric field.

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
wl0	1[um]	1E-6 m	Center wavelength
lda0	wl0	1E-6 m	Wavelength
f0	c_const/lda0	2.9979E14 1/s	Frequency
a	wl0/2	5E-7 m	Hexagon side length
h0	3*wl0	3E-6 m	Air height
a1	a/2	2.5E-7 m	Sphere radius
theta	pi/3	1.0472	Elevation angle
phi	12[deg]	0.20944 rad	Azimuth angle

Notice that the azimuth angle phi above measures the angle for the wave vector of the incident wave from the x-axis.

Here, c_const is a predefined COMSOL constant for the speed of light in vacuum.

Study 1

Step 1: Frequency Domain

In the window, under click .

In the window for , locate the section.

In the text field, type f0.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

The geometry consists of an extruded hexagon, with a semisphere removed from it at the bottom.

Work Plane 1 (wp1)

In the toolbar, click

In the window for , click

Work Plane 1 (wp1) > Plane Geometry

In the window, click .

Work Plane 1 (wp1) > Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


xw (µm)	yw (µm)
a	0
a/2	sqrt(3)/2*a
-a/2	sqrt(3)/2*a
-a	0
-a/2	-sqrt(3)/2*a
a/2	-sqrt(3)/2*a

Extrude 1 (ext1)

In the window, right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Distances (µm)
h0

Sphere 1 (sph1)

In the toolbar, click

In the window for , locate the section.

In the text field, type a1.

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the object only.

Click

Click the

button in the toolbar.

Click the

button in the toolbar.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Electromagnetic Waves, Frequency Domain (emw)

In the window for , locate the section.

From the list, choose .

Periodic Structure 1

In the toolbar, click

and choose .

In the window for , locate the section.

Clear the checkbox, as the bottom boundary should use the default node.

Locate the section. In the α1 text field, type theta.

In the α2 text field, type phi, as this angle is measured from the subnode of the .

In the window, click .

Locate the section. From the list, choose , as no angle sweep will be done in this example model.

Click .

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the text field, type wl0/6.

Click

Study 1

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

From the list in the column, choose .

Click

In the dialog, type 0.9[um] in the text field.

In the text field, type 0.01[um].

In the text field, type 1.1[um].

Click .

In the window for , locate the section.

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
lda0 (Wavelength)	range(0.9[um],0.01[um],1.1[um])	um

In the toolbar, click

Results

Add Arrow Surface plots showing the wave vector directions for the incident field, the reflected field and the diffracted fields. Notice that the diffracted fields come in pairs, where each pair have the same wave vector. Thus, only two wave vectors for the diffraction orders need to be added in this case.

Arrow Surface 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > > .

Locate the section.

Select the checkbox. In the associated text field, type Incident wave (red).

Locate the section. In the text field, type 12.

Locate the section.

Select the checkbox. In the associated text field, type 3e-8.

Selection 1

Right-click and choose .

Select Boundary 4 only, to select only the port boundary.

Arrow Surface 2

In the window, under > right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > > .

Locate the section. From the list, choose .

Locate the section. In the text field, type Reflected wave (blue).

Arrow Surface 3

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > > .

Locate the section. From the list, choose .

Locate the section. In the text field, type Mode m = -1, n = -1 (green).

Arrow Surface 4

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose > > > > .

Locate the section. From the list, choose .

Locate the section. In the text field, type Mode m = 0, n = -1 (yellow).

Electric Field (emw)

Select the wavelength closest to the critical wavelength for the mode m = 0, n = -1.

In the window, click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Click the

button in the toolbar.

Your plot should look the same as Figure 4.

Global 1

In the window, expand the > node, then click .

In the window for , locate the section.

In the table, enter the following settings:


Expression	Unit	Description
abs(emw.S11)^2	1	Reflectance
abs(emw.S31)^2	1	m = -1, n = -1, in-plane
abs(emw.S41)^2	1	m = -1, n = -1, out-of-plane
abs(emw.S51)^2	1	m = 0, n = -1, in-plane
abs(emw.S61)^2	1	m = 0, n = -1, out-of-plane
abs(emw.S71)^2	1	m = 0, n = 0, orthogonal