Periodicity enclosed in boundaries : local density of states in photonic clusters

We have studied the effect of the boundaries on the local density of states (LDOS) in one-, two-, and three-dimensional finite size photonic structures. The LDOS was calculated with the help of the local perturbation method (LPM) and a new LPM-Bloch method using periodicity of the system. The methods are applicable for the clusters made of small (relative to the incident wavelength) particles or for the clusters which can be considered as made of such particles. It was demonstrated that the LPM-Bloch method is an accurate numerical tool for calculation of the LDOS in the finite size photonic structures with weak interference. © 2008 Optical Society of America OCIS codes: (350.4238) Nanophotonics and photonic crystals; (290.5850) Scattering, particles References and links 1. A. A. Asatryan, L. C. Botten, N. A. Nicorovici, R. C. McPhedran, and C. Martijn de Sterke, “Frequency shift of sources embedded in finite two-dimensional photonic clusters,” Waves in Random and Complex Media, 16, 151-165 (2006) 2. A. F. Koenderink, A. Lagendijk, and W. L. Vos, “Optical extinction due to intrinsic structural variations of photonic crystals,” Phys. Rev. B, 72, 153102 (4 pp.) (2005) 3. F. G. Bass, V. D. Freilikher, and V. V. Prosentsov, “Electromagnetic wave scattering from small scatterers of arbitrary shape,” J. of Elm. Wav. and Appl., 14, 269-283 (2000) 4. V. Prosyentsov and A. Lagendijk, “The local density of states in finite size photonic structures, small particles approach,” Photonics and Nanostructures, 5, 189-199 (2007) 5. S. Kole, New methods for the numerical solution of Maxwell’s equations, Ph. D. thesis, University of Groningen (2003) 6. M. Wubs and A. Lagendijk, “Local optical density of states in finite crystals of plane scatterers,” Phys. Rev. E, 65, 046612 (12 pp.) (2002) 7. R. C. McPhedran, L. C. Botten, J. McOrist, A. A. Asatryan, C. M. de Sterke, and N. A. Nicorovici, Density of States Functions for Photonic Crystals, Phys. Rev. E, 69, 016609 (16 pp.) (2004)


Introduction
It is well known in physics that sometimes it takes much less effort to study an ideal infinite system than a finite one.This fact has deep roots and far-reaching consequences.One of them is that infinite-theory formalism is frequently used to study real systems which are finite ones.
Consider for example an infinite photonic crystal containing point sources placed periodically with the same period as the system.Obviously, the field generated by these sources should be also periodic and it is sufficient in this case to consider the fields in only one unit cell.The same should be true for the identical system of finite size when boundary effects can be neglected.The important question is: are the boundary effects really weak in finite photonic crystals?One can not answer on this question in advance.There are photonic systems where the boundary effects are important for any size of the system because of the strong interference between boundary and volume.Recently, the impact of the boundary on the local density of states (LDOS) and the frequency (Lamb) shift was investigated in the work [1] for two-dimensional photonic crystals.It was shown that the LDOS and the Lamb shift are sensitive functions of the shape of the cluster.In reality, the boundary effects can be screened due to multiple imperfections and loses existing copiously in photonic crystals [2].Without the losses, when scattering by particles is strong, the boundaries should affect the fields inside the cluster.This is not the case for an ideal (without defects) infinite system where the particles interact but never scatter as a boundary does.
In this paper we calculate the local density of states (LDOS) in a photonic cluster by using two methods.The first method is the classical local perturbation method (see for example works [3] and [4]) ideally suited for finite structures made of small particles.The second one is a new method taking into account the periodicity of the scatterers in the cluster.In addition it also uses elements of the local perturbation method (LPM).In the following we will call the second method the LPM-Bloch method.The method is based on the assumption that the scatterers (or particles forming the scatterers) are small compared to the incident wavelength and that fields inside a cluster are almost periodic, implicitly implying what the boundary effects are weak.The advantage of the new method is that a much larger cluster can be calculated because the computer memory demand is much less then for the LPM.
In this paper we will compare the two methods for the case when the results of the LPM are exact (with the high numerical accuracy).For this purpose we study LDOS in the clusters made of particles which are small compared to the incident wavelength (see also [4] for comparison between LPM and exact results).We come to the surprising conclusion that the results of the LPM and the LPM-Bloch method do not converge when the size of the cluster grows.
It is interesting to note that in the work [5] the LDOS was studied for the finite photonic structures with the help of the periodic boundary conditions (PBC).It was shown that in band gap region the LDOS tends to the LDOS of an infinite system when the size of the cluster increases.From the physical point of view, the application of the PBC means that one considers an infinite system where the calculation accuracy grows with the density of k vectors (implicitly taken into account) and the density of the vectors is proportional to the cluster's size.From a calculation with the PBC it is difficult to extract the properties of a realistic cluster.

