## Abstract

This work describes a 3-D Finite-Difference Time-Domain (FDTD) computational approach for the optical characterization of an opal photonic crystal. To fully validate the approach we compare the computed transmittance of a crystal model with the transmittance of an actual crystal sample, as measured over the 400 ÷ 750 nm wavelength range. The opal photonic crystal considered has a face-centered cubic (FCC) lattice structure of spherical particles made of polystyrene (a non-absorptive material with constant relative dielectric permittivity). Light-matter interaction is described by numerically solving Maxwell’s equations via a parallelized FDTD code. Periodic boundary conditions (PBCs) at the outer edges of the crystal are used to effectively enforce an infinite lateral extension of the sample. A method to study the propagating Bloch modes inside the crystal bulk is also proposed, which allows the reconstruction of the *ω*-*k* dispersion curve for *k* sweeping discretely the Brillouin zone of the crystal.

© 2014 Optical Society of America

## 1. Introduction

The Finite-Difference Time-Domain (or simply FDTD) method [1,2] is a technique in Computational Electromagnetics which solves Maxwell’s equations directly in the time-domain with a minimum of numerical complexity and computer resources. In this method a 3-D mesh of sampling points in space and a discrete sampling of the time axis for the electric and magnetic vector component values is used, thus getting a space-time grid. A characteristic staggered positioning of the electric and magnetic fields on the grid (i.e. in both space and time) provides second-order accuracy in the partial derivatives, allowing an explicit updating scheme of the field values at any instant of time, resulting in a complete description of their time evolution, starting from a zero field initial configuration. Frequency-domain data are easily obtained by performing an in-line discrete Fourier transform (DFT) of the fields, if the excitation of the system is done by means of signals of finite duration.

In a previous paper [3] we used the FDTD method to characterize the optical behavior of an opal photonic crystal (PC) at a fundamental level, by numerically solving Maxwell’s equations fully accounting for the electromagnetic interaction with its material structure on a macroscopic scale. In that article we considered a small finite crystal sample made of face-centered cubic (FCC) packed polystyrene spheres, illuminated with a light beam impinging perpendicularly to its (1, 1, 1) reticular planes, where use is made here of the Miller indices notation characterizing a given family of planes [4]. The impinging light beam was simulated through a linearly polarized plane wave pulse, injected on the FDTD grid by means of the so called total-field/scattered-field (TF/SF) method. A special boundary is constructed in space on the FDTD mesh and field corrections are applied at any instant of time to the field components on this boundary to maintain the contrast between the interior (*total-field* region) and the exterior (*scattered-field* only region) of the boundary.

Analyzing the scattered far-field obtained by means of the Kirchhoff formula [5, 6] applied to the time Fourier transformed near-field data on a surface surrounding the sample, we numerically determined its transmittance intensity spectrum and traced out the bandgap arising from the crystal symmetry and the refractive index of the polystyrene spheres from which it is assembled. Due to edge effects of the small sample used and to the low dielectric contrast, although the numerical data compared fairly well with the experimental ones for what concerns the wavelength occurrence of the bandgap, this was reproduced with a very small depth. To get a well indented transmittance spectrum a lot of (1, 1, 1) reticular planes along the impinging direction have to be stacked in the numerical model, similarly to what happens in the experimental situation in which a macroscopic crystal sample is used. Actually, it is likely that secondary irradiation from the lateral sample surfaces bends the planar wavefronts thus preventing a full Bragg destructive interference in the forward direction. However, increasing the number of modeled reticular planes requires a lot of computational resources. One of the most efficient solutions consists in the use of periodic boundary conditions (PBCs) [2] on surfaces with normals perpendicular to the beam propagation direction, allowing the modelization of an ideal crystal slab of infinite extension and of given thickness. Moreover, to evaluate its transmittance and reflectance it is sufficient to perform a Poynting vector integration on parallel surfaces before and after the slab (at a sufficient distance to neglect the actual discontinuous crystal structure) of the time Fourier transformed data without the need of a near-to-far-field step.

