1. Introduction

During the past several years research activities in the field of photonic crystals (PCs) have rapidly grown, primarily driven by the eminent advances in microstructuring and modeling. Although technologically challenging, 3D PCs are the most intriguing manifestation of the PC concept. Examples are the opal [1, 2], the inverse opal [1–3] or the woodpile structure [4–7]. From bandstructure diagrams [1,2,4,5] it is well-known that both latter structures may exhibit complete photonic bandgaps for a certain minimal index contrast. However, bandstructure diagrams provide only information about maxima and minima of the bands as well as derived quantities, e. g., the group velocity, in certain directions. This is because the bands are only computed and visualized along 1D paths in k space.

Recently the focus of research has appreciably shifted from looking at bandgap effects to the exploration of the engineerable propagation characteristics of Bloch waves in PCs [8–13]. In order to determine the propagation characteristic of Bloch waves it is indispensable to have access to their dispersion properties in the 3D k space. For example, based on this knowledge angular resolved diffraction of finite beams in a 3D PC can be calculated, because the energy velocity and the group velocity (GV) v _g = ∇_k ω are equivalent in periodic media. In films with two-dimensional PCs such isofrequency curves have been used to design custom-made PCs exhibiting superprism [8, 14, 15], self collimation [8, 10 , 12, 13, 16], and negative refraction effects [8, 17, 18].

These film PCs have a 2D Brillouin zone (BZ) with 1D curves as isofrequency surfaces (IFSs). Recently, a kind of self-collimation in a 3D simple cubic PC has been studied by means of sections through the 3D BZ [19] leading again to 1D isofrequency curves. Using this approximate approach light propagation can only be studied in the plane of these sections. Moreover, the full IFS of a 3D body centered cubic (bcc) crystal was used to confirm a negative index behavior [20].

In this work we calculate the dispersion relations of low- and high-index woodpile-structures [4–7] and for a high-index inverse opals [2,3] in the complete 3D BZ of these 3D PCs. We show the resulting IFSs of different bands of two high-index PCs. Inspecting these 2D manifolds regions of self-collimation are identified for the three different crystals. These predictions are numerically verified by tracing the light path with finite-difference time-domain (FDTD) calculations. The method used can also be applied to search for frequency domains where negative refraction and/or diffraction appears.

2. Computational methods and visualization

The bandstructures presented here were obtained from the calculation of fully-vectorial eigen-modes of Maxwell’s equations with periodic boundary conditions by preconditioned conjugate-gradient minimization of the block Rayleigh quotient in a plane wave basis, using a freely available software package [21].

This method was likewise used to obtain the IFSs by solving the vectorial eigenvalue equation for ω as a function of k. For the calculation of the dispersion relation in the complete 3D BZ this approach is more efficient than that which calculates k_z for a given set (ω, k_x , k_y ) [22]. We restricted k to half an octant in reciprocal space, where the vectors in the plane of the rod layers were restricted to a triangle and no restriction was applied in the perpendicular direction. Then, the complete 3D dispersion relation can be constructed by taking advantage of the symmetry of the crystal.

For sufficiently accurate results 64 × 64 × 64 plane waves have to be taken into account, corresponding to a mesh of 64 points in each of the fcc lattice directions.

The dispersion relation ω(k) is visualized as IFSs in the normalized cartesian cube k/Λ ∈ [-1, 1]³ with 2Λ = 4π/a being the lattice constant of the conventional cubic cell of the bcc reciprocal lattice. Note that for better visibility all depicted IFSs are clipped by the plane located at k/Λ = (1/2,1/2,1/2) in (1,1,1)-direction.

