Ch. Hafner, Computational Optics Group, IFH, ETH, 8092 Zurich, Switzerland

Back to MaX-1, MaX-1 Examples Overview

**Download ZIP file containing the
MaX-1 projects described here**

Photonic crystals are artificial, periodic structures consisting of at least 2 different materials with high contrast of the refractive index. As in natural crystals, one can describe the behavior of photonic crystals by the electromagnetic modes of the structure. For obtaining the modes, an eigenvalue problem is to be solved, i.e., the situation resembles the computation of waveguide modes and of resonators. Because of the periodicity of the photonic crystals one also has the same situation as for the grating computation. As a consequence, one must select both eigenvalue problem and periodic symmetry in the Project dialog of MaX-1 when one solves such problems with MaX-1. In the following, the main steps are explained. In order to keep the model as simple as possible, we consider the following model.

Let us consider the 2D photonic crystal consisting of circular dielectric rods with diameters 0.3 microns and relative permittivity 8.41, lattice constants dx=dy=a=1 micron in both x and y directions. For this crystal, the E waves exhibit band gaps that shall be analyzed.

First, set the project data as for periodic structures, i.e.,
select time-harmonic fields, specify the E wave type and specify the lattice
constants in the **Project** dialog. Remember that you should not set any
symmetry numbers for periodic problems. Now, specify the eigenvalue data. Note
that - for computing the Bloch modes - you search for frequencies. Since there
are no losses in the model, the search range is real. With increasing
frequencies, the number of modes is also increased. In order to obtain not too
many modes, select the upper bound of the frequency range not too high. A
reasonable value might be in the order of 3E14 Hz, which corresponds to 1
micron wavelength, i.e., the lattice constant.

Now, specify the material properties in the Domain dialog and specify three boundaries. Let the first boundary be a circle with center at x=y=0.5 microns (in the center of the original cell of the periodic structure) and radius 0.15 microns. Impose standard boundary conditions on this boundary. The second and third boundaries are periodic boundaries of the original cell. You best use straight lines of 1 micron length, located at x=1 micron and y=1 micron respecively. Impose the appropriate periodic boundary conditions on these boundaries and make sure that the domain numbers on both sides of these boundaries are equal, i.e., model the boundaries exactly as for any structure that is periodic in x and y directions.

Specify a set of expansions along the boundaries and along the two boundaries of the original cell that were not explicitly defined - exactly as for any structure that is periodic in x and y directions. Since you have a circular rod at x=y=0.5 microns you best use a 2D Bessel expansion at x=y=0.5 microns for modeling the field inside the rod and a multipole with the same maximum order at the same position for modeling the field outside the rod. Now, add multipole expansions outside the original cell for modeling the field inside the original cell (outside the rod). These expansions take the interactions with the rods outside the original cell into account. You best place 4 multipole expansions at 1) x=1.5 microns , y=0.5 microns, 2) x=-0.5 microns , y=0.5 microns, 3) x=0.5 microns , y=1.5 microns, and 4) x=0.5 microns , y=-0.5 microns. Use identical orders for all expansions and start with moderate maximum orders. Once you have obtained your first results, check the errors along the boundaries and increase the maximum orders if necessary.