In the present paper we are thus able to numerically reproduce, simulating the basic principles of the crystal interaction with the light electromagnetic field, the optical behavior of an opal PC as a function of the wavelength in the 400 ÷ 750 nm range in a very satisfactory agreement with experimental measurements. Thanks to the FDTD data results on the crystal bulk and by using spatial Fourier transforms, we are also able to reconstruct quantitatively parts of the *ω* − *k* crystal dispersion curve, a generalizable result which cannot be obtained analytically for arbitrary structures.

For the present work we wrote and used a parallelized version of an in-house FDTD code based on the Message Passing Interface (MPI) library [7–9]. All the calculations where performed on the BGQ, a Southern Ontario Smart Computing Innovation Platform (SOSCIP) BlueGene/Q supercomputer, located at the University of Toronto’s SciNet HPC facility.

The rest of the paper is organized as follows. In Section 2 we give details of the crystal model and on the FDTD simulations. In particular, we describe the way in which parallelization is achieved. In Section 3 we give examples of the calculated field distributions inside the crystal, both in the time and frequency domains. We also compare reflectance and transmittance numerically calculated spectra with a corresponding experimental measurement. In Section 4 we describe the procedure we introduce to reconstruct the dispersion curve of the crystal and show a graph of it. Finally, we summarize our work in Section 5.

## 2. Modeling and simulation

The crystal model is made of polystyrene spheres of relative dielectric permittivity *ε _{r}* = 2.4 which are assumed non-absorptive (

*σ*= 0 Siemens/m). The cell edge of the FCC packing is

*a*= 333.75 nm and the spheres’ radius is

*R*= 118 nm. This gives a structure with a maximal filling factor of 0.74 (i.e., the ratio of the void to the filled space). In the (1,1,1) planes the centers of the spheres lie at the vertices of equilateral triangles of edge length $\ell =a/\sqrt{2}$ and height $h=\sqrt{3}\ell /2$. The FCC packing is obtained by stacking a sequence ABCABC ... of such planes along the

*y*axis. The (1,1,1) planes are then parallel to the

*xz*Cartesian coordinate system. In each plane the triangles are incrementally shifted by a vertical distance of

*d*=

*h*/3 (in B planes the triangles are also laterally shifted by

*ℓ*/2). The various (1, 1, 1) planes are at a distance $D=a/\sqrt{3}$ apart from each other. Fig. 1 depicts examples of sectional views of the crystal in the

*xy*,

*yz*(lateral) and

*xz*(front) coordinate planes, to convey an idea of how the stacking of reticular (1,1,1) planes creates the final 3-D structure of the crystal sample modeled.

For the simulation we use cubic Yee unit cells [1, 2] of edge size *δ* = 2 nm and a time step *δt* = *δ*/2*c _{o}*, to satisfy the Courant-Friedrichs-Lewy stability condition [10],

*c*being the vacuum light velocity. The computational domain truncation is got by means of convolutional perfectly matched layers (CPML) absorbing boundary conditions [11]. However, these are applied only to the upstream and downstream external faces, with respect of the impinging beam propagation direction. On the lateral faces PBCs are applied, which consist in reinjecting the field values in a wraparound fashion. To this end, the interacting crystal structure has to be represented in a consistent manner: the part left out at one side must be mirrored at the other opposite side. To excite the system, we used a “raised cosine” plane wave pulse type with a time profile having compact support like the ones described in [12–14]. The “plane wave injector” no more has lateral and ending guiding surfaces, but only the starting one (the vertical dashed line in Fig. 2). Behind it only the reflected field is present. In that region, if sufficiently far from the crystal slab, reflectance can be evaluated by integrating the real part of the complex Poynting normal vector. The transmittance can be similarly evaluated by considering downstream planes beyond the crystal slab in the total field region (incident plus transmitted ones). By means of the PBCs, the numerical solution obtained is the one corresponding to an infinite extension of the slab in planes perpendicular to the impinging beam. Frequency domain field component values are obtained by updating a DFT made in-line of the time evolving data at each mesh sampling point until excitation has exhausted. The FDTD algorithm is inherently well suited for parallelization. To this end, we decomposed the overall spatial computational FDTD domain into blocks each containing many Yee unit cells, like in Fig. 3. These blocks run concurrently, but not independently from other blocks, on distinct CPU cores. The tangential electric and magnetic field components on the six outer surfaces of each such domain have to be passed among MPI processes running topologically neighbouring blocks, according to Fig. 4, showing schematically a single interface between two contiguous domains, numbered