It is known that the curvature of the IFSs and the modulus of the GV at a certain point in k space determine the diffraction properties of a finite beam [9,11]. The curvature at a point of the surface varies between a minimum and a maximum value, called the principal curvatures κ ₁ and κ ₂, which are in perpendicular directions. Thus, in order to evaluate the IFS curvature, we calculated the principal curvatures κ ₁ and κ ₂ at every point of the IFS from the implicit derivatives. To improve accuracy the first derivatives (v _g) were calculated with the Hellmann-Feynman theorem from the fields [21]. Self-collimation in two dimensions requires that both curvatures tend to zero [11]. Thus, for visualization the isolines of the square average curvature $\overline{\kappa }$ = [(${\kappa }_1^2$ + ${\kappa }_2^2$)/2]^1/2 were mapped onto the IFS.

The light propagation in time domain was simulated with 3D FDTD [23] calculations. Here Maxwell’s equations are discretized directly in space and time without further approximation and are propagated in time.

3. Propagation in woodpile photonic crystals

Although effects such as self-collimation do not require a bandgap it will prove useful to operate the crystal close to such a gap because then usually only a few bands are involved for a given frequency. Therefore as a first representative structure we choose a face centered cubic (fcc) woodpile structure with material filling fraction f = 28% and a refractive index of n = 3.4, which is typical for semiconductor PCs. In the fcc woodpile crystal the height of four layers of rectangular rods equals a = √2d, where a is the lattice constant of the conventional cubic cell and d is the in-plane rod distance. The filling fraction is given by f = w/d, where w is the rod width. This structure is known to have a large complete photonic bandgap (gap-to-midgap ratio 19%, 0.48 < a/λ < 0.58) between the second and the third band [5] (see Fig. 1).

We restricted ourselves to the four lowest order bands, because the bandgap appears in this frequency domain (see Fig. 1), and overlaps of different bands, which would certainly deteriorate self-collimation, are minimized. By using the iterative eigenvalue solver and locating the origin at the inversion center of the woodpile structure CPU time and memory requirements have been considerably reduced.

 

Fig. 1. Bandstructure for the fcc woodpile structure with f = 28%, n = 3.4 (left) and f = 40%, n = 1.6 (middle) and for the inverse opal with n = 3.4 (right). X′, U′, K′ and W′ are the high symmetry points with the larger z components of k, obtained by exchanging k_y and k_z of X, U, K and W.

Download Full Size | PDF

The normalized frequencies Ω = ω/cΛ = a/λ (c is the velocity of light in vacuum) selected in this work are 0.44 for the first two bands, 0.69 and 0.73 for the third and fourth band, respectively. At these frequencies the IFSs exhibit regions of very small $\overline{\kappa }$, allowing for self-collimation in the respective direction of an quasi-isotropic radiation of a point source [11]. In Fig. 2 these IFSs with the isolines of smallest normalized average curvature K̅ = $\overline{\kappa }$Λ are displayed. Obviously, the IFSs reflect the crystal symmetry, i.e., ω($C_4^z$ k) = ω(k), where $C_n^{\mathbf{a}}$ denotes a counter-clockwise rotation by an angle 2π/n around the axis a. Similar symmetries do not hold for the other axes because only x and y axis are equivalent.

In the long-wavelength limit (Ω ≪ 1) the crystal acts as uniaxial homogeneous medium (z being the optical axis), where the principal refractive indices can be derived from Ω(k) and amount to n _o= 1.7 and n _e = 1.6. Extending these results to Ω = 0.44 the minimum average curvature of this homogeneous medium is K̅ = Λ/|k| = 1/0.44/1.7 = 1.34. This value will serve as a reference.

Now, inspecting the IFS for case (a) in Fig. 2 we obtain the points of smallest average curvature, which are situated in the x-y plane, namely k/Λ = (0.81,0.41,0) and its symmetric copies. Detailed calculations provide then the average curvature K̅ = 0.134 and GV v _g/c = (0.25,0.005,0), where the latter value indicates that the GV will play only a minor role in the angular-dependent diffraction because it varies not so much compared to its maximum |v _g|/c = 0.37.

