Accurate calculation of the local density of optical states in inverse-opal photonic crystals
aa r X i v : . [ phy s i c s . op ti c s ] F e b Accurate calculation of the local density of optical states ininverse-opal photonic crystals
Ivan S. Nikolaev , Willem L. Vos , , and A. Femius Koenderink , ∗ Center for Nanophotonics, FOM Institute for Atomic and Molecular Physics (AMOLF), Kruislaan 407, NL-1098SJAmsterdam, The Netherlands Complex Photonic Systems, MESA + Institute for Nanotechnology, University of Twente, The NetherlandsCorresponding author: [email protected]
Submitted to OSA Dec 23, 2008We have investigated the local density of optical states (LDOS) in titania and silicon inverse opals –three-dimensional photonic crystals that have been realized experimentally. We used the H-field plane-waveexpansion method to calculate the density of states and the projected local optical density of states, whichare directly relevant for spontaneous emission dynamics and strong coupling. We present the first quantitativeanalysis of the frequency resolution and of the accuracy of the calculated local density of states. We havecalculated the projected LDOS for many different emitter positions in inverse opals in order to supply atheoretical interpretation for recent emission experiments and as reference results for future experiments andtheory by other workers. The results show that the LDOS in inverse opals strongly depends on the crystallattice parameter as well as on the position and orientation of emitting dipoles.
OCIS codes:
1. Introduction
Photonic crystals are metamaterials with periodic varia-tions of the dielectric function on length scales compara-ble to the wavelength of light. These dielectric compos-ites are of a keen interest for scientists and engineersbecause they offer exciting ways to manipulate pho-tons [1, 2]. Of particular interest are three-dimensional(3D) photonic crystals possessing a photonic bandgap, i.e. a frequency range where no photon modes exist atall. Photonic bandgap materials possess a great potentialfor drastically changing the rate of spontaneous emissionand for achieving localization of light [3–6]. Control overspontaneous emission is important for many applicationssuch as miniature lasers [7], light-emitting diodes [8] andsolar cells [9]. According to Fermi’s golden rule, the rateof spontaneous emission from quantum emitters such asatoms, molecules or quantum dots is proportional to the‘local radiative density of states’ (LDOS) [10, 11], whichcounts the number of electromagnetic states at given fre-quency, location and orientation of the dipolar emitters.In addition, nonclassical effects can occur that are be-yond Fermi’s Golden Rule [2, 11]. These strong couplingphenomena, such as fractional decay, rely on coherentcoupling of the quantum emitter to a sharp feature ina highly structured LDOS. Whether these effects are infact observable in real photonic crystals depends on howrapid the LDOS changes within a small frequency win-dow [12]. Therefore local density of states calculationsfor experimentally realized photonic crystals are essen-tial to assess current spontaneous experiments and fu-ture strong coupling experiments alike, provided that thecalculations are accurate and have a controlled frequencyresolution. LDOS effects on spontaneous emission in photoniccrystals have been experimentally demonstrated in avariety of systems. Since only three-dimensional (3D)crystals promise full control over all optical modes withwhich elementary emitters interact, many groups havepursued 3D photonic crystals. Fabrication of such pe-riodic structures with high photonic strengths is, how-ever, a great challenge [2, 13]. Inverse opals are amongthe most strongly photonic 3D crystals that can be fabri-cated relatively easily using self-assembly methods. Suchcrystals consist of fcc lattices of close-packed air spheresin a backbone material with a high dielectric constant[14–18]. In these inverse opals, continuous-wave exper-iments on light sources with a low quantum efficiencyrevealed inhibited radiative emission rates [19, 20]. En-hanced and inhibited time-resolved emission rates havebeen observed from highly-efficient emitters in 3D in-verse opals [21]. Simultaneously, several groups have re-alized that emission enhancement and partial inhibitioncan also be obtained in 2D slab structures [22–26].Following the first calculations by Suzuki and Yu [27],several papers have reported on calculations of the lo-cal density of optical states in photonic crystals usingboth time domain [28–30] and frequency domain meth-ods [31–34]. Unfortunately, analysis of experimental datafor emitters in 3D crystals in terms of these LDOS calcu-lations is compounded by several problems in literature:1. Most prior calculations were performed for modelsystems that don’t correspond to structures used inexperiments, and for emitter positions that are notprobed in experiments.2. The accuracy of the reported LDOS has never been1iscussed, hampering comparison to experiments.3. The frequency resolution of the reported LDOS hasremained unspecified. Therefore sharp features ofrelevance for non-classical emission can not be as-sessed [11, 12].4. Many previously reported LDOS calculations are er-roneous for symmetry reasons, as pointed out byWang et al. [33].In this paper we aim to overcome all these problems.We benchmark the accuracy and frequency resolution forour LDOS calculations to allow quantitative comparisonwith experiments and with quantum optics strong cou-pling requirements. Since prior LDOS calculations arescarce (and partly erroneous, item 1 above) we presentsets of LDOS’s for experimentally relevant structures,and for spatial positions where sources can be practi-cally placed. Specifically, we model the spatial distribu-tion of the dielectric function ǫ ( r ) in such a way that itclosely resembles ǫ ( r ) in titania (TiO ) [14, 15] and sili-con (Si) inverse-opal photonic crystals [17, 18]. For thesetwo structures we calculated the LDOS at various posi-tions in the crystal unit cell and for specific orientationsof the transition dipoles. The results on the TiO inverseopals are relevant for interpreting recent emission exper-iments [21,35]. To aid other workers to interpret their ex-periments and to benchmark their codes, we make datasets that we report in figures throughout this manuscriptavailable as online material.The paper is arranged as follows: in Section 2, wepresent a detailed description of the method by whichwe have calculated the photonic band structures and theLDOS. We discuss the accuracy and frequency resolutionof our calculations. In Section 3 we compare our compu-tations with the known DOS in vacuum and with earlierresults on the DOS and LDOS [33] in 3D periodic struc-tures. Section 4 describes the LDOS in inverse opals fromTiO , and in Section 5 we present results of the LDOSin Si inverse opal photonic band gap crystals.
2. Calculation of local density of states
A. Introduction
The local radiative density of optical states is defined as N ( r , ω, e d ) = 1(2 π ) X n Z BZ d k δ ( ω − ω n, k ) | e d · E n, k ( r ) | , (1)where integration over k -vector is performed over thefirst Brillouin zone, n is the band index and e d is theorientation of the emitting dipole. The total density ofstates (DOS) is the unit-cell and dipole-orientation av-erage of the LDOS defined as N ( ω ) = P n R BZ d k δ ( ω − ω n, k ). The important quantities that determine theLDOS are the eigenfrequencies ω n, k and electric fieldeigenmodes E n, k ( r ) for each k -vector. Calculation ofthese parameters will be discussed below in Section B. The expression for the LDOS contains the term | e d · E n, k ( r ) | that depends on the dipole orientation e d . It isimportant to realize that in photonic crystals the vec-tor fields E n, k ( r ) are not invariant under the latticepoint-group operations α , as first reported in Ref. [33].Explicitly, this means that the projection of the fieldof a mode at wave vector k on the dipole orientation e d ( i.e. | e d · E n, k ( r ) | )is not identical to the projection | e d · E n,α [ k ] ( r ) | of the symmetry related modes with wavevectors α [ k ] on the same e d . As a consequence one cannot calculate the LDOS for a specific dipole orientationby restricting the integral in Eq. (1) over the irreduciblepart (1/48th) of the Brillouin zone, since symmetry re-lated wave vectors do not give identical contributions.Unfortunately, in many previous reports on the LDOS,this reduced symmetry for vector modes as compared toscalar quantities was overlooked, resulting in erroneousresults [33]. In general, the only symmetry that can beinvoked to avoid using the full Brillouin zone for LDOScalculations is inversion symmetry, which corresponds totime reversal symmetry. Consequently, correct results re-quire that exactly half of the Brillouin zone is consideredfor LDOS calculations, rather than the irreducible partof the Brillouin zone that was used in most previous lit-erature. We have explicitly verified that our implemen-tation (using the k z ≥ B. Plane-wave expansion
We use the H-field inverted plane wave expansionmethod [31, 36, 37] to solve for the electromagnetic fieldmodes in photonic crystals. For nonmagnetic materials,it is most convenient to solve the wave equation for the H ( r ) field [38] ∇ × (cid:2) ǫ ( r ) − ∇ × H ( r ) (cid:3) = ω c H ( r ) . (2)because the operator ∇ × ǫ ( r ) − ∇× is Hermitian, andconsequently has real eigenvalues ω /c [1,31,36,37]. Be-cause of the periodicity of the dielectric function ǫ ( r ) inphotonic crystals, the field modes H k ( r ) of the eigen-value problem Eq. (2) satisfy the Bloch theorem [39]: H k ( r ) = e i k · r u k ( r ) . (3)These Bloch modes are fully described by the wavevector k and the periodic function u k ( r ), which has the period-icity of the crystal lattice so that u k ( r ) = u k ( r + R ). Tosolve the wave equation, the inverse dielectric functionand the Bloch modes are expanded in a Fourier seriesover the reciprocal-lattice vectors G : ǫ ( r ) − = η ( r ) = X G η G e i G · r and (4)2 k ( r ) = X G u kG e i ( k + G ) · r , (5)where η G and u kG are the 3D Fourier expansion coeffi-cients of respectively η ( r ) and u k ( r ). Substituting theseexpressions into the H-field wave equation in Eq. (2), weobtain a linear set of eigenvalue equations: − X G ′ η G − G ′ ( k + G ) × [( k + G ′ ) × u n, kG ′ ] = ω n ( k ) c u n, kG . (6)This infinite equation set with the known parameters G and η G − G ′ determines all allowed frequencies ω n ( k ) foreach value of the wave vector k , subject to the transver-sality requirement ∇ · H k ( r ) = 0. Due to the periodicityof u k ( r ), we can restrict k to the first Brillouin zone. Foreach wave vector k , there is a countably infinite num-ber of modes with discretely spaced frequencies. All themodes are labeled with the band number n in order ofincreasing frequency and are described as a family ofcontinuous functions ω n ( k ) of k .To compute the eigenfrequencies ω n ( k ) and the ex-pansion coefficients of the eigenmodes u n, kG , the infiniteequation set is truncated. By restricting the number ofreciprocal-lattice vectors G to a finite set G with N G ele-ments, Eq. (6) is limited to a 3 N G dimensional equationset. In our implementation we choose the truncated set G to correspond to the set of all reciprocal lattice vectorswithin a sphere centered around the origin of k-space.The transversality of the H-field gives an additional con-dition on the eigenmodes: ( k + G ) · u n, kG = 0, which elim-inates one vector component of u n, kG . Following Ref. [31],for each k + G one needs to find two orthogonal unit vec-tors e , k + G that form an orthogonal triad with k + G . Byexpressing the eigenmode expansion coefficients in theplane normal to k + G as u n, kG = u n, kG , e k + G + u n, kG , e k + G ,we remove one third of the unknowns. Then, Eq. (6) be-comes X G ′ ∈G η G − G ′ | k + G || k + G ′ | · "(cid:18) e k + G · e k + G ′ − e k + G · e k + G ′ − e k + G · e k + G ′ e k + G · e k + G ′ (cid:19) u n, kG ′ , u n, kG ′ , ! = ω n ( k ) c u n, kG , u n, kG , ! , ∀ G ∈ G . (7)To find the matrix of Fourier coefficients η G − G ′ , we usedthe method of Refs. [31, 36]. The coefficients are com-puted by first Fourier-transforming the dielectric func-tion ǫ ( r ), and then truncating and inverting the resultingmatrix. As first noted by Ho, Chan and Soukoulis [36],using the inverse of ǫ G − G ′ rather than η G − G ′ dramati-cally improves the (poor) convergence of the plane-wavemethod that is associated with the discontinuous natureof the dielectric function [37]. A rigorous explanationfor this improvement was put forward by Li [40], who studied the presence of Gibbs oscillations in the trun-cated Fourier expansion of products of functions withcomplementary jump discontinuities. Using the H-fieldinverted matrix plane wave method, the frequencies ob-tained with N G = 725 (for fcc structures) deviate byless than 0.5 % from the converged band structures [31].Solving Eqs. (7) gives the frequencies ω n ( k ) and theFourier expansion coefficients for the H-field eigenmodes H n, k ( r ) needed to calculate the LDOS in the photoniccrystal. The required E-fields E n, k ( r ) are obtained usingthe Maxwell equation ∂ D /∂t = ∇ × H : E n, k ( r ) = 1 ω n ( k ) ǫ X G , G ′ ∈G η G ′ − G | k + G | · (cid:20)(cid:16) u n, kG , e k + G − u n, kG , e k + G (cid:17) e i ( k + G ′ ) · r (cid:21) . (8)From the orthonormality of the eigenvectors of Eq. (7)it follows that the Bloch functions H n, k ( r ) and E n, k ( r )defined above satisfy the orthonormality relations: Z BZ H n, k ( r ) · H ∗ n ′ , k ′ ( r ) d r = δ ( k − k ′ ) δ n,n ′ , (9) Z BZ ǫ ( r ) E n, k ( r ) · E ∗ n ′ , k ′ ( r ) d r = δ ( k − k ′ ) δ n,n ′ . (10)It should be noted in the definition of E n, k that the mul-tiplication by 1 /ǫ ( r ) to calculate E from D is not in thedenominator in front of the Fourier expansion. Rather itappears as a matrix multiplying the D -field in Fourierspace, i.e., within the sum over reciprocal lattice vec-tors. This ordering ensures that complementary jumpsin D and 1 /ǫ ( r ) cancel, even for the truncated Fourierseries, as can be easily checked by plotting calculatedmode profiles for TM modes in 2D crystals. C. Frequency resolution and accuracy of LDOS
We are not aware of any previous report that bench-marks the accuracy of the calculated local radiative den-sity of states or that specifies the frequency resolution.Motivated by the requirements for accurate results forcomparison with experiments and for judging the utilityof crystals for non-classical emission dynamics, we con-sider the accuracy and resolution of our approximationto the LDOS integral in Eq. (1). The main approxima-tion is to replace the integration over wave vector k byan appropriately weighted summation over a discrete setof wave vectors on a discretization grid [41]. Either in-terpolation schemes [31] or simple histogramming (‘root-sampling’) methods are used to compute the LDOS. Thek-point grid density sets the number of k-points in theBrillouin zone, N k . The accuracy of the resulting LDOSapproximation is set by the the density of grid pointsthat is used to discretize the wave vector integral. Dueto the transparent relation between the accuracy of theDOS, the frequency resolution and the k -vector samplingresolution, we will focus on the simple histogramming3ig. 1. (Color online) DOS per volume in units 4 /a c forvacuum modeled as an empty fcc crystal. The calculatedDOS shown by red histogram bars is compared to theanalytically derived ω behavior (curve). In vacuum theDOS per volume equals the dipole-averaged LDOS.Fig. 2. Average absolute deviation of the calculated DOSfrom the exact total DOS N ( ω ) for an ‘empty crystal.’The average runs over the frequency range 0 < ω < πc/a and the deviation is in units of the DOS N ( ω )per volume at ω a/ πc = 1, i.e., in units 4 /a c . In accor-dance with Eq. (11), the error is inversely proportionalto the ratio ∆ ω/ ∆ k of the histogram bin width ∆ ω tothe integration grid spacing ∆ k . Symbols correspond tointegration using the number of k-points N k = 280, 770,1300, 2480, 2992 and 3570 ( (cid:3) , (cid:4) , ◦ , • , ♦ , (cid:7) ) respectively,with various ∆ ω .approach. For the LDOS computations one chooses acertain frequency bin width ∆ ω to build an LDOS his-togram. For a desired frequency resolution ∆ ω and adesired accuracy for the LDOS content N ( r , ω, e d )∆ ω ineach frequency bin, one needs to choose an appropriate k -vector spacing ∆ k that depends on the steepness ofthe sampled dispersion relation ω ( k ). Indeed, the usefulfrequency resolution ∆ ω of a histogram of the DOS andLDOS is limited by the resolution ∆ k of the grid in k space to be ∆ ω ≈ ∆ k |∇ k ω | , (11)as detailed in Ref. [42] for the electronic DOS. This cri-terion relates the separation between adjacent k -gridpoints to their approximate frequency spacing via thegroup velocity. If histogram bins are chosen too narrowcompared to the expected frequency spacing betweencontributions to the discretized LDOS integral, unphysi-cal spikes appear in the approximation, especially in thelimit of small ω , where the group velocity |∇ k ω | is usu-ally largest. Apart from full gaps in the LDOS, photoniccrystals also promise sharp lines at which the LDOS isenhanced, which are important for non-classical emis-sion dynamics. Hence it is especially important to dis-tinguish sharp spikes that are due to histogram binningnoise from true features. Unfortunately, many reports inliterature feature sharp spikes that are evidently binningnoise (wave vector undersampling errors), as they occurin the long wavelength limit, below any stop gap.To improve the resolution without adding time-consuming diagonalizations, several interpolationschemes have been suggested [42]. Within the his-togramming approach, an interpolation scheme essen-tially improves binning statistics by adding histogramcontributions from intermediate grid points on theassumption that quantities vary linearly between gridpoints. While interpolation decouples the binning noisefrom the frequency bin width ∆ ω , resulting in arbitrarilysmooth LDOS approximations, it has the disadvantageof obscuring the intrinsic relation between frequencyresolution and wave vector resolution in Eq. (11).A good benchmark for the binning noise of the k -spaceintegration, independent of the convergence of the plane-wave method, is to calculate the LDOS or DOS of an‘empty’ crystal, with uniform dielectric constant equalto unity. Such an ‘empty’ crystal represents a limit ofzero photonic strength and the maximum possible groupvelocity |∇ k ω | . In Figure 1 we show the DOS in an empty fcc crystal. As expected for a crystal with zero dielectriccontrast, the calculated DOS is independent of dipole po-sition and orientation. In agreement with the well-knownDOS in vacuum the calculated DOS increases paraboli-cally with frequency. Fluctuations around the parabolaare due to binning noise associated with the finite k-space discretization.We have calculated the relative root-mean-square er-ror in the calculated density of states (DOS) for an fcc ‘empty’ crystal averaged over all histogram binsin the frequency range 0 < ωa/ πc < ω/ ∆ k . In most cases, one wants to calculate the LDOSwith a given frequency resolution ∆ ω , i.e. N ( r , ω, e d )∆ ω ,to within a predetermined accuracy. For instance, calcu-lating the vacuum DOS for frequencies 0 < ωa/ πc < . · /a c (1% of the vacuum DOS at ωa/ πc = 1) and a desiredfrequency resolution of ∆ ω = 0 . πc/a ), requires using∆ k ∼ ∆ ω/ . c . The number of k-points correspondingto this wave vector sampling equals 2480 k-points in theirreducible wedge of the fcc Brillouin zone, 59520 in therequired half Brillouin zone, or equivalently N k = 119040in the full Brillouin zone. Photonic crystals with nonzeroindex contrast cause a pronounced frequency structureof the LDOS. Due to the flattening of bands comparedto the dispersion bands of vacuum-only crystals, the k-space integration itself is at least as accurate, assumingthat the eigenfrequencies and the field-mode patterns areknown with infinite accuracy. In practice, one needs toadjust the number of plane waves to obtain all eigen-frequencies to within the desired frequency resolution∆ ω . In our computations for fcc crystals, we representedthe k-space of a half of the first Brillouin zone by anequidistant grid consisting of N k / ω = 0 . · πc/a , and the spatial resolution of theLDOS is approximately a/
40 judging from the maximum | G − G ′ | ≈ /a involved in the expansion with N G =725 plane waves in Eq. (4). D. Computation time required for LDOS
Calculating the LDOS in 3D periodic structures is atime-consuming task. The chosen degree of k-space dis-cretization ( N k / N G = 725) are the result of a trade-offbetween the desired accuracy and tolerated duration ofthe calculations. Essentially the computation consists ofsolving for the lowest n eigenvalues and eigenvectors ofa real symmetric 2 N G × N G matrix for each of the N k k-points independently. To accomplish this calculation,we use the standard Matlab ”eigs” implementation [43]of an implicitly restarted Arnoldi method (ARPACK)that takes approximately 2.2 seconds to find the low-est 20 eigenvalues and eigenvectors of a single Hermitian1500x1500 (i.e., N G ≈
750 plane waves) matrix on a3 GHz Intel Pentium 4 processor, or on a 2.4 GHz In-tel Core Duo processor. For 145708 k-points the result-ing computation time is hence on the order of 90 hours(4 days). Since the Matlab ARPACK routine is alreadyhighly optimized, we do not expect that the computationtime per k-point can be significantly reduced for LDOScalculations based on the standard H-field inverted ma-trix plane wave method. In terms of the number of planewaves N G , the computation time scales as N G . As inRef. [44], the algorithm can be accelerated by realizingthat the iterative eigensolver does not require the fullmatrix H multiplying the vector u in Eq. 7, but rathera function that quickly computes the image H u of anytrial vector u . Since calculating H u repeatedly involvesa slow matrix-vector multiplication with η G − G ′ , the al-gorithm is only accelerated for N G > fcc crys-tal consisting of spheres with ǫ = 7.35 in a mediumwith ǫ = 1.77 with a filling fraction of the spheres of25 vol %. The solid dotted curve represents calculationsfrom Ref. [31]. Our result is plotted as a histogram (red).
3. Comparison with previous results
To test the computations, we compare our results withearlier reports. We have calculated the DOS and LDOSin an fcc crystal consisting of dielectric spheres with ǫ = 7.35 (TiO ) in a medium with ǫ = 1.77 (water) –the same structure as was analyzed in Refs [31, 33]. Thespheres occupy 25 vol% of the crystal. Figure 3 shows thetotal DOS in this photonic crystal calculated by us andby Busch and John [31] . We reproduce the earlier cal-culations of the total DOS: both results shown in Fig. 3are in excellent agreement, with deviations less than 2%throughout the frequency range 0 < ωa/ πc <
1. In Fig-ure 4 (Media 1)we demonstrate the LDOS in the samephotonic crystal at a specific location: at a point equidis-tant from two nearest-neighbor spheres. In this calcu-lation, we used the same number of reciprocal-latticevectors N G = 965 as in the only available benchmarkpaper by Wang et al. [33] that does not contain sym-metry errors. We find that our calculations (empty cir-cles) are in good agreement up to a/λ = 0 .
85 with theLDOS reported previously (solid curve): the deviationsare smaller than 1%. Deviations at higher frequenciesare either due to a difference in accuracy of k-spaceintegration (frequency binning resolution and k-spacesampling density), or to a difference in accuracy of theplane wave methods (eigenmodes and eigenfrequencies).Unfortunately, neither the accuracy nor the frequencyresolution is specified in Ref. [33]. The k-space samplingdensity in our work is approximately twice the densityspecified in Ref. [33]. Based on the excellent agreementof our total DOS calculations with those of Busch andJohn [31], and on the fact that Wang et al. used a k-spacedensity and plane wave number comparable to that ofBusch and John, it is unlikely that inaccuracy of the k-space integration in either calculation is the source of dis-crepancy. We therefore surmise that deviations are dueto a difference in the method of evaluation of E n, k ( r ).Conversion of D to E by multiplication with 1 /ǫ ( r ) in5ig. 4. (Color online) Dipole-averaged LDOS in the samephotonic crystal as in Fig. 3 at a position ( , , a/λ = 0.495.real space as proposed in [33] can cause incorrect modeamplitudes due to Gibbs oscillations, as opposed to theFourier space conversion using Eq. (8). In general ourcalculations confirm the result by Wang et. al. [33] thatthe LDOS is only correctly calculated by integration overhalf the Brillouin zone, rather than over the irreduciblepart as was incorrectly used in earlier reports.
4. LDOS in TiO inverse opals In recent time-resolved experiments, enhanced and in-hibited emission rates were demonstrated for quantumdots embedded inside strongly photonic TiO inverseopals [19,21,35]. In the framework of these experiments,it is highly relevant to calculate the LDOS inside suchinverse-opal photonic crystals, especially for the sourcepositions occupied in experiments. We model the posi-tion dependence of the dielectric function ǫ ( r ) as shownin Figure 5. This model assumes an infinite fcc latticeof air spheres with radius r = 0 . √ a ( a is the cubiclattice parameter). The spheres are covered by overlap-ping dielectric shells ( ǫ = 6 .
5) with outer radius 1.09 r .Neighboring air spheres are connected by cylindrical win-dows of radius 0.4 r . The resulting volume fraction ofTiO is equal to about 10.7%. These structural param-eters are inferred from detailed characterization of theinverse opals using electron microscopy and small an-gle X-ray scattering [14, 15]. Moreover, the stop gaps inthe photonic band structure (Fig. 6 (Media 2 and Me-dia 3)) calculated using this model agree well with re-flectivity measurements in the ranges of both the first-order ( a/λ = 0 .
7) and second-order Bragg diffraction( a/λ = 1 .
2) [45].Henceforth, we will consider the relative LDOS, whichis the ratio of the LDOS in a photonic crystal to that ina homogeneous medium with the same volume-averageddielectric function. This scaling is motivated by exper-imental practice, in which emission rate modificationsare judged by normalizing measured rates to the emis-sion rate of the same emitter in crystals with muchsmaller lattice constant a [19, 21, 35]. In these reference Fig. 5. Rendering of the dielectric function in one fcc unitcell that models the TiO inverse-opal structure in sec-tion 4: an fcc lattice of air spheres of radius r = 0 . √ a with a being the cubic lattice parameter. The spheresare covered by shells with ǫ = 6 . r . Neighboring air spheres are connected by win-dows of radius 0.4 r . The letters ( a – d ) indicate four dif-ferent positions at the TiO -air: a = (1 , , / (2 √ b =(1 , , / (4 √ c = (1 , , / (2 √
6) and d = (0.33,0.13,0)(points shown are symmetry-equivalents). The dash-dotted line shows the main diagonal of the cubic unitcell.Fig. 6. Left : Photonic band structure (Media 2)for theTiO inverse opal shown in Fig. 5. The grey rectanglesindicate stopgaps in the ΓL direction and one stopgapin the ΓX direction for the inverse opal. The stopgapsresult in the decreased DOS ( right : circles, Media 3) atcorresponding frequencies compared to the DOS in a ho-mogeneous medium with n av = 1 .
27 ( right : solid line).6ig. 7. Relative LDOS in the inverse opal shown in Fig 5at three different positions: (a, Media 4) r = (0,0,0) [thecenter of an air sphere, solid curve], r = (1 , ,
1) [amongthree air spheres, dash-dotted curve] and r = ( ,0,0)[midway between two spheres along [1,0,0] direction, dot-ted curve]; (b, Media 5) r = (1 , ,
0) [in the windowbetween two spheres] projected on [1,1,0], [-1,1,0] and[0,0,1] directions shown by solid, dash-dotted and dot-ted curves, respectively.systems, emission frequencies correspond to the effectivemedium limit a/λ ≪ . n av = √ ǫ av = 1 .
27 for the TiO crystals [19]. In unitsof 4 /a c , the LDOS in a homogeneous medium is equalto n av ( a/λ ) /
3. In Figure 7a (Media 4) we plot the re-sulting LDOS at three positions in the unit cell: at r =(0,0,0), (1 , ,
1) and ( ,0,0). Due to the high symmetryof these points, the LDOS does not depend on the dipoleorientation, as we verified explicitly. A first main obser-vation from Fig. 7a (Media 4) is that the LDOS differsconsiderably between these three positions at all reducedfrequencies. This observation illustrates the well-knownstrong dependence of the LDOS on position within theunit cell of photonic crystals [10, 27, 31].A second main observation in Fig. 7a (Media 4) is thatthe LDOS strongly varies with reduced frequency, reveal-ing troughs and peaks caused by the pseudogap near a/λ = 0 .
7, which is related to 1 st -order stop gaps suchas the L-gap. The effects of 2 nd -order stopgaps appearbeyond a/λ > .
15 [45]. In the middle of the air regionat position r = (0,0,0) there is a sharp, factor-of-threeenhancement at a/λ ≈ .
35 within a narrow frequencyrange. This feature could be probed by resonant atomsinfiltrated in the crystals [46]. At position r = ( ,0,0) inan interstitial, the mode density has a broad trough near a/λ = 1 .
25 that will lead to strongly inhibited emission.A third main observation is that at spatial positionswith low symmetry, the LDOS clearly depends on theorientation of the transition dipole moment. Figure 7b(Media 5) shows the frequency dependent LDOS forthree perpendicular dipole orientations at a position inthe center of a window that connects two neighboringair spheres (see Figure 5). The LDOS differs for all threeorientations, and is thus anisotropic. At low frequency( a/λ = 0 . , ,
0) direction, with increasing fre-quency the highest rate shifts to the (0 , ,
1) and then tothe ( − , ,
0) orientation. We emphasize that for emit-ters with fixed or slowly varying dipole orientations, such Fig. 8. Relative LDOS at four key frequencies, a/λ =0.535, 0.725, 0.865 and 1.295, in the inverse opal as afunction of position r from (0,0,0) in the [1,1,1] direction.The hatched boxes indicate the position of the dielectricTiO shell. The LDOS is projected on two dipole ori-entations: (a) [1,1,1] perpendicular to the dielectric-airinterface, and (b) [-1,1,0] parallel to the interface. TheLDOS projected on the [-1,-1,2] and [-1,1,0] directions isequal. For r/a ∈ [ √ , √
3] the LDOS is mirror-symmetricto that in the region [0 , √ ].as dye molecules or quantum dots on solid interfaces, theemission rate is determined by the optical modes that areprojected on the dipole orientation. Therefore, knowl-edge of the projected LDOS is important for controllingspontaneous-emission rates as well as for interpreting thedata from experiments on emitters in photonic metama-terials.A fourth main observation from Fig. 7a and b (Media4 and 5) is that at low frequencies a/λ < .
5, the rela-tive LDOS hardly varies with frequency, which meansthat the mode density is proportional to ω , as in homo-geneous media. Interestingly, there is a clear dependenceof the LDOS on both the position and the dipole orienta-tion even at these low frequencies, i.e. , long wavelengthsrelative to the crystal periodicity. The reason for theseeffects is that the photonic Bloch modes exhibit localvariations of the electric field related to local variationsof the dielectric function in order to satisfy the conti-nuity equations at dielectric boundaries for the parallel E or perpendicular D field, respectively [34, 47]. Conse-quently, the LDOS strongly varies on length scales muchless than the wavelength, i.e., even in electrostatic or ef-fective medium limit. While such behavior may appearsurprising, its origin in electrostatic depolarization ef-fects has been discussed before [47, 48].To gain more insight in the spatial dependence of theLDOS in the inverse opals, we have performed calcula-tions for dipoles positioned along a characteristic axis inthe unit cell. Figure 8 shows the LDOS at four key fre-quencies for dipoles placed on the body diagonal of thecubic unit cell, that is, from (0 , ,
0) in the [1 , ,
1] direc-tion. The LDOS has clear maxima and minima along thediagonal, varying by more than 10 × . Figure 8(a) showsthat for a dipole oriented in the [1 , ,
1] direction (per-pendicular to the dielectric-air interface), the LDOS isstrongly ( > shell( r/a ≈ . .00.51.01.52.02.53.0 e d = [1,1,2]e d = [-1,-1,1]e d = [-1,1,0]Position b R e l . L D O S e d = [-1,1,0]e d = [1,1,1] R e l . L D O S a/ Position c e d = [0,1,0]e d = [1,0,0] R e l . L D O S Position a e d = [0.33,0.13,0]e d = [-0.13,0,0.33]e d = [0,0,1] a/ Position d R e l . L D O S (d)(c) (b)(a) Fig. 9. Relative LDOS in the inverse opal at four differentpositions on the TiO -air interface shown in Fig. 5. Ateach position the LDOS is projected on three mutuallyorthogonal dipole orientations e d . (a, Media 6) point a for e d = [1,0,0] and [0,1,0]. LDOS at e d = [0,0,1] and[0,1,0] is the same. (b, Media 7) point b for e d = [1,1,2],[-1,1,0] and [-1,-1,1]. (c, Media 8) point c for e d = [1,1,1]and [-1,1,0]. LDOS at e d = [-1,-1,2] is equal to that at e d = [-1,1,0]. (d, Media 9) point d for e d = [0.33,0.13,0],[-0.13,0,0.33] and [0,0,1].suppression or enhancement occurs for dipole orienta-tions parallel to the dielectric interface, see Fig. 8(b).The strongly differing LDOS for dipoles near the dielec-tric rationalize the broad distributions of emission ratesthat were recently observed for quantum dots in inverseopals [35], see below. Enhancements and inhibitions alsooccur in the air regions: the LDOS is enhanced at r/a ≈ r/a = 1 / √ ≈ .
57 that lies in the(111) lattice plane in the air region, the LDOS is inhib-ited at all frequencies for the dipole orientations [-1,1,0]and [-1,-1,2] that are perpendicular to the (111) plane(see Fig. 8(b)). Finally, for dipoles parallel to the bodydiagonal (Figure 8(a)) the mode density shows pseudo-oscillatory behavior on length scales much less than thewavelength, e.g. , a period of 0 . a at frequency 0 . a/λ (period corresponds to 1 / th of a wavelength). This ob-servation confirms that photonic crystals are bona fide metamaterials where optical properties strongly vary onlength scales much less than the wavelength.In recent time-resolved experiments [21, 35], emissionfrom quantum dot light sources distributed at the in-ternal TiO -air interfaces of the inverse opals was in-vestigated. To analyze the experimental data, we havecalculated the LDOS at several symmetry-inequivalentpositions on the TiO shells: at points a , b , c and d (see Fig. 5) for three mutually orthogonal orientationsof the emitting dipole, where the first orientation is cho-sen along the vector pointing from (0,0,0) toward thecorresponding point. Figures 9(a) through 9(d) (Media6 through 9) show that at all these positions the LDOS strongly varies with reduced frequency and position, asexpected, and also with orientation of the dipole. Inbroad terms, for a dipole parallel to the interface theLDOS is near 1 at low frequency, increases to a peakenhancement of 2 . × at a/λ = 1 .
22 before stronglyvarying at high frequencies. For reference, the frequency a/λ = 1 .
22 is in the range where a band gap is expectedfor more strongly interacting crystals. The plots also re-veal that for dipoles perpendicular to the TiO -air in-terface, the LDOS is strongly inhibited (more than 10 × )over broad frequency ranges. For instance, Figure 9(a)shows that the emission rate for a dipole perpendicularto the interface is 16-fold inhibited to a level of 0.06, witha maximum of 0.22 (5-fold inhibition) at a/λ = 1 . i.e. , more thantwo octaves in frequency. While the inhibition is notcomplete, the bandwidth is much larger than the maxi-mum bandwidth of 12% for a full 3D bandgap in inverseopals [49] and far exceeds the bandwidth of the 2D gapfor TE modes in membrane photonic crystals, that allowsa seven-fold inhibition [26, 30] over a 30% bandwidth.In the frequency range up to a/λ = 1 .
1, the depen-dence of the LDOS on frequency is quite similar at allthe positions and dipole orientations. This remarkableresult agrees with the experimental observation fromRef. [35]: the complex decay curves of light sources inthe inverse opals are described by one and the same func-tional shape (log-normal) of the decay-rate distributionfor all reduced frequencies studied.For the interpretation of the time-resolved experi-ments on ensembles of emitters in the inverse-opal pho-tonic crystals, the results presented above mean that1. the decay rate of an individual emitter is determinedby its frequency, position and also by the orientationof its transition dipole in the photonic crystal;2. the measured spontaneous-emission decay dependson how the emitters are distributed in the crystal;3. even in the low-frequency regime, ensemblemeasurements will reveal non-exponentional decaycurves;4. the similar shape of the reduced-frequency de-pendence of the LDOS allows modeling of thenon-exponentional decay curves with a single typeof decay-rate distributions. In other words, onedistribution function can be successfully used tomodel the multi-exponentional decay curves meas-ured from crystals with different lattice parameters.
5. LDOS in silicon inverse opals
A complete inhibtion of spontaneous emission may onlybe achieved in photonic crystals with a 3D photonicband gap. Therefore, there has been much effort to fabri-cate inverse opals from silicon rather than titania, sincethe higher index of silicon allows for a photonic band8ig. 10.
Left : Photonic band structure (Media 10) foran inverse opal from silicon ( ǫ = 11.9). Right : The to-tal DOS in the Si inverse opal (Media 11). The DOS isstrongly depleted for frequencies near ΓL and ΓX stop-gaps (grey rectangles). A photonic band gap (grey bar)occurs between bands 8 and 9, as also reflected in thevanishing DOS.Fig. 11. Relative LDOS in a Si inverse opal at: (a, Me-dia 12) r = (0,0,0) [the center of an air sphere, solidcurve], r = (1 , ,
1) [among three air spheres, dash-dotted curve] and r = ( ,0,0) [midway between twospheres along [1,0,0] direction, dotted curve]; (b, Media13) r = (1 , ,
0) [in the window between two air spheres]projected on [1,1,0], [-1,1,0] and [0,0,1] directions shownby solid, dash-dotted and dotted curves, respectively.gap [31, 37]. For the calculation of the LDOS in such Siinverse opals, the dielectric function ǫ ( r ) was modeledsimilarly to that in the inverse opals shown in Fig. 5.From SEM observations we inferred the following struc-tural parameters [17, 18]: the outer radius of the over-lapping dielectric shells with ǫ = 11 . r (where r = 0 . √ a is the air-sphere radius, a is thelattice parameter). The cylindrical windows connectingneighboring air spheres have a radius of 0.2 r . The largerouter radius and the smaller window size compared tothe TiO inverse opals are commensurate with a highervolume fraction of about 23% Si [18].The band structure and total DOS for this system areshown in Fig. 10 (Media 10 and Media 11). The DOS isstrongly depleted in a pseudo-gap between the 2 nd and3 rd bands [36], at frequencies near 0.55. Near frequency0.85, both the band structures and the DOS reveal a3D photonic band gap of relative width ∆ ω/ω ≈ th and 9 th bands [31]. Compared to theTiO structure, both the lowest-order L-gap and the 8 th and 9 th bands are shifted to lower reduced frequencies.This shift is due to the higher effective refractive index ofthe Si inverse opals ( n av = 1 . r =(0,0,0), (1 , ,
1) and ( , , inverse opals,as a result of the larger dielectric contrast that leads tostrongly modified dispersion relations and Bloch modeprofiles. While the maxima up to 3.2 are not much higherthan in TiO inverse opals, the minima in the mode den-sity are reduced and the slopes are steeper. As expected,the LDOS is completely inhibited at all positions in thefrequency range of the band gap. Previously, it has beensuggested that the LDOS could be inhibited at salientpositions in the unit cell for frequencies outside a com-plete gap. While Figure 11(a) (Media 12) reveals thatthe mode density is strongly reduced above the bandgap, e.g., at r = ( , , r = (1 , , a/λ < . − , ,
0) orientation, intermediate forthe (0 , ,
1) orientation, and lowest for the (1 , ,
0) ori-entation. With increasing frequency, the highest LDOSalso changes to other orientations, with the (0 , ,
1) ori-entation having the highest LDOS above the pseudo-gapand even a narrow peak at frequency 0.95. Therefore, ifit is possible to orient a quantum emitter with its dipoleparallel to (0 , ,
6. Conclusions
We have performed intensive calculations of the localdensity of states in TiO and Si inverse opals with exper-imentally relevant structural parameters. Since conflict-ing and incorrect reports have appeared on the LDOS inphotonic crystals, we have set out to validate our methodof choice, i.e., the H-field plane-wave expansion method.This validation relied on comparison to literature results,on the explicit verification of required symmetries thatprevious reports failed to satisfy, and on quantitativeconsiderations of resolution and accuracy. Results for9ach structure are made available for other workers inthe field, both as benchmarks and for comparison withexperimental data. With the help of these computationswe have obtained quantitative insight in the LDOS rel-evant for time-resolved ensemble fluorescence measure-ments on photonic crystals, such as obtained in recentexperimental work [35]. The results of our numerical cal-culations reveal a surprisingly strong dependence of theLDOS on the orientation of the emitting dipoles. Acknowledgments
We thank Dries van Oosten for careful reading of themanuscript. This work is part of the research programof the “Stichting voor Fundamenteel Onderzoek der Ma-terie (FOM),” which is financially supported by the“Nederlandse Organisatie voor Wetenschappelijk Onder-zoek (NWO)”. AFK and WLV were supported by VENIand VICI fellowships funded by NWO. WLV also ac-knowledges funding by STW/NanoNed.
References
1. J. D. Joannopoulos, S. G. Johnson, J.N. Winn and R. D.Meade,
Photonic Crystals: Molding the Flow of Light,2nd edition , (Princeton University Press, 2008).2. C. M. Soukoulis, editor,
Photonic Crystals and Light Lo-calization in the st Century , (Kluwer, Dordrecht 2001).3. V. P. Bykov, “Spontaneous emission in a periodic struc-ture,” Sov. Phys. JETP , 269-273 (1972).4. “Spontaneous emission from a medium with a bandspectrum,” Sov. J. Quant. Electron. , 861-871 (1975).5. E. Yablonovitch, “Inhibited spontaneous emission insolid-state physics and electronics,” Phys. Rev. Lett. ,2059-2062 (1987).6. S. John, “Strong localization of photons in certain dis-ordered dielectric superlattices,” Phys. Rev. Lett. ,2486-2489 (1987).7. O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D.O’Brien, P. D. Dapkus and I. Kim,“Two-dimensionalphotonic band-gap defect mode laser,” Science ,18191821 (1999).8. H.-G. Park, S.-H. Kim, S.-H. Kwon, Y.-G. Ju, J.-K.Yang, J.-H. Baek, S.-B. Kim, Y.-H. Lee, “Electricallydriven single-cell photonic crystal laser,” Science, ,1444-1447 (2004).9. M. Gr¨atzel, “Photoelectrochemical cells,” Nature ,338-344 (2001).10. R. Sprik, B.A. van Tiggelen, and A. Lagendijk, “Opti-cal emission in periodic dielectrics,” Europhys. Lett. ,265-270 (1996).11. N. Vats, S. John and K. Busch, “Theory of fluorescencein photonic crystals,” Phys. Rev. A , 043808 (2002).12. P. Kristensen, A. F. Koenderink, P. Lodahl, B. Trom-borg, and J. Mørk, “Fractional decay of quantum dots inreal photonic crystals,” Opt. Lett. , 1557-1559 (2008).13. K. Busch, S. L¨olkes, R. B. Wehrspohn and H. F¨oll, ed-itors, Photonic Crystals. Advances in Design, Fabrica-tion, and Characterization ,(Wiley-VCH Verlag GmbH& Co., Weinheim 2004). 14. J.E.G.J. Wijnhoven and W. L. Vos, “Preparation of pho-tonic crystals made of air spheres in titania,” Science , 802-804, (1998).15. J.E.G.J. Wijnhoven, L. Bechger and W. L. Vos, “Fab-rication and characterization of large macroporous pho-tonic crystals in titania,” Chem. Mater. , 4486-4499(2001).16. A. A. Zakhidov, R. H. Baughman, Z. Iqbal, C. Cui, I.Khayrullin, S. O. Dantas, J. Marti, and V. G. Ralchenko,“Carbon structures with three-dimensional periodicityat optical wavelengths,” Science , 897-901 (1998).17. A. Blanco, E. Chomski, S. Grabtchak, M. Ibisate, S.John, S. W. Leonard, C. L´opez, F. Meseguer, H. M´ıguez,J. P. Mondia, G. A. Ozin, O. Toader, and H. M. vanDriel, “Large-scale synthesis of a silicon photonic crys-tal with a complete three-dimensional bandgap near 1.5micrometres,” Nature , 437-440 (2000).18. Y. A. Vlasov, X. Z. Bo, J. C. Sturm, and D. J. Norris,“On-chip natural assembly of silicon photonic bandgapcrystals,” Nature , 289-293 (2001).19. A.F. Koenderink, L. Bechger, H. P. Schriemer, A. La-gendijk and W. L. Vos, “Broadband fivefold reductionof vacuum fluctuations probed by dyes in photonic crys-tals,” Phys. Rev. Lett. , 143903 (2002).20. S. Ogawa, M. Imada, S. Yoshimoto, M. Okato and S.Noda, “Control of light emission by 3D photonic crys-tals,” Science , 227-229 (2004).21. P. Lodahl, A. F. van Driel, I. S. Nikolaev, A. Irman, K.Overgaag, D. Vanmaekelbergh, and W. L. Vos, “Control-ling the dynamics of spontaneous emission from quan-tum dots by photonic crystals,” Nature , 654 (2004).22. A. Badolato, K. Hennessy, M. Atat¨ure, J. Dreiser, E. Hu,P. M. Petroff, A. Imamo˘glu, “Deterministic coupling ofsingle quantum dots to single nanocavity modes,” Sci-ence , 1158-1161 (2005).23. A. Kress, F. Hofbauer, N. Reinelt, M. Kaniber, H. J.Krenner, R. Meyer, G. B¨ohm, and J. J. Finley, “Manip-ulation of the spontaneous emission dynamics of quan-tum dots in two-dimensional photonic crystals,” Phys.Rev. B , 241304 (2005).24. M. Fujita, S. Takahashi, Y. Tanaka, T. Asano, S. Noda,“Simultaneous inhibition and redistribution of sponta-neous light emission in photonic crystals,” Science ,1296-1298 (2005).25. D. Englund, D. Fattal, E. Waks, G. Solomon, B. Zhang,T. Nakaoka, Y. Arakawa, Y. Yamamoto, J. Vu˘ckovi´c,“Controlling the spontaneous emission rate of singlequantum dots in a two-Dimensional photonic crystal,”Phys. Rev. Lett. , 013904 (2005).26. B. Julsgaard, J. Johansen, S. Stobbe, T. Stolberg-Rohr,T. S¨unner, M. Kamp, A. Forchel and P. Lodahl, “Decaydynamics of quantum dots influenced by the local den-sity of optical states in two-dimensional photonic crystalmembranes,” Appl. Phys. Lett. , 094102 (2008).27. T. Suzuki and P. K. L. Yu, “Emission power of an elec-tric dipole in the photonic band structure of the fcc lat-tice,” J. Opt. Soc. Am. B , 570-582 (1995).28. R. K. Lee, Y. Xu, A. Yariv, “Modified spontaneous emis-sion from a two-dimensional photonic bandgap crystalslab,” J. Opt. Soc. Am. B , 1438–1442 (2000).29. C. Hermann and O. Hess, “Modified spontaneous emis-sion rate in an inverted-opal structure with complete hotonic bandgap,” J. Opt. Soc. Am. B , 3013–3018(2002).30. A. F. Koenderink, M. Kafesaki, C.M. Soukoulis andV. Sandoghdar, “Spontaneous emission control in two-dimensional photonic crystal membranes,” J. Opt. Soc.Am. B , 1196-1206 (2006).31. K. Busch and S. John, “Photonic band gap formation incertain self-organizing systems”, Phys. Rev. E , 3896-3908 (1998).32. Z.-Y. Li, L.-L. Lin, and Z.-Q. Zhang, “Spontaneous emis-sion from photonic crystals: full vectorial calculations,”Phys. Rev. Lett. , 4341 (2000).33. R. Wang, X.-H. Wang, B.-Y. Gu and G.-Z. Yang, “Localdensity of states in three-dimensional photonic crystals:calculation and enhancement effects,” Phys. Rev. B ,155114 (2003).34. D. P. Fussell, R. C. McPhedran, and C. M. deSterke, “Three-dimensional Green’s tensor, local den-sity of states, and spontaneous emission in finite two-dimensional photonic crystals composed of cylinders,”,Phys. Rev. E , 066608 (2004).35. I. S. Nikolaev, P. Lodahl, A. F. van Driel, A. F. Koen-derink, and W. L. Vos, “Strongly nonexponential time-resolved fluorescence of quantum-dot ensembles in three-dimensional photonic crystals,” Phys. Rev. B , 115302(2007).36. K. M. Ho, C. T. Chan and C. M. Soukoulis, “Existence ofa photonic gap in periodic dielectric structures”, Phys.Rev. Lett. , 3152-3155 (1990).37. H.S. S¨oz¨uer, J.W. Haus, and R. Inguva, “Photonicbands: convergence problems with the plane-wavemethod,” Phys. Rev. B , 13962 (1992).38. J.D. Jackson, Classical Electrodynamics , (John Wiley &Sons, New York 1975).39. N. W. Ashcroft and N. D. Mermin,
Solid State Physics ,(Holt, Rinehart and Winston, New York, 1976).40. L. Li, “Use of Fourier series in the analysis of discon-tinuous periodic structures,” J. Opt. Soc. Am. A ,1870-1876 (1996).41. H. J. Monkhorst and J. D. Pack, “Special points forBrillouin-zone integration,” Phys. Rev. B , 5188-5192(1976).42. G. Gilat, “Analysis of methods for calculating spec-tral properties in solids,” J. Comput. Phys. , 432-465(1972).43. R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACKUsers Guide: Solution of Large-Scale Eigenvalue Prob-lems with Implicitly Restarted Arnoldi Methods , (SIAMPublications, Philadelphia 1998).44. S. G. Johnson and J. D. Joannopoulos, “Block-iterativefrequency-domain methods for Maxwell’s equations in aplanewave basis,” Opt. Express , 173-190 (2001).45. W. L. Vos and H. M. van Driel, “Higher order Braggdiffraction by strongly photonic fcc crystals: onset of aphotonic bandgap,” Phys. Lett. A , 101 (2000).46. P.J. Harding, Photonic crystals modified by opti-cally resonant systems , 1736-1738 (2003). 48. H. Miyazaki and K. Ohtaka, Phys. Rev. B , 6920(1998).49. A.F. Koenderink, Emission and transport of light inphotonic crystals