_{o}*ρ*and

*ρ*+ 1, with a splitted plane

*π*in common. Suitable MPI structured data types are defined to simplify the code accounting for these exchanges. The PBCs correspond to a torus topology, because the first and last blocks of the decomposition along the directions perpendicular to the beam propagation direction have a coinciding face.

## 3. Results

Figure 5 depicts the interaction dynamics of the *z*-polarized light pulse, impinging from the left, with the crystal. It shows animated, sectional, two-dimensional maps of the normalized electric field module distribution when the beam of light propagates down the sample, which in this case is made of only 8 stacked (1, 1, 1) reticular planes.

Numerical results for the calculated transmittance/reflectance are shown in the graph of Fig. 6, which refers to 20 stacked (1, 1, 1) planes, in a sequence ABCABC .... The transmittance *T* and reflectance *R* curves in the graph are calculated independently, by integrating the calculated Poynting vector upstream in the scattered field region and downstream in front of the slab respectively, and not by using the relationship: *R* + *T* = 1, which holds without absorption. It is apparent a bandgap that opens around 540 nm.

The numerical result for the transmittance spectrum of a bigger sample, comprising 100 stacked (1, 1, 1) reticular planes, is reported in Fig. 7 along with the experimental transmittance behavior [15] for a macroscopic sample of the same FCC crystal.

As can be seen, there is a good agreement between the numerical predicted result and the observed one, meaning that there is a complete and effective physical description of the electromagnetic wave interaction with the crystal material structure, including the multiple interference effect from the various reticular planes giving raise to the bandgap at 540 ÷ 550 nm. A ripple appears outside the bandgap, due to minor interference effects, which is tracked by the numerical simulation but not by the instrument, due to its limited resolution bandwidth. The simulation has however the further advantage of being able to calculate the field values inside the crystal on a very small length scale, giving a detailed map of its distribution and accounting for the discontinuous crystal structure.

The use of PBCs avoids any lateral diffractive effects which could affect, weakening it, the transmittance bandgap, as was previously observed using a crystal sample of finite extension. Further, the PCBs allow a saving in the computational resources, because only a small portion in planes perpendicular to the impinging direction has to be considered. We used the minimum transverse extension which could be periodically prolonged. Larger extensions did not give appreciable different results even if they are used in the maps for illustrative purposes. Using the minimum transverse dimensions, in fact, allows time saving.

We also tried the two different linear polarization orientations in the (1, 1, 1) reticular planes (here the *x* and *z* axes) without observing changes in the response, as confirmed also by experimental observations.

Figures 8 and 9 depict sectional, two-dimensional maps, extracted from inside the crystal structure, of the Fourier transformed *E⃗* field distribution. The Fourier transforms have been obtained by means of a DFT on all the spatial values, at the angular frequencies corresponding to the vacuum wavelengths of 500 nm, 550 nm and 600 nm respectively. The data shown refer to the amplitude module, normalized to that of the same frequency component in the incident pulse. The light pulse impinges from the left, and the middle maps in the figures confirm the transmittance bandgap around 550 nm because, as can be inferred from them, the field in the downstream part of the crystal is low compared with the other maps. The field distribution is seen to be confined on the upstream side, i.e., reflected. These graphs were obtained by considering 20 stacked reticular (1, 1, 1) planes.

## 4. Crystal dispersion curve reconstruction