Thus, here the woodpile-structure PC reduces the average curvature,and thus the effective diffraction, by a factor of 10. At this point in k space the curvature K̅ = 0.134 appears exclusively in z-direction, whereas it vanishes in the x-y plane. Thus, for a point source exciting band 1 at Ω = 0.44 we expect self-collimation in the x-y plane and a small divergence in z-direction. On the contrary, for band 2 we find small curvatures at reduced GV in direction close to the space diagonals of the conventional bcc cube.

To verify these predictions we simulated the field propagation of a continuous wave excitation at Ω = 0.44 by means of 3D FDTD calculations. The computing cell consisted of a woodpile photonic crystal with 24 rods in every layer and a height of 40 layers. The rods were oriented along the cartesian axes of the computing cell, i. e., rotated by 45° with respect to the conventional fcc cell and to the IFSs in Fig. 2. In order to avoid the field enhancement of a point source we used a thin wire with Gaussian current distribution along the wire of width w ₀ = 300nm. After a steady state was reached the energy density distribution was monitored and rendered in the 3D cell (see Fig. 3).

 

Fig. 2. IFSs of the high-index woodpile crystal for (a) band 1 at Ω = 0.44, (b) band 2 at Ω = 0.44, (c) band 3 at Ω = 0.69 and (d) band 4 at Ω = 0.73. The curvature K̅ is mapped onto the IFSs (see color bars). The values of the black isolines of lowest curvatures K̅ are 0.2, 0.3, 0.4 for (a), 0.5, 0.6, 0.7 for (b), 0.2, 0.3, 0.4 for (c) and 0.4, 0.5, 0.6 for (d). The red lines show the outline of a square and a hexagonal face of the BZ.

Download Full Size | PDF

 

Fig. 3. Rendered 3D energy density in the woodpile PC with n = 3.4 for Ω = 0.44 (left) and Ω = 0.35 (right). The wire source is located at the origin and oriented along z.

Download Full Size | PDF

Self-collimation in the x-y plane due to the excitation of band 1 can clearly be seen. A closer inspection of the data also reveals a symmetric divergence in z-direction. Other directions of self-collimation out of the plane (see Fig. 3) correspond to the flat IFS regions in band 2. This radiation pattern can be better excited by placing the wire source parallel to x or y.

To verify that this collimation is an effect of the particularly shaped IFS we performed the same calculation at Ω = 0.35 (see Fig. 3, right). Here, instead of self-collimating in definite directions the light diffracts quasi-isotropically.

For the particular diffraction properties do not require a bandgap, it is worthwhile to study the behavior of Bloch waves in 3D PCs made of low-index material suitable for processing via direct laser writing [24]. Thus we performed analogous calculations for n = 1.6 and f = 40%, where a stop band for certain directions just opens up between the second and third band (see Fig. 1). Then the IFS at Ω = 0.73 of the first band has a minimum curvature K̅ = 0.06 at k/Λ = (0.89,0,0.24) and equivalent points (see Fig. 4), with v _g/c = (0.68,0,0.16) leading to a significant splitting of the beam in z-direction (see Fig. 4). However, beams with reduced diffraction can be identified.

4. The inverse opal photonic crystal

A second 3D PC structure of technological impact is the inverted opal structure, because it can be fabricated by inversion of a self-assembled fcc lattice of close packed silica spheres [2,3]. The material filling fraction implied by the crystal’s geometry is f = 1 - π/3 √2 ≈ 26%. The higher symmetry of this structure (full symmetry of the fcc point group) results in a smaller irreducible BZ. For an index of n = 3.4 a bandgap exists between bands 8 and 9 with a size of 4.6%, as can be seen from the bandstructure in Fig. 1. For this crystal we calculated the 3D dispersion relation up to the tenth band. In the following we focus on two IFSs where selfguiding should occur. A guess for the first selfguiding region is already possible from the bandstructure because band 3 is rather flat between X and W as well as X and U in the range 0.53 < Ω < 0.54. This means that there should exist an IFS in band 3, which approximately connects the four corners (W) of the square face of the BZ, the four edge centers (U) and the square center (X) simultaneously. By inspecting the IFS at Ω = 0.54 (see Fig. 5) it can be confirmed that the IFS is indeed almost a flat square connecting these points.

 

