SQRxxx.PRO: Photonic Crystal with Square Lattice

Band diagram of a perfect photonic crystal, Photonic band gap, 2D MMP, Advanced eigenvalue search, Periodic structure

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.

Circular dielectric rods on a square lattice

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 modeling steps

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.

Eigenvalue definition

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).

Eigenvalue search functions

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.

Tracing eigenvalues

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.

Back to MaX-1, MaX-1 Examples Overview