This Section describes an important postprocessing analysis of the calculated numerical data values, aimed to the quantitative reconstruction of the *ω*-*k* dispersion curve for the FCC crystal structure considered, relatively to the GL direction [4] — which in our simulation geometrical setup corresponds to the impinging direction of the light beam — in the *k⃗*-space of its first Brillouin zone.

To this end, an arbitrarily normalized spatial Fourier transform is *numerically* performed on the *frequency domain*, bulk, FDTD calculated data (a spatial Fourier transform for each analyzed frequency), according to the approximation:

*k*= |

*k⃗*| varies discretely from 0 to $\beta =\sqrt{3}(2\pi /a)$, while its direction is held fixed at the GL one. The width of the interval for

*k*corresponds to the length of the L path, i.e. the GL segment for the FCC lattice with edge size

*a*. Other directions could as well be set, to analyze the S or GK path, the D or GX path, or others.

Here the ∼ symbol denotes frequency domain values, *r⃗ _{ℓ,m,n}* denotes the position vector of a point in space according to the FDTD discrete sampling with indices

*ℓ*,

*m*,

*n*along the Cartesian axes, and

*i*in the exponential is the imaginary unit.

Only the electric field component corresponding to the polarization direction (here the *z* axis) of the incident exciting signal has been considered, although components along other Cartesian axes are also excited within the crystal, but with considerable less intensity.

The necessarily finite integration domain is at the center of the crystal sample, embracing some FCC cells along the beam propagation direction, but with the boundaries well apart from those of the sample, to avoid contributions from unwanted reflected components. This implies some negligible truncation errors in the tails of the calculated spectra. The integration domain may instead be extended along the directions perpendicular to the faces at which periodic boundary conditions are imposed, until they are meet, without any truncation effects.

The mode (Bloch wave) supported for a given angular frequency *ω* is identified from the corresponding *k*-spectrum as a dominant peak (Fig. 10), which is a normalization independent operation. What we are really interested in is not the peak height, but the *k* value of the peak occurrence. Such peak detection has been made finer by fitting the points around the identified peak through a spline interpolation, to better localize the value at the maximum.

By comparing the relative height of the peaks allows us also to recognize the bandgap along the abscissas *k*-axis, because a non-propagating mode has a lower amplitude spectrum, compared with those of the propagating ones. On the other hand, the (angular) frequency *ω* occurrence of the bandgap is also already known from the transmittance curve, like the one in Fig. 6 or Fig. 7. There, the abscissas are labelled as vacuum wavelengths for more practical insight, but actually correspond univocally to frequency domain data: *ω* = 2*πc/λ*. Examples of numerically calculated amplitude *k*-spectra are given in Fig. 10 and Fig. 11, the latter showing together the graphs of the previous one. Note that, as expected, at the angular frequency corresponding to *λ* = 540 nm (the middle graph), i.e. the frequency of the first bandgap (see Fig. 6), the magnitude of the peak is less than one half the magnitudes of the other peaks (left and right graphs in the figure) which correspond to fully propagated modes.

We eventually end up with a sequence of (*k*, *ω*) pairs which, when graphed, allows the quantitative reconstruction of the dispersion curve searched for. That in Fig. 12 is for a stack of 20 reticular (1, 1, 1) planes, in a sequence ABCABC .... The *k* and *ω* axes of the graph are normalized: the horizontal one corresponds to *k/β* values, the vertical one to *ωa*/2*πc* values. To get the second branch of the curve above the first bandgap (the so called “air-line”), the relevant *k* data obtained by means of the spatial Fourier transforms have been mirrored in the vertical line through the value *k _{bg}* = 0.499

*β*corresponding to the first bandgap.

It is apparent from the figure that the slope flattens at the interval extrema on the *k*-axis — except at the origin — as expected. The four dots standing at the most right position in the graph of Fig. 12 correspond to the bandgap and should be neglected. The width of the bandgap here corresponds to that of Fig. 6.