Fig. 4. IFS (left) and rendered 3D energy density of the FDTD calculation (right) for the woodpile PC with n = 1.6 for Ω = 0.73. The curvature K̅ is mapped onto the IFS. The values of the black isolines of lowest curvatures K̅ are 0.2, 0.5, 0.6. Note that in the FDTD simulation the structure is rotated by 45° in the x-y plane compared to the conventional fcc cell used for the IFS.

Download Full Size | PDF

Again a 3D FDTD calculation with the same wire source along z as in the woodpile crystal was performed, but for Ω = 0.54. The computing cell consisted of 25 × 25 × 25 conventional cubic cells of the fcc lattice. Four well-collimated beams (see Fig. 6) which propagate in ΓX directions in the x-y plane can be identified. There is no propagating field into z direction because the wire source does not radiate in this direction. A movie of the temporal behavior of H_z in the plane z = 0 in Fig. 7 visualizes the propagation.

On the other hand, it is evident that for Ω = 0.58 the regions of small curvature are in the square face of the BZ for band 3 but in the hexagonal one for band 4, entailing a reduced diffraction around the ΓL direction for band 4. Along the degeneracy lines ΓX and ΓL both IFSs touch.

Because the wire source does not selectively excite a single band we expect a superposition of Bloch waves from both bands in the FDTD simulation. By solely changing the excitation frequency in the FDTD code the volume rendered image shown in Fig. 6 was obtained. The four beams in ΓX directions with eight additional beams in direction of the space diagonals of the fcc cube (ΓL directions) can clearly be recognized.

The IFSs for bands 8 and 9 close to the bandgap are depicted in Fig. 8 as prominent examples for IFSs where even the minimum curvature is large. The individual sheets of the IFSs contract to points located at W for band 8 and at X for band 9 when the frequency approaches the band edges, as can be anticipated already from the bandstructure (see Fig. 1, right). Simultaneously the minimum average curvature of the IFSs tends to infinity and the maximum |v _g| goes to 0, corresponding to standing waves at these exceptional points (W and X). In addition, no other bands exist in the respective frequency range, which otherwise could interfere with the standing waves.

 

Fig. 5. IFSs of the inverse opal structure for (a) band 3 at Ω = 0.54, (b) band 4 at Ω = 0.54, (c) band 3 at Ω = 0.58 and (d) band 4 at Ω = 0.58. The curvature K̅ is mapped onto the IFSs (see color bars). The values of the black isolines of lowest curvatures K̅ are 0.2, 0.3, 0.4. The red lines show the outline of a square and a hexagonal face of the BZ.

Download Full Size | PDF

 

Fig. 6. Rendered 3D energy density in the inverse opal with n = 3.4 for Ω = 0.54 (left) and Ω = 0.58 (right). The wire source is located at the origin and oriented along z.

Download Full Size | PDF

 

Fig. 7. Movie (2.4 MByte) of H_z in the plane z = 0 for the inverse opal at Ω = 0.54.

Download Full Size | PDF

 

Fig. 8. IFSs of the inverse opal structure for (a) band 8 at Ω = 0.78 and (b) band 9 at Ω = 0.84. The curvature K̅ is mapped onto the IFSs (see color bars).

Download Full Size | PDF

5. Conclusions

The dispersion relation in the complete 3D BZ of 3D woodpile-structure for the four lowest order bands and of an inverse opal PCs for the ten lowest order bands has been calculated. By inspecting the resulting IFSs at appropriate frequencies we found domains in k space where the curvature can be minimized and self-collimation in both structures can be expected. These predictions have been verified by numerical experiments using 3D FDTD calculations. Although self-collimation is more pronounced in high-index materials it can also be observed in the low-index counterpart. We also verified that due to the strong dependence of the IFS shape on the frequency the propagation behavior can be drastically altered by slightly changing the operation wavelength.