Theory of the LPM-Bloch method
The local density of states ρ in a dielectric medium is proportional to the imaginary part of the trace of the Green's tensor in polarization space and defined as (see, for example [6] and interesting discussions therein about the definition of the LDOS) where ω is the circular frequency, c is the velocity of light in vacuum, ε 0 is the permittivity of the host medium, and G(r, r ) is the Green's tensor at observation point r with source placed at point r .Im and Tr denote imaginary part and trace respectively.Note that the formula (1) holds in any dielectric medium without gain.The expression (1) shows that the LDOS can be easily computed when the Green's tensor G(r, r ) is known.In the following we will calculate the Green's tensor by using the LPM and by using the periodicity of the system.The equation for the quasiperiodic Green's tensor G(r, r , k) in the finite periodic medium with periodically placed sources is defined by the following equation (see for example [7] where an infinite system was considered) Here ⊗ defines tensor product, ε(r) is the permittivity of the medium containing scatterers, I is the 3 × 3 unitary tensor in polarization space, and δ is the Dirac delta function.The vector r n ≡ ∑ 3 i=1 m i a i , where m i are integers and a i are the basis vectors of the system having a Brillouin zone with volume V BZ .Note that the vector k defines the quasiperiodicity or Bloch property of the tensor G(r, r , k) such that in the case of the infinite system or a finite system with periodic boundary conditions and the Green's tensor G(r, r ) is found from the formula Thus, solving the equation (2) relative G(r, r , k) one can find the Green's tensor G(r, r ) with the help of (4).We assume that the condition (3) is valid for finite size photonic structures.In this case the tensor G(r, r , k) can be found with the help of the LPM as it is shown in Appendix.
Below we present the final result for the Green's tensor (the details of the solution of the Eq. (2) can be found in Appendix) where , and f n (q) is the Fourier transform of the function describing the shape of the scatterer (see also Appendix) and is the quasiperiodic Green's tensor of the free space.The formula ( 5) is the main result of this section.It explicitly shows that the total Green's tensor is a sum of the free space term G 0 and term due to the scattering by particles.The later is proportional to the contrast of the particles ε n − ε 0 and depends on the particle's shape f n .The interference between the scatterers as well as their resonance properties are taken into account by only one term G(0, r , k) that should be calculated by solving the system of 9 linear equations obtained by substituting r by r n into (5).In contrast, the LPM requires to solve the system of 9 × N 1 × N 2 × N 3 linear equations that is computer memory and time consuming.The formula (5) can be essentially simplified for identical ( case the second integral in ( 5) can be explicitly calculated with the help of the of method of residues.Finally we obtain where is the distance between the observation point and n-th scatterer.Note that f n (k 0 ) V n /8π 3 , where V n is volume of the n-th scatterer.Formula (7) gives the Green's tensor of a medium periodically filled with N identical spherically symmetric particles which are small compared to the incident wavelength.Note that, when the scatterers are not sufficiently small (k 0 L ≥ 0.2, for example), they should be subdivided into smaller particles with size L sub such that k 0 L sub ≤ 0.1.It is interesting to note that the expression (7) describes the Green's tensor of infinite photonic crystal when N → ∞.Moreover, the formula explicitly shows that even small losses in the host medium (small positive imaginary part of k 0 ) limit the number of interacting particles in infinite photonic crystal.
It should be emphasized here that the main assumption of the LPM-Bloch method is the existence of the quasiperiodicity (expression (3)) for a finite size photonic system.Obviously, the assumption is not correct for the lossless systems where light freely propagates from one boundary to another.Weak interference between scatterers (due to small contrast or large distance, for example) and small losses allow to treat finite system as infinite one because of the effective absence of the boundaries.In this case the expression (3) should be correct.When there is a deviation from the quasiperiodicity, the error of the LPM-Bloch method increases with the number of particles (and the size of the cluster as well) as can be seen from the formula (7).Additional discussion about influence of the quasiperiodicity on the error of the LPM-Bloch method can be found at the end of the Sec.

Numerical examples
In this section we present the numerically calculated LDOS for photonic clusters by using the LPM and the LPM-Bloch formalism.The error of the LPM-Bloch method, on one hand, increases with number of the scatterers and their contrast.On the other hand, the dimensionality of the system and the distance between the particles effectively decrease the influence of the boundaries (as well as interference between the scatterers) and as a result, the error of the LPM-Bloch method drops.In the case of the vector wave situation becomes more complicated due to the strong near field interference between the particles.When the distance between the particles grows, the interference sharply goes down and the accuracy of the LPM-Bloch calculation effectively increases.

Scalar case: 1-3D systems
Consider first the scalar wave propagation in simplest one-dimensional (1D) system formed by periodically placed layers.The Fig. 1(a) and 1(b) show the LDOS calculated for 1D systems with the help of the LPM and the LPM-Bloch formalism.The Fig. 1(a) shows the normalized LDOS versus size parameter ωd/2πc for the system consisting of 21 layers at observation point x = 0.1 μm.The layers are centered at x = 0 and have width L = 0.0125 μm and permittivity ε sc = 1.5.The period of the system is d = 1 μm and the permittivity of the host medium is ε 0 = 1.The red curve shows the LPM result and the blue dashed one represents the LPM-Bloch result.One can see the very good agreement between the results until ωd/2πc = 1.5.When we increase the permittivity, the width, or the number of layers the accuracy of the LPM-Bloch method decreases.As an example, the Fig. 1(b) shows the LDOS for the same system consisting of layers with ε sc = 3.Our calculations show that period does not affect the accuracy.This is due to infinite interference length between layers in 1D system.For 2D and 3D systems, in contrast, the interference between particles decreases with distance between them.The decreasing of the accuracy of the LPM-Bloch method with increasing the contrast or the number of layers 0.02 0.04 0.06 0.08 0.1 0.5 Fig. 3.The normalized scalar LDOS of the finite 3D system versus size parameter ωd/2πc at the point x = 0.01 μm, y = z = 0.The system consists of 13 × 13 × 13 (Fig. 3(a)) and 21 × 21 × 21 (Fig. 3(b)) spheres with radius L = 0.01 μm and permittivity ε sc = 9.The period of the system is d = 0.05 μm and the permittivity of the host medium is ε 0 = 1.The particles are arranged into simple cubic lattice centered at r = 0.The red curve shows the LPM result ("exact result").The blue dashed and solid lines represent the LPM-Bloch results calculated by using 11 and 81 subdivisions in k-space (along each axis) for numerical integration.
indicates that lossless finite one-dimensional system has no LDOS quasiperiodicity.
Another interesting example is the two-dimensional (2D) system consisting of the small circles arranged into square lattice.The Fig. 2(a) and 2(b) show the LDOS calculated for the 2D systems with the help of the LPM and the LPM-Bloch formalism.The Fig. 2(a) shows the normalized LDOS versus size parameter ωd/2πc for the system consisting of 13 × 13 circles at observation point x = 0.01 μm, y = 0.The circles are centered at x = y = 0 and have radius L = 0.01 μm and permittivity ε sc = 9.The period of the system is d = 0.05 μm and the permittivity of the host medium is ε 0 = 1.The red curve shows calculations made with the help of the LPM and the blue ones (dashed and solid) represent the LPM-Bloch result.One can see relatively good agreement between the results despite the high index contrast.The Fig. 2(b) shows the LDOS for the same system consisting of 21 × 21 circles.One can see that the accuracy of the calculations deteriorates with the number of particles (as in 1D case), however less pronounced.The reason is that 2D systems are more open than 1D ones and the interference in 2D systems decreases as r −1/2 while in 1D systems it oscillates.
The last example is the three-dimensional (3D) system consisting of the small spheres arranged into cubic lattice.The Fig. 3(a) and 3(b) show the LDOS calculated for the 3D systems with the help of the LPM and the LPM-Bloch formalism.The Fig. 3(a) shows the normalized LDOS versus size parameter ωd/2πc for the system consisting of 13 × 13 × 13 spheres at observation point x = 0.01 μm, y = z = 0.The spheres are centered at r = 0 and have radius L = 0.01 μm and permittivity ε sc = 9.The period of the system is d = 0.05 μm and the permittivity of the host medium is ε 0 = 1.The red curve shows the result calculated with the help of the LPM and the blue ones (dashed and solid) represent the LPM-Bloch result.One can see relatively good agreement between the results.The Fig. 3(b) shows the LDOS for the same system consisting of 21 × 21 × 21 spheres.One can see that there is no obvious loss of accuracy due to the number of the particles (for this period of the system).We would like to emphasize that we Fig. 4. The normalized vector LDOS of the finite 3D system versus size parameter ωd/2πc at the point x = 0.01 μm, y = z = 0.The system consists of 9 × 9 × 9 spheres with radius L = 0.01 μm and permittivity ε sc = 9.The periods of the systems are d = 0.025 μm (Fig. 4(a)) and d = 0.05 μm (Fig. 4(b)) and the permittivity of the host medium is ε 0 = 1.The particles are arranged into simple cubic lattice centered at the point r = 0.The red curve shows the LPM result ("exact result").The blue dashed and solid lines represent the LPM-Bloch results calculated by using 11 and 81 subdivisions in k-space (along each axis) for numerical integration.used typical values for the parameters and other values bring us to the same conclusions.

Vector case: 3D systems
The Fig. 4(a) and 4(b) show the vector LDOS calculated for 3D system with the help of the LPM and the LPM-Bloch formalism.It show the normalized LDOS versus size parameter ωd/2πc for the system consisting of 9 × 9 × 9 spheres at observation point x = 0.01 μm, y = z = 0.The spheres are centered at r = 0 and have radius L = 0.01 μm and permittivity ε sc = 9.The periods of the systems are d = 0.025 μm and d = 0.05 μm respectively.The permittivity of the host medium is ε 0 = 1.The red curves show the result calculated with the help of the LPM and the blue ones (dashed and solid) represent the LPM-Bloch result.One can see relatively good agreement between the results.However, there is distinctive difference between the LPM and the LPM-Bloch curves.The difference does not vanish with increasing of the number of subdivisions in k-space and it becomes practically constant (offset to some extent) for low frequencies.When the period grows the difference decreases considerably.The Fig. 5 shows the relative error Δ defined as versus normalized period of the cluster d/2L at wavelength λ = 14.5 μm (static regime).Solid line shows the relative error Δ while the dashed one shows cubically decreasing function ∼ (d/2L) −3 for comparison.The Fig. 5 shows that the results predicted by the LPM and the LPM-Bloch formalism become practically identical (Δ = 0.015) for values d/2L 1.5.This indicates that the photonic cluster made of spheres can be treated as infinite one (at least for the LDOS calculations) when d 3L.
The best way to understand the nature of the difference at low frequencies is to compare the Green's tensors calculated with the help of the LPM and with the LPM-Bloch method.For our purposes it will be sufficient to consider the cluster made of identical spheres.In this case the Green's tensor calculated with the help of the LPM takes the following form (the result can be easily obtained form the formula (11) in work [4]) where G(r n , r ) is the Green's tensor inside n-th scatterer.The Green's tensor calculated with the help of the LPM-Bloch method is given by the formula (7).Comparing the formulas (7) and (8) we can see that the expressions will coincide when G(r n , r ) = BZ G(0,r ,k) V BZ e ik•r n dk, i.e. when the condition of quasiperiodicity (3) will be satisfied.For the static regime when the incident wavelength is much larger then the cluster's size (k 0 N 1/3 d 1) the tensor G(0, r , k) = G(r 1 , r , k) , satisfies the following equation while the tensor G(0, r ) can be obtained from the following system of equations where one should use G(0, r ) = G(r 1 , r ) . Note that the sums in the expressions (9) and (10) are due to interference between the particles forming the cluster and in the absence of the interference the expressions are identical (after the integration in k space required for the formula (9)).The solution of the equation (10) can be found by substituting of the tensors G(r n , r ) in this equation by the expression (10) for the points r n .For the tensor G(0, r ) we have where we wrote down only three first terms.Comparing the expression (11) with the expression (9) (after integration in k space) one can see that they differ by the term proportional to α 2 /r 6 jn .The factor α 2 r −6 jn is of the order of 1 for smallest distance between the particles (d = 2L) when the error is maximal.When the period between the particles grows to d = 3L, the factor α 2 r −6 jn decreases in 10 times that roughly corresponds to the situation presented in the Fig. 5.This analysis shows that the LPM-Bloch formalism does not take into account the second order interference effects which are very strong for short periods and exist only in finite size structures (for infinite structures the formalism is exact).

Conclusions
We have introduced a new approximate LPM-Bloch method calculating fields in the photonic clusters with weak interference between particles.We have demonstrated that the LPM-Bloch method is an accurate numerical tool for the calculation of the LDOS in such photonic structures.
The boundaries affect the LDOS inside a photonic cluster when interference between scatterers is strong.These boundary effects should be taken into account when the LDOS is calculated in low loss and high contrast photonic clusters.

Appendix: Computation of the Green's tensor by using the periodicity and the local perturbation method
Before solving the equation (2) we note that for simple lattices permittivity of the medium can be represented in the following form where ε n is the permittivity of the n-th scatterer, N is the total number of the scatterers, r n = ∑ 3 i=1 m i a i is the radius vector of the n-th particle, and m i are integers.The function f n (r − r n ) is the unit step function describing the shape of the n-th particle, such that Note, that the complex lattices containing more than one particle in unit cell can be described the similar way.Taking into account the formula (12) we rewrite the equation ( 2) for the quasiperiodic Green's tensor in the following form where k 0 = ω c √ ε 0 = 2π/λ 0 is the wave number in the host medium.
An approximate solution of the equation ( 14) can be found with the help of the local perturbation method (see for example work [3]) when the characteristic size (L) of the scatterer is small compared to the incident wavelength, i.e. when k 0 L 1.The key hypothesis of the local perturbation method (LPM) is the assumption that the solution G(r, r , k) is a relatively slow function of coordinates, so that .The equation ( 16) can better be solved by using the spatial Fourier transformation on the coordinate r.The Fourier transform G(q, r ) of the quasiperiodic Green's tensor takes the form , where

Fig. 1 .
Fig.1.The normalized scalar LDOS of the finite 1D periodic system versus size parameter ωd/2πc for observation point x = 0.1 μm.The system consists of 21 layers centered at x = 0.The permittivities of the layers are ε sc = 1.5 (Fig.1(a)) and ε sc = 3 (Fig.1(b)) and their width is L = 0.0125 μm.The period of the system is d = 1 μm and the permittivity of the host medium is ε 0 = 1.The red curve shows the LPM result ("exact result") and the blue dashed one represents the LPM-Bloch result calculated by using 201 subdivisions in k-space for numerical integration.

Fig. 2 .
Fig.2.The normalized scalar LDOS of the finite 2D system versus size parameter ωd/2πc at the point x = 0.01 μm, y = 0.The system consists of 13 × 13 (Fig.2(a)) and 21 × 21 (Fig.2(b)) circles with radius L = 0.01 μm and permittivity ε sc = 9.The period of the system is d = 0.05 μm and the permittivity of the host medium is ε 0 = 1.The scatterers are arranged into square lattice centered at x = y = 0.The red curve shows the LPM result ("exact result").The blue dashed and solid lines represent the LPM-Bloch results calculated by using 11 and 81 subdivisions in k-space (along each axis) for numerical integration.

ΔFig. 5 .
Fig. 5.The relative error Δ versus normalized period of the cluster d/2L at wavelength λ = 14.5 μm.The solid line shows Δ while the dashed one shows cubically decreasing function ∼ (d/2L) −3 for comparison.The LDOS was calculated at the point x = 0.01 μm, y = z = 0 for the system consisting of 9 × 9 × 9 spheres with radius L = 0.01 μm and permittivity ε sc = 9.The permittivity of the host medium is ε 0 = 1.The particles are arranged into simple cubic lattice centered at the point r = 0.