In eigenvalue problems obtained from resonators you search
for the resonance frequencies of all modes. that only depend on the geometry and
material properties of the model. For photonic crystals, the situation is more
complicated, because of the periodicity of the structure. Note that you must
define not only the lattice constants, but also the** CX **and **CY**
values in the **Project** dialog for structures that are periodic in X and Y
directions. These values depend on the lattice constants and the x and x
components of the wave vector. For each pair (kx, ky), i.e., for each pair (CX,CY)=(kx*dx,ky*dy)=(kx*a,ky*a),
you obtain a new eigenvalue problem. Because of the high symmetry of the
problem, it is sufficient to study the so-called Bloch modes in the first
Brioullin zone only. The first Brillouin zone of our square lattice is a
triangle in the k (Wavevector) space with the three corners Gamma (kx=ky=0), X (kx=Pi/a,
ky=0, where a=dx=1 micron), and M (kx=ky=Pi/a). Usually, one only considers the
eigenvalues along the border of the first Brillouin zone, i.e., from Gamma to X,
form X to M, and from M to Gamma. According to that, it is best to split a
project into projects, one for each side of the triangle mentioned above. In
order to trace the eigenvalues along one side of the triangle, you can take
advantage of the eigenvalue estimation technique (check the **PET** box in
the **Eigenvalue search** group of the **MMP** dialog. Furthermore,
setting the CX and CY values requires some extra computation - especially for
non-square lattices. In order to overcome this problem, you may use the movie
directives "SET PERiod Phase 0 0 0", "SET PERiod Phase 180 0
0", and "SET PERiod Phase 180 180 0" for setting the Gamma, X,
and M point respectively. In this directive, Phase means the the value of kx*a,
ky*a, (kz is not used because there is no periodicity in Z direction) in
degrees. Along the kx axis, from Gamma to X, kx*a varies from 0 to Pi, i.e., the
"phase" varies from 0 to 180. When you want to obtain the Bloch modes
along this line in 180 steps, 1 degree per step, you first set the Gamma point
with the directive "SET PERiod Phase 0 0 0" and increase the
"phase" using "INCrease PERiod Phase 1 0 0" 180 times in a
loop.

Since we are solving an eigenvalue problem, the last
expansion plays the role of the excitation. When you know something about the
mode you want to compute, you should place the expansion that is dominant at the
last position. When you want to compute the band diagram, you would like to
obtain all modes rather than a special one. Therefore, it is reasonable to
introduce a "fictitious excitation" as the last expansion. The
fictitious excitation may be an arbitrary monopole that should be located at a
place where it can excite all modes. For example, you can put it at a random
position inside the rod. Note that the domain number of the fictitious
excitation is the same as the domain number of the interior of the rod.
Therefore, when you plot the field in the rod, you will obtain a singularity at
the origin of the fictitious excitation. Since the field of the mode you compute
does not contain the fictitious excitation, you best exclude the fictitious
excitation by setting the value -6 (you have defined 6 expansions for modeling
the field of the mode, the fictitious excitation is expansion number 7) in the **use
expansion** box of the **Field formula** dialog.

When the frequency is close to the resonance frequency of a
mode, the field will become huge - except at special points, where the field of
the mode is zero. When you evaluate the field at an arbitrary point in the
original cell, the field will be a function of the frequency and this function
will exhibit peaks that represent the resonances. In order to obtain a scalar
function, you can either evaluate only a component of the field or a derived,
scalar quantity. The most natural quantity is the total energy density of the
electromagnetic field. When you select **Rect. Amplitude definition** in the **MMP**
dialog and specify a small **Rectangular area** with **nx**=1 and **ny**=1
in the **Integral** dialog, MaX-1 will evaluate the field in a single point
and use it as the Amplitude for the eigenvalue search. Since the standard
eigenvalue search function is defined as f0=Residual/Amplitude, f0 should have a
minimum when Amplitude has a maximum. This is fine because MaX-1 searches for
minims. Unfortunately, the Residual also exhibits sharp peaks near the resonance
frequencies. Therefore, you better replace the search function by f1=1/Amplitude
by not checking the **Use residual **box in the **MMP** dialog. At first
sight, this function seems to be very nice. It exhibits nice valleys at the
positions of the resonances.

Animation: Eigenvalue search function f1 from the point Gamma
to the point X of the first Brillouin zone. Generated with the example projects **SQR000**,.
For the corresponding solutions from X to M and M to Gamma, run the projects **SQR001**,
and **SQR002**.

When you have a closer look at the animation above, you see some strange needles in some of the minims. In order to understand this, you should have a closer look. Before you do this, you should modify the fictitious excitation in order to make sure that all eigenvalues are detected. Try the following: set the excitation at the center of the rod. Because of the symmetry, it is very likely that some modes cannot be excited by such an excitation. If this happens, some of the minims should be suppressed. Try also to use several fictitous monopole excitations by embedding them in a connection.

Animation: Eigenvalue search function f1 from the point Gamma
to the point X of the first Brillouin zone. Red: single monopole excitation at a
non-symmetric point, Green: monopole excitation at the center of the rod, Black:
4 monopoles in a connection used as fictitious excitation. Generated with
the example projects **SQR100**,. For the corresponding solutions from X to M
and M to Gamma, run the projects **SQR101**, and **SQR102**. Note: the
four monopole excitations of the Black curve are on a regular grid, parallel to
the x and y axes, with origins in the center of the rod. Such a regular
distribution of excitations causes additional problems (unphysical minims at
high frequencies).

When you zoom in and let MaX-1 draw the eigenvalue search function f1 very close to a resonance, you can see that there are two (downward) peaks rather than a single one.

This "twin peak" behavior is even more pronounced
when you watch f0 rather than f1. From this, you can guess that it is related to
the behavior of the residual. When you want to inspect the function f2=Residual,
you select **1** type of the **Amplitude definition** in the **MMP**
dialog and check the **Use residual **box.

As you can see, f2 has a nice single peak. This peak
specifies the eigenvalue very precisely. In order to obtain a minimum rather
than a maximum, you may invert the peak by selecting -1 in the **Exponent**
box of the **MMP** dialog. You will obtain f3=1/f2=1/Residual. Unfortunately,
the minims of f3 are so sharp that a very fine grid for the eigenvalue rough
search is required in order to make sure that all existing eigenvalues are
detected. This leads to undesired, high computation times of the rough search.
Therefore, you might use different search function definition for the rough and
fine search routines. Note that this is only important when you want to compute
the eigenvalues with a high precision - the twin-peak phenomenon of the
functions f0 and f1 only disturbs the fine search when a high accuracy is
desired. An alternative is to define a more complicated eigenvalue search
function using the Amplitude and Residual data. The following definition seems
to be useful: f4=((1/Amplitud)**(1/Exponent))/(Residual**Exponent). You obtain
this eigenvalue search function when you you select **Rect.** type of the **Amplitude
definition** in the **MMP** dialog , leave the **Use residual **box
unchecked and specify the **Exponent** value in the corresponding box. When
Exponent is too small, the twin peaks will be there, whereas a big Exponent can
cause numeric underflow problems. For our example, Exponent 5 is a good
compromise. This leads to the following search function:

The example projects **SQR010**, **SQR011**, and **SQR012**
contain all data for obtaining the three figures shown above.

When you want to draw the dependence of the Bloch modes along
the border of the first Brillouin zone, you might start the rough eigenvalue
search at the Gamma point and use the eigenvalue estimation technique to trace
all eigenvalues from there to the X point, to the M point and back to the Gamma
point. From the animations above, you can see that some of the modes are
degenerated at the Gamma (and X and M) point of the first Brillouin zone. Thus,
you MaX-1 cannot distinguish all modes when you start the rough search in the
corners of the first Brillouin zone. Moreover, the eigenvalue estimation might
switch from one of the degenerated modes to the other one when you run through
such a point. Therefore, a more complicated procedure is required: Start between
the Gamma and X point and trace once from the start point to the X point and
once to the Gamma point. For doing this, you can run the rough search only once.
The project **SQR200** with the directives SQR200.DIR (search from middle
point to X point) and SQR210.DIR (search from middle point to Gamma point - note
that SQR210.DIR is loaded automatically from SQR200.DIR) illustrates this.

Figure: Band structure of the first 5 Bloch modes from the Gamma point (at x=0) to the X point (at x=180), from X to M and from M back to Gamma..

The projects **SQR220** and **SQR240** compute the
modes from the X to the M point and from the M to the Gamma point.