The behavior of the curve toward the origin is, from a linear regression analysis, that of a straight line with zero intercept. The slope of the straight line allows to calculate the effective refractive index of the crystal. We get a numerically predicted value of *n _{eff}* = 1.393; the experimental value of

*n*= 1.436 reported by [15] and that of

_{eff}*n*= 1.427 from geometric mean consideration in effective medium approximation there cited are in fairly good accordance with ours.

_{eff}This is a general procedure usable for different packing geometries and which may be also extended to the higher branches of the dispersion curve and for the entire Brillouin zone.

## 5. Conclusion

In this paper we demonstrate the effectiveness and robustness of a parallel FDTD code with periodic boundary conditions, which we developed to characterize numerically the optical physical behavior of an opal photonic crystal. In particular, we applied this analysis to the quantitative reconstruction of the dispersion curve, a task which cannot be accomplished analytically. The described procedure can be generalized and we are confident, in future, to apply this simulation tool also to the random close packing sphere configuration for the analysis of photonic glasses.

## Acknowledgments

This work is supported by IBM Canada Research and Development Center.

## References and links

**1. **K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. **14**, 302–307 (1966). [CrossRef]

**2. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3rd ed. (Artech House, 2005).

**3. **A. Vaccari, L. Cristoforetti, A. Calà Lesina, L. Ramunno, A. Chiappini, F. Prudenzano, A. Bozzoli, and L. Calliari, “Parallel FDTD modeling of an opal photonic crystal,” Opt. Eng. **53**(7), 071809 (2014). [CrossRef]

**4. **J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, *Photonic Crystals. Molding the Flow of Light*, 2nd ed. (Princeton University, 2008), Appendix B.

**5. **J. A. Stratton, *Electromagnetic Theory* (McGraw-Hill, 1941), Chap. VIII.

**6. **J. D. Jackson, *Classical Electrodynamics*, 3rd ed. (Wiley, 1998), Chap. 9.

**7. **W. Gropp, E. L. Lusk, and A. Skjellum, *Using MPI. Portable Parallel Programming With the Message Passing Interface*, 2nd ed. (MIT, 1999).

**8. **W. Gropp, E. L. Lusk, and R. Thakur, *Using MPI-2. Advanced Features in the Message Passing Interface* (MIT, 1999).

**9. **P. Pacheco, *Parallel Programming with MPI* (Morgan Kaufmann, 1997).

**10. **A. Taflove and M. E. Brodwin, “Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell’s equations,” IEEE Trans. Microwave Theory Tech. **23**, 623–630 (1975). [CrossRef]

**11. **J. A. Roden and S. D. Gedney, “Convolution PML (CPML): An efficient FDTD implementation of the CFS-PML for arbitrary media,” Microw. Opt. Tech. Lett. **27**(5), 334–339 (2000). [CrossRef]

**12. **R. Pontalti, L. Cristoforetti, and L. Cescatti, “The frequency dependent FD-TD method for multi-frequency results in microwave hypertermia treatment simulation,” Phys. Med. Biol. **38**, 1283–1298 (1993). [CrossRef]

**13. **R. Pontalti, J. Nadobny, P. Wust, A. Vaccari, and D. Sullivan, “Investigation of static and quasi-static fields inherent to the pulsed FDTD method,” IEEE Trans. Microwave Theory Tech.50(8), (2002). [CrossRef]

**14. **A. Vaccari, R. Pontalti, C. Malacarne, and L. Cristoforetti, “A robust and efficient subgridding algorithm for finite-difference time-domain simulations of Maxwell’s equations,” J. Comput. Phys. **194**, 117–139 (2004). [CrossRef]

**15. **A. Chiappini, C. Armellini, A. Chiasera, M. Ferrari, L. Fortes, M. C. Gonçalves, R. Guider, Y. Jestin, R. Retoux, G. N. Conti, S. Pelli, R. M. Almeida, and G. C. Righini, “An alternative method to obtain direct opal photonic crystal structures,” J. Non-Cryst. Solids **355**, 1167–1170 (2009). [CrossRef]