Table Of ContentA flexible Bloch mode method for computing complex band structures and
impedances of two-dimensional photonic crystals
Felix J. Lawrence,1,a) Lindsay C. Botten,2 Kokou B. Dossou,2 R. C. McPhedran,1 and C. Martijn de Sterke1
1)CUDOS and Institute of Photonics and Optical Science (IPOS), School of Physics, University of Sydney,
NSW 2006, AUSTRALIA
2)CUDOS and Department of Mathematical Sciences, University of Technology, Sydney, NSW 2007,
AUSTRALIA
(Dated: 17 January 2012)
We present a flexible method that can calculate Bloch modes, complex band structures, and impedances of
two-dimensional photonic crystals from scattering data produced by widely available numerical tools. The
method generalizes previous work which relied on specialized multipole and FEM techniques underpinning
transfermatrixmethods. Wedescribethenumericaltechniqueformodeextraction,andapplyittocalculatea
complexbandstructureandtodesigntwophotoniccrystalantireflectioncoatings. Wedothisforfrequencies
at which other methods fail, but which nevertheless are of significant practical interest.
2 I. INTRODUCTION the first Wood anomaly, which ranges from a /λ = 1/n
x
1 fornormallyincidentlighttoa /λ=1/2nforlightatthe
x
0 When modeling photonic crystals (PCs), it is impor- Brillouin-zone edge, where ax is the length of the lattice
2
tant to consider all the relevant Bloch modes. Light at vector parallel to the interface, λ is the free space wave-
n a fixed frequency, polarization, and incident angle exists length and n is the PC’s background index. Above this
a in a PC as a superposition of a set of propagating and frequency, generally several Bloch modes must be simul-
J
evanescent Bloch modes, the PC’s eigenstates. At low taneously considered in each PC, regardless of whether
6 frequencies, only one mode generally needs to be con- these modes are propagating or evanescent. Reflection
1
sidered. For light at frequencies above the first Wood at a PC/PC interface is well-described by a matrix that
] anomaly1,eachrowofholesinthePCdiffractslightinto maps incident modes to reflected modes, as we have
s several propagating orders, so the PC may support mul- shown previously4,12. In our experience, the minimum
c
i tiple propagating Bloch modes. At the PC’s front and acceptable dimension of this reflection matrix, as argued
pt back interfaces, some of its modes couple via reflection, in Sec. IIA, is usually
o affecting the overall reflection and transmission through
(cid:106)a (cid:107) (cid:106)a (cid:107)
. the PC, so it is important to model all relevant modes. M = x(1+sinθ ) + x(1−sinθ ) +1, (1)
s min nλ i nλ i
c It is often important to include evanescent modes2. If
si the PC is not long—for example, if it is a layer in a thin where θ is the incident angle from a uniform dielectric
i
y antireflection coating—then evanescent modes can play with the PC’s background index, and (cid:98)x(cid:99) denotes the
h a role in energy transport3. Evanescent modes can also floor of x. We have previously achieved accurate results
p
play a role in field matching across an interface between modeling PC stacks using impedance matrices of this di-
[
PCs4 or PC waveguides5. The propagative qualities of mension and higher4,12,13.
2 an evanescent mode are well-represented by its complex AnumberofmethodsforfindingmultipleBlochmodes
v band structure6, which augments the traditional band and complex band structures have been demonstrated.
5 structure,conveyinginformationabouttherateatwhich Transfer-matrix14 and scattering-matrix15 based meth-
3
4 the mode accumulates phase together with information ods were developed to derive a PC’s Bloch modes from
3 about the mode’s decay rate. the properties of a single grating layer. The plane wave
. There have been a number of studies seeking to de- expansion method has also been extended to include
1
rive impedance-like quantities to characterize reflection evanescent modes16. Finally, Ha et al. presented a
1
1 atPCinterfacesbyascalar7,8. Furthermore,anumberof methodforextractingBlochmodesfromtheoutputofan
1 studieshaveadaptedmetamaterialparameterextraction EM solver17, or even near-field measurements18,19. We
: techniques9tophotoniccrystals,andusedthemtodesign improve the accuracy, stability and efficiency of Ha et
iv antireflection coatings10,11. However, since these tech- al.’s method and extend it to calculate PC impedances
X niques characterize reflection and transmission by a sin- for two-dimensional (2D) PCs, which can be used to cal-
r gle complex number each, they cannot handle problems culatereflectionandtransmissionatinterfaces4,12. These
a involvingmultiplemodes,whereeverymodereflectsinto PCimpedancesandthereflectionandtransmissionoper-
every other mode. Scalar-based methods generally give ators are represented by matrices; our method supports
manifestlyincorrectresultsforlightatfrequenciesabove thepresenceandinteractionofmultipleBlochmodesand
so it can work well both above and below the first Wood
anomaly.
Wehavemadesoftwareavailablethatusesthemethod
a)[email protected] described in this paper to calculate PCs’ Bloch modes,
2
complex band structures, and impedances. The soft- a isthelengthofthelatticevectorparalleltothex-axis.
x
ware,calledBlochCode,canthenusethesecomplexband The wavevector component in the direction perpendicu-
(cid:113)
structures and impedances to calculate reflection and lar to the grating is k(p) = k2−k(p)2 where k is the
transmissionmatricesandcoefficientsforarbitrarystacks y x
wavenumber in the medium. Evanescent grating orders
ofPCs. BlochCodeisopen-sourceandisavailableonthe
internet20. haveimaginaryky(p),soforagivenkandkx(p),thenumber
In Sec. II, we present our method for finding Bloch of propagating grating orders is the number of solutions
modesfromtheelectricfieldE andthemagneticfieldH to Eq. (2) with real ky(p), or Mmin in Eq. (1). In our
in a PC structure. Sec. IIA recaps some useful results experience, M also provides an upper bound on the
min
fromourpreviouswork12 andprovidessomebackground number of propagating Bloch modes, and at non-normal
theory. Sec. IIB details our improvements to Ha et al.’s incidenceisalowerboundonthenumberofBlochmodes
method17 of finding Bloch factors and modal fields, and requiredtomodelaPCaccurately. Atnormalincidence,
Sec. III outlines our procedure for successfully applying symmetry allows odd modes to be ignored, so in this
thismethodtominimizetheresidualderivedinSec.IIB. casegoodresultsmaybeobtainedwithfewerthanM
min
Sec.IICexplainshowwecalculatePCimpedancematri- modes—seeSec.IVC.UsingBlochmodesfoundfromac-
cesfromthemodalfields. InSec.IVweapplyourmethod curatemultipoleandFEMtransfermatrixmethods23,24,
to demonstrate its utility. In Sec. IVA we calculate the we have consistently had success modeling PCs with no
complex band structure for light normally incident on a more than M +2 Bloch modes.
min
triangular lattice PC. In Sec. IVB we reproduce the de-
sign process of a known antireflection coating for a PC, Bloch’stheoremrelatestheelectricandmagneticfields
at a frequency and incident angle for which it is critical associated with each mode at equivalent points in differ-
to include at least two Bloch modes in the calculations. ent unit cells of a PC. The ratio of each mode’s field
Finally, in Sec. IVC we use our method to design an at points separated by the lattice vector e = (a ,0)
1 x
all-polarization antireflection coating for a square lattice is eikxax. For the PC’s other lattice vector e2, this ra-
self-collimating PC, at a high frequency where a scalar tio is different for each mode and is the mode’s Bloch
method cannot find a coating for the PC21. factor, denoted by µ. Calculating µ for each mode is
the goal of Sec. IIB. For square and rectangular lattices,
e2 = (0,ay) and µ = eikyay, where ky is the y compo-
II. THEORY nent of the mode’s wavevector. For triangular lattices,
thelatticevectore is(a /2,a )andsotheBlochfactor
2 x y
Our method uses a two-step process to extract a PC’s may be written µ=ei(kxax/2+kyay).
modes and impedance from the field in a finite length of
the PC. The PC is assumed to be two-dimensional, loss-
Bloch modes come in forward/backward pairs. Popov
less, and to have relative permeability µ = 1. Like Ha
r et al. provide a useful discussion of symmetry
et al.’s method17, we could use data generated by FEM properties25. We assume mirror symmetry in each unit
or FDTD simulations, or even experimentally measured
cell,whichmeansthateachbackwardmode’sfieldprofile
by a near-field probe such as a SNOM19, although the
in a unit cell is the reflection on the x-axis of its forward
impedance part of our method is not valid for SNOM
partner’s. TheBlochfactorsofapairarerelatedbecause
data,whichisderivedfroma3Dobject. First,theBloch
of this: for square and rectangular lattices, µ = 1/µ ,
b f
factors and the Bloch modal fields are found (Sec. IIB),
where µ and µ are respectively the Bloch factors of
f b
then these modes are analyzed to calculate the PC’s
the forward and backward modes. For triangular-like
impedance (Sec. IIC).
lattices, the symmetry is more complicated since the re-
flection of e is not −e , the translation corresponding
2 2
to the field ratio 1/µ , but (a /2,−a ); these vectors
f x y
A. Background Theory
differ by −e . Accounting for this discrepancy, we find
1
µb =e−ikxax/µf for triangular lattices.
Two-dimensional PCs in the x−y plane may be de-
scribed as a stack of gratings parallel to the x axis22,
A PC’s impedance is defined in terms of two matrices,
each of which diffracts incident light into an infinite set
E and H12. For E = E polarized light, each matrix
z
of grating orders. At the edge of each unit cell, the PC’s
mapsavectorofforwardBlochmodeamplitudesc toa
+
Bloch modes may be written as a superposition of the
vectoroftheE orH fieldsassociatedwitheachgrating
underlying grating orders15. Their directions are given z x
diffractionorder. Specifically,E ,the(p,m)thelement
p,m
by the grating equation
ofE,istheE fieldofnormalizedmodemduetoforward
z
2πp 2πp andbackwardplanewavesingratingorderp,atthecen-
k(p) =k + =ksinθ + , (2)
x x a i a tre(x=0)ofaunitcell’sedge. Thus,forasetofforward
x x
propagating/decaying Bloch modes c , the field compo-
+
wherekx isthexcomponentoftheincidentplanewave’s nents along the edge of the unit cell, i.e., the quantities
wavevector, k(p) is that of the pth diffraction order, and that are continuous across an interface between PCs or
x
3
dielectrics, are
Ez(x)=(cid:88)Epc+eikx(p)x, Hx(x)=(cid:88)Hpc+eikx(p)x,
p p l=0 1 2 3 4
(3)
whereE andH aretherowsofEandHcorresponding
p p
to grating order p. In the H = H polarization, E and
z
H map to E and H fields, and these quantities replace
x z
Ez and Hx in Eq. (3). l=0
Previously12, we defined PC impedances in terms of y 1
these matrices. For E polarized light, the impedance of
z 2
a PC is
x 3
Z =H T(I+Q)E+E T(I−Q)H, (4) 4
0 0
FIG. 1: Schematic of L=5 PC structures for a square
and for H polarized light it is
z and a triangular PC lattice. The squares with solid
(cid:16) (cid:17) edges are the unit cells used by our method. For the
Z =− H0T(I−Q)E+E0T(I+Q)H , (5) triangular lattice PC, the field in the solid-edge unit
cells are calculated from the unit cells of the simulated
where E and H are calculated for the PC, and E and structure (dashed edges) using Bloch’s theorem, with
0
H0 are calculated for a reference material, usually free the ratio eikxax between adjacent cells’ fields.
space. Q is a diagonal matrix that takes into account
thehalf-periodshiftofgratingsintriangularlatticePCs:
forsquarelatticesQ=I, andfortriangularlatticesQ= cells separated by a lattice vector (solid edges); we ap-
diag((−1)p), where p is the grating order. ply Bloch’s theorem with integer multiples of the lattice
Given impedances Z and Z for two PCs, it is sim- vector (a ,0).
1 2 x
ple to calculate the reflection and transmission matrices We seek to write these electric field components as a
across their interface12: superposition of forward and backward Bloch modes. So
we want to express every U (r), i.e., the E or E field
(cid:96) z x
T12 =(AT12A12+I)−12AT12, (6a) component for sampled point r in unit cell (cid:96), as
R =(A AT +I)−1(A AT −I), (6b)
12 12 12 12 12 U (r)=(cid:88)µ(cid:96) A (r)+(cid:88)(1/µ L−1−(cid:96))A (r)+w((cid:96),r),
(cid:96) m m m(cid:48) m(cid:48)
where A12 =Z1−1Z2. m m(cid:48)
(7)
where A (r) and µ are respectively the modal field
m m
B. Finding modes and the Bloch factor of forward mode m; m(cid:48) denotes
backward modes, and w((cid:96),r) is the residual error. More
specifically,forforwardmodes,A (r)isthefieldcompo-
OurmethodoffindingtheBlochmodesandBlochfac- m
nentofmodematpointrofthefirstunitcell,(cid:96)=0. The
tors is based on the method presented by Ha et al.17, al-
Bloch factor µ is the ratio of the field in cells (cid:96)+1 and
though our method offers some significant improvements m
(cid:96), so µ(cid:96) A (r) is the field component of forward mode
in accuracy and efficiency. We take field data for several m m
m at point r of unit cell (cid:96). To avoid ill-conditioning, the
unitcellsofaPC,andtrytowriteitasasuperpositionof
field A (r) at point r of each backward mode m(cid:48) is de-
Blochmodes,thusfindingthemodalfieldsandBlochfac- m(cid:48)
finedinthelastunitcell,(cid:96)=L−1. Thismeansthatthe
tors. The final steps of our mode-finding method impose
coefficients of A (r) and A (r) in Eq. (7) have moduli
symmetry relationships between forward and backward m m(cid:48)
no greater than 1. As noted in Sec. IIA, the Bloch fac-
modal fields, increasing accuracy by almost halving the
tor µ of each backward mode is related to that of its
numberofunknownsintheproblem. Wenowoutlineour m(cid:48)
forward partner; we enforce this relationship in practice,
method.
thereby halving the number of Bloch factors that must
In an EM solver, we simulate a section of 2D PC
be found.
withBloch-Floquetperiodicboundaryconditionsontwo
Equation(7)forall(cid:96)andallsampledrmaybewritten
boundaries,anduniformdielectricontheothers(Fig.1).
in matrix form as:
WesampletheE orE (dependingonpolarization)field
z x
component at many (Np) points in unit cell (cid:96) = 0, and U=CA+W, (8)
then at the equivalent points in each of the other unit
cells. If desired, E , H , H , or H may be used in whereUcontainstheE orE fieldcomponentsfromthe
y x y z z x
place of or in addition to E and E . For triangular EM solver, A is a matrix of modal fields, C is a matrix
z x
lattice PCs, we use the field in the simulated unit cells constructed from Bloch factors, and W is a matrix of
(dashed edges in Fig. 1) to calculate the field in the unit residuals w((cid:96),r) that must be minimized. U is a L×
4
N matrix: the field in its (cid:96)th row and rth column is improve the quality of our results by enforcing this rela-
p
U = U (r), the field component at point r in unit cell tionship in the minimization process.
(cid:96),r (cid:96)
(cid:96). Similarly, A is a M ×N matrix; the field in its mth We commence by partitioning the forward (f) and
p
rowandrthcolumnisA =A (r),thefieldofmodem backward (b) modes, and the points in the left (L;
m,r m
atpointrincell(cid:96)=0forforwardmodes,orcell(cid:96)=L−1 y ≤ a /2) and right (R; y ≥ a /2) halves of the unit
y y
forbackwardmodes. CisaL×M matrix. Foraforward cell:
mode m, the ((cid:96),m)th element of C is µ (cid:96), and for a
m
backwardmodem(cid:48),the((cid:96),m(cid:48))thelementis1/µ L−1−(cid:96). U=(UL,UR), C=(Cf,Cb), (12a)
m(cid:48)
Ifmultiplefieldcomponents(e.g. E , H andH )areto
z x y
be used to find the modes, then the additional data can (cid:18) (cid:19)
A A
be added as extra columns in U. A= L,f R,f . (12b)
A A
L,b R,b
We start the optimization process knowing U, and
with information about the structure of C, and no di- After normalization, the field of a backward mode is the
rect information about A. In our method, we first find field of its forward partner reflected about the x-axis,
the Bloch factors that determine C, a relatively difficult thus
problem. OnceCisknown,solvingEq.(8)forthemodal
fields A becomes a pure least-squares problem that can (A ,A )=(cid:0)γA P,γA P−1(cid:1), (13)
L,b R,b R,f L,f
be solved accurately and efficiently using standard tech-
where P is the permutation matrix that maps points
niques.
(x,a −y) to (x,y), and γ is a normalizing diagonal ma-
To find the modes, we seek to minimize the difference y
trix whose elements are the ratio of backward and for-
between the observed field U and the superposition of
ward mode amplitudes. The columns of A and A ,
Bloch mode fields CA. That is, we seek to minimize R,f R,b
||W||2 in Eq. (8), the sum of squared moduli of the el- corresponding to points in the right half of the unit cell,
F
can easily be ordered so that P = I; from now on we
ements of W. Constraining the problem by dividing by
thesquaredFrobeniusnorm||U||2 ofU,thequantitywe assume this ordering. Eq. (8) can now be written with
F
roughly half as many unknowns,
minimize is
(cid:18) (cid:19)
w2 = ||U|−|UC||2A||2F, (9) (UL,UR)=(Cf,Cbγ) AARL,,ff AARL,,ff +W. (14)
F
where w2 = ||W||2/||U||2. First we eliminate A from Cbγ represents each backward mode’s amplitude in each
F F cell, relative to that of the corresponding forward mode
Eq. (9) in order to find C with a numerical minimizer.
in cell 0.
We use an alternative representation of the Frobenius
(cid:112) The constraints on A (Eq. (13)) mean that Eq. (14)
norm, ||U|| = tr(UHU), to write
F
doesnothavealeast-squaresform, somaynotbeimme-
tr((UH −AHCH)(U−CA)) diatelysimplifiedinthewaythatEq.(9)ledtoEq.(11).
w2 = . (10) TotransformEq.(14)intoamoreusefulform, weblock-
||U||2
F diagonalize A and right-multiply by the matrix (cid:0)I I (cid:1),
I−I
Finding A for arbitrary C is a standard least-squares to show
problem; the optimal A satisfies CHCA = CHU. We
(U , U )=(C A , C A )+W(cid:48). (15)
expandEq.(10),twiceapplythisrelation,andrearrange + − + + − −
to get
Here we have introduced the symmetric and antisym-
tr(UHCC+U) metric forms U± = UL ±UR, C± = Cf ±Cbγ, and
w2 =1− , (11) A =A ±A .
||U||2 ± L,f R,f
F Eq. (15) takes the form of two independent least-
squares equations, each with half the dimension of
where C+ = (CHC)−1CH is the Moore-Penrose pseu-
Eq. (14). The two equations must be satisfied simul-
doinverse of C.
taneously, so to find the Bloch factors we can minimize
Using Eq. (11) and a numerical minimizer, the Bloch
factors that determine C may often be found to a useful ||U −C A ||2 +||U −C A ||2
levelofaccuracy(seeSec.IIIforimplementationdetails). w2 = + + + F − − − F, (16)
||U ||2 +||U ||2
In order to improve the accuracy and reliability of the + F − F
results, we impose further physical constraints. or equivalently
The PC impedance method4,12 assumes the unit cell
to be up-down symmetric, which causes the forward and tr(UHC C+U )+tr(UHC C+U )
w2 =1− + + + + − − − − . (17)
backward modes to be related. So far, we have only im- ||U ||2 +||U ||2
+ F − F
posed a relationship between the forward and backward
Bloch factors, not the modal fields within each unit cell. Again, this quantity may be minimized by a numerical
WecanhalvethenumberofunknownsinAandstrongly optimizer. The residual w2 for any solution to Eq. (17)
5
isequaltotheresidualobtainedbyinsertingthesolution III. NUMERICAL PROCEDURE
intoEq.(11): thetwoequationsdifferonlyinthesymme-
tryconstraintonbackwardmodalfields(Eq.(13)). Com-
Having outlined the theoretical basis of our method
paredtoEq.(11),wehaveremovedN M unknownsfrom
p forfindingtheBlochfactorsandimpedanceofaPCata
A (where N (cid:29) M is the number of sampled points in
p givenfrequency,incidentangle,andpolarization,wenow
eachunitcell),halvingitsdimensionatthecostofadding
provide some practical detail about our implementation
M unknowns to C as γ. These new unknowns must be
± ofthemethod. WeoutlinetheprocedureforM =3pairs
found simultaneously with the Bloch factors using a nu-
of Bloch modes.
merical minimizer, so it is important to supply a good
InCOMSOLMultiphysics4.2,wesimulatea1×8unit
starting estimate; our method for doing so is detailed in
cellsampleofPC,embeddedinitsbackgrounddielectric,
Sec. III.
with Bloch-Floquet periodic boundary conditions along
the two long boundaries (Fig. 1 shows a 1 × 5 struc-
ture). Eq. (15) is a set of LN equations, with 2M
C. Calculating impedance p
and MN unknowns in C and A respectively. To be
p ± ±
overspecified, the method requires LN > MN +2M;
p p
Once the Bloch factors and γ are known, the modal
thus L = 8 periods and a large N is sufficient to find
p
fields can be reconstructed and analyzed to determine
M = 3 modes. A deeper structure with more unit cells
thePC’simpedance. Theessentialquantitiesforthiscal-
doesnotnecessarilyprovideusefulinformationaboutad-
culation are the E and H field components in the plane
ditionalevanescentmodes,astheiramplitudedeepinside
of the PC interface (i.e., E and H , or E and H , de-
z x x z the structure may be negligible. From COMSOL we ex-
pendingonpolarization)ofeachBlochmodemalongthe
porttherelevantE andH fieldcomponentsintheL=8
left edge (y =0) of a unit cell (see Fig. 1). These quan-
unit cells, sampled over a 101×(50L+1) grid.
tities, E (x) and H (x), may be found from Eq. (15)
m m In order to compute a mode, it must be present in the
usingtheknownvaluesforC andC andinsertingthe
+ − structure with sufficient amplitude to be detected. Light
appropriate E or H fields into U and U .
+ − atnormalincidenceoftenfailstoexciteoddBlochmodes;
To calculate the impedance, we find the E and H ma-
these uncoupled modes26 consequently cannot be found
trices for the PC, as defined in Sec. IIA. Inserting mul-
byanoptimization,whichlosesaccuracyinsearchingfor
tiples of unit vectors c into Eq. (3), we can show that
+ modesthatarenotpresent. Atfrequenciesabovethefirst
Woodanomaly,thefrequenciesatwhichthehigherorder
Em(x)=Am(cid:88)Ep,m eikx(p)x, (18a) modesaremostimportant, thisproblemmaybeavoided
by exciting the PC slab not with a normally incident
p
plane wave, but with the first grating diffraction order.
Hm(x)=Am(cid:88)Hp,m eikx(p)x, (18b) TunhcisoutpeclehdnimquoedeisisusneodtirnelSeevca.nItVtoAaanpdarStieccu.laIVr pCr.obIfletmhe,
p it may instead be ignored.
where A is the amplitude of the normalized mode m, If we seek to find M = 3 Bloch modes, then find-
m
and E and H are the elements of E and H. It is ing a global minimum of Eq. (17) involves searching
p,m p,m
straightforward to exploit the orthogonality of the plane for 2M = 6 complex numbers. This is a hard prob-
wave grating diffraction orders to show that lem if attacked directly, but we use an algorithm that
gives more consistent success by providing a good start-
ing estimate. We start by minimizing the residual w2
AmEp,m =1/ax(cid:90) ax/2Em(x)e−ikx(p)xdx, (19a) in Eq. (11), which forces a relationship between forward
−ax/2 and backward Bloch factors but not the modal fields.
This involves finding only M complex numbers. As a
AmHp,m =1/ax(cid:90) ax/2Hm(x)e−ikx(p)xdx. (19b) statakretitnhgeersetsimulattoeffoarntehieghfobrowrainrgdBsimlouchlaftaiocnto,ros,rwtheeeiathnear-
−ax/2 lytically calculated Bloch factors for the dielectric back-
Eqs. (19) let us calculate each element of the E and H ground of the PC. At every step of the minimization,
matrices, up to a normalization constant A per col- evanescent modes are sorted into forward and backward
m
umn. We remove the constants by calculating the PC’s decaying modes, based on the moduli of their Bloch fac-
impedance(Eq. (4)or(5))withthePCitselfastherefer- tors. The minimization can be done by any standard
encematerial: byreciprocity-derivedBlochmodeorthog- numerical minimizer, such as SciPy’s27 fmin, which is a
onality relations12, this quantity should be the identity modifiedNelder-Meadoptimization28. Atthispoint, the
matrix. ThediagonalentriesofthismatrixaretheA 2; results are equivalent to those from the method of Ha et
m
the off-diagonal terms, which should be zero, provide an al.17,exceptthatwehavelessenedthelikelihoodofCbe-
error estimate. After normalizing the E and H matrices ing ill-conditioned by renormalizing the backward Bloch
forthePC,wecalculateitsimpedancematrixZ fromEq. factors µ in Eq. (7) and setting their phase origin to
m(cid:48)
(4) or (5) using a reference medium such as free space. the end of the PC.
6
Occasionally, we encounter an instability in which a
1.0 1.0
pair of modes have very large equal and opposite field
0.8 0.8
amplitudes and very small Bloch factors. When this oc-
curs, we follow a Gram-Schmidt-like process: we sub- |µ|00..46 |µ|00..46
tractthefieldofnon-problematicmodes(i.e.,modeswith 0.2 0.2
|µ|>10−3) from U and repeatedly minimize Eq. (11) to
(a)0.0 0.1 0.2 0.3 0.4 0.5 (b)0.0 0.1 0.2 0.3 0.4 0.5
find each of the remaining modes individually. a/λ a/λ
Using the solution to Eq. (11) as our estimate for the 0.5
Blochfactors,themodalfieldsmaybefoundwithaleast- 0.6
0.4 0.5
soqfubaarcekswoaprtdimainzdatfioornw.aTrdhemaovdereasggeivfieesldursaatnioeosftiemacahtepfaoirr λa/0.3 000...234Im(k)y
γ. Wenowhaveaplausibleestimateforγ andtheBloch 0.2 0.1
0.0
factors, which we can use as a starting estimate to mini- 0.1 0.4
miTzeoEfuqr.t(h1e7r).refine the estimates, we repeatedly iterate (c0).00.0 0.2 |a0.r4g(µ)0|/.6π 0.8 1.0 (d)0.|2Re0(.k4y)|0/.π6 0.8 0.10.20a./3λ
through the modes, fixing all but one µ and the corre-
FIG. 2: (Color online) Complex band structure for the
sponding element of γ, minimizing Eq. (17) to find the
PC. The Wood anomaly (a/λ=0.333) is marked. The
two variables. After this process, we finally minimize
modes are sorted into colors by |µ|; where two modes
Eq. (17) across all 6 complex dimensions simultaneously
are propagating (i.e., have |µ|=1), they are sorted by
toobtainthecorrectBlochfactorsandmodalfieldsfrom
|arg(µ)|. (a) Magnitude of Bloch factors |µ|, with three
which we calculate impedances. Forward and backward
Bloch modes found at all frequencies. (b) |µ| with two
propagating modes are sorted based on their flux15, be-
Bloch modes found below the Wood anomaly, three
fore impedances are calculated as outlined in Sec. IIC.
above. (c) Argument of Bloch factors. (d) Complex
band structure in 3D.
IV. APPLICATIONS
low frequencies (where the p=1 grating order is evanes-
We now apply our method to a range of typical prob-
cent and so may not be used to excite the structure,
lems. Each of these problems involves frequencies above
as mentioned in Sec. III) increases the accuracy of the
the first or second Wood anomaly—frequencies at which
other two modes (Fig. 2(b)). The complex arguments of
scalar methods fail and multiple modes are required to
the Bloch factors, which quantify how phase is acquired
describe the system. BlochCode, software that imple-
through propagation, are shown in Fig. 2(c), and the in-
ments our method in Python, using SciPy27 and Sage29,
formation about amplitude and phase is summarized in
is freely available on the internet20; we use it here.
a single plot in Fig. 2(d). Aside from slight errors in the
phase of strongly evanescent modes in Fig. 2(c), there is
good agreement between Fig. 2 and Bloch factors calcu-
A. Complex band structure
lated by highly accurate multipole techniques.
Figure 2 shows that at frequencies below the Wood
The first application of our method is to calculate the
anomaly there is at most one propagating Bloch mode,
complex band structure of a PC. The PC is a triangular
which becomes evanescent in the first bandgap with a
lattice of circular air holes with radius r = 0.3 a and
decay factor |µ| of no less than 0.5; it still decays far
lattice constant a = a in a dielectric background with
x more slowly than the other evanescent Bloch modes at
n=3. Wecalculatethebandstructureforlightpolarized
that frequency. Fig. 2(c) shows that for the evanescent
withtheHfieldoutofthePCplane(H polarization)at
z modes, either 0 or π phase is acquired across each unit
frequencies a/λ∈(0,0.5) in the Γ−M direction, i.e., at
cell.
normalincidence. UsingCOMSOL,wecalculatethefield
in an 8 period slab of the PC, and we apply our method
to find the largest three Bloch factors. w2 varies: it is
less than 10−8 at low frequencies and less than 10−4 at B. Antireflection coating
high frequencies.
Fig. 2 summarizes the propagation properties of the Our next application is to reproduce the design of an
two/three most dominant modes. The moduli of the antireflection coating we presented previously12, found
Bloch factors |µ|, which quantify how the modes’ ampli- using PC impedances calculated with a specialized
tudesvarywithpropagation,areshowninFigs.2(a)and transfer-matrixmethod24. Asinthispreviouspaper,our
2(b). Below the Wood anomaly, an inspection of A and design strategy is to try out a very large number of po-
γ showsthatthethirdmodeisbarelyexcitedbythenor- tential coatings, and choose the coating that gives the
mallyincidentplanewave,andthisreducestheaccuracy lowestreflectanceoffthecoatedstructure. TheuseofPC
oftheresults(Fig.2(a)). Ignoringtheuncoupledmodeat impedances makes this a feasible problem, as the evalu-
7
ation of each coating is quick, involving a few operations
on M ×M (here 3×3) matrices. 1.8
The target PC is a triangular lattice with lattice con-
stant a = a, consisting of air holes in a dielectric back-
x 1.6
ground with n = 2.86. The holes are cylinders with ra- 0.
7
dius r = 0.25 a. We seek to coat the PC to minimize 0.5
reflection for light with frequency a/λ = 0.38, incident 1.4 95
from air at an angle of 30◦ in the E polarization. At 3
z p
thisfrequencyandincidentangle,Mmin =2;weconsider / 1.2
atotalof3modestoensureaccuracy. Asinourprevious 2
work12, we seek a two-layer coating, where the degree ay2 0 0.7
of freedom is ay, the lattice vector component perpen- 1.0 .95 0.03.5 5
dicular to the air/PC interface. For a regular triangular
alantWdticeseimc,hauoyloa=stee√182231ap.ecrainoddsidoafteeaPcChsinwiCthOaMyS∈O[L0..6W,1e.8]ap√2p3lay 0.8 0.5 0.5 0.15 0.5
0.8 1.0 1.2 1.4 1.6 1.8
ourmethodtotheresultingdata,usingtheBlochfactors
of the previous PC as the starting estimate for the next.
a 2/ 3
BlochCode processes the 121 PCs in approximately 13 y1
minutes on a 3.06 GHz Intel Core 2 Duo desktop com- p
puter. An equivalent approach that only requires one FIG. 3: (Color online) Reflectance of the coated PC as
PC to be evaluated is detailed in Sec. IVC; we do not a function of ay1 and ay2, the relative thicknesses of the
use it here since the purpose of this section is to demon- two coating layers, calculated using PC impedances
strate the reliability and consistency of the optimization from BlochCode. The minimum reflectance is marked.
procedure.
Wethencalculatethereflectancesoffthe1212 =14641
coated stacks (Fig. 3), which takes 34 seconds on a sin-
con (n = 3.518) PC with r = 0.45 a is self-collimating
gle core of the desktop computer. The optimal coat-
√ forbothpolarizationsatnormalincidence. Thelargera-
ing is found√ to have thicknesses ay1 = 1.53 23a and dius is an extreme case that is challenging to simulate
ay2 =0.65 23a, and reduces the reflectance of the struc- accurately. At this frequency Mmin =3, so for Ez polar-
ture from R =0.945 to R =1.96×10−4. The results in ized light we include M = 3 modes in our calculations,
Fig. 3 agree well with data calculated by a highly accu- with light incident from the p = 1 grating order so that
rate multipole scattering matrix method: the RMS dif- the otherwise uncoupled mode is excited. For Hz light,
ference is 3.4×10−3, and the only noticeable differences thisproceduredoesnotyieldaccurateresults—Blochfac-
occur on the two sharp resonant features near the lower tors are calculated accurately, but the calculated reflec-
edgeofthefigure. Specifically,themultipole-basedcalcu- tion coefficients differ from those calculated directly in
lationsshowthatthecoatingreducesthePC’sreflectance COMSOL. The calculated impedances prove sufficiently
from R=0.943 to R=4.29×10−4. accuratetodesignaneffectiveantireflectioncoating, but
the inaccuracies mean that the coating is not optimal.
To avoid these inaccuracies in H polarization, we ex-
z
C. All-polarization antireflection coating ploitthesymmetrythatcausestheuncoupledmode. The
physical structure and normally incident field are both
symmetric about the y-axis, and so modes without even
Finally, we apply our methods to find an all-
symmetry are not coupled to. Therefore we formally ig-
polarizationantireflectioncoatingforasilicon-basedself-
nore the uncoupled odd mode, in each PC and in the
collimating square-lattice photonic crystal presented by
Park et al.21. They investigated this class of structures reference medium, setting M = 2. In our Hz COMSOL
simulations for this structure, light is normally incident.
using a scalar treatment of reflections, and were able to
design an all-polarization coating at a/λ = 0.28, be- In Fig. 2 of Park et al.’s paper21, they state that
low the first Wood anomaly. Since their scalar treat- R (cid:39) 0.28 for Ez polarized light, and R (cid:39) 0.35 for Hz
ment does not support multiple propagating or evanes- light. We calculate with BlochCode that a semi-infinite
cent Bloch modes, it generally does not work above the slab of the PC has R = 0.284 for Ez, and R = 0.354 for
Wood anomaly. Our method does not have this limita- Hz polarized light at this frequency, when incident from
tion and we demonstrate this by designing an antireflec- silicon. Specialized FEM-based transfer-matrix calcula-
tion coating for both polarizations at a frequency well tions agree, showing R = 0.284 for Ez polarization, and
above the Wood anomaly, using more than one Bloch R=0.357 for Hz polarization.
mode. At a/λ=0.368, normally incident light is reflected by
Park et al.21 showed that at a/λ = 0.368, a 2D sili- the PC into three propagating diffraction orders. Due to
8
the symmetries of the problem, the ±1 orders are only thelargenumberofcandidatecoatings(∼1.3×107),the
excitedinanevensuperposition, solightisreflectedinto embarrassingly parallel problem was split over 16 cores
twomodes. Asuccessfulcoatingneedstosuppressreflec- of the workstation, taking approximately 80 minutes per
tion into both these modes simultaneously, and so must polarization.
balance two modes’ amplitudes and two modes’ phases The best E coating reduces R from 0.284 to 9.56×
z
simultaneously for each polarization. Thus the design of 10−5, and the best H coating reduces R from 0.354 to
z
a perfect all-polarization coating requires 8 continuous 3.33×10−4. The best all-round coating is taken to be
degrees of freedom. Rather than trying to search an 8- the one with the lowest total reflection in the two po-
dimensional parameter space, which is computationally larizations. This coating has r = 0.13 a, d = 0.89 a,
1 1
expensive even when the evaluation of each point is effi- r =0.17 a, and d =0.90 a (Fig. 4). In E it reduces R
2 2 z
cient, we consider coatings with four degrees of freedom to0.0141,andinH itreducesRto0.0197. Calculations
z
andacceptthatweareunlikelytofindanall-polarization from a specialized transfer matrix method24 agree with
coatingwithzeroreflectance. Nevertheless, thisisapar- these results, giving R = 0.0142 in E polarization and
z
ticularly difficult problem: not only do we need many R=0.0211 in H .
z
degrees of freedom to find a satisfactory coating, but if To verify these results without the aid of our special-
either of the Bloch factors in a PC is incorrect or any el- izedmethods, implementationsofwhicharenotpublicly
ement of the PC’s impedance matrix is wrong, then the available,wesimulatethestructureusingCOMSOLMul-
calculated net reflection off the structure is incorrect as tiphysics. Since COMSOL cannot directly calculate re-
well. flection coefficients off semi-infinite PCs, we simulate a
Tolimitthecoating’sthickness,weembedthefourde- 20-period section of the uncoated PC surrounded by the
grees of freedom into two rows of holes by varying both backgrounddielectric,andcomparetheresultstoasimu-
theholeradii,r1andr2,andthespaceafterthelayers,d1 lationwiththeantireflectioncoatingonbothsidesofthe
andd2(Fig.4). Increasingd1andd2issimilartoincreas- PC section. BlochCode calculates the reflectance of the
ing ay, as in Sec. IVB, but because the candidate PCs uncoated and coated structures to be 0.407 and 0.0124
areindependentofd,onlyonePCperradiusneedstobe respectivelyintheE polarization,and0.574and0.0074
z
simulated in COMSOL. Furthermore, the properties of in the H polarization. The COMSOL simulations agree
z
the layers of silicon with thickness di may be calculated with these results, showing that the coating reduces R
analytically. We consider 36 possible hole radii in the from 0.407 to 0.0129 in the E polarization, and from
z
range ri ∈ [0.10,0.45] a and 99 values of di ∈ (0,1) a. 0.585 to 0.0055 in the Hz polarization.
To allow a thin coating, we set a = 2r+0.1a for each
y
PC. If necessary, additional degrees of freedom could be
added to find a coating with even lower reflectances.
V. DISCUSSION & CONCLUSION
y
We have detailed a method for calculating the com-
plexbandstructureandimpedanceofPCs. Themethod
takes into account structural symmetries in the PC, and
x enforces relationships between the fields of forward and
backward modes, thus improving the method’s accuracy
r r by eliminating ill-conditioning and constraining modal
1 2
fields. We have applied the method to three cases, and
have demonstrated that it works for a variety of square
and triangular lattice 2D photonic crystals, for light in
both polarizations and at different incident angles. We
d d
1 2 havedemonstratedthatourmethodworksatfrequencies
both above and below the first Wood anomaly, the fre-
FIG. 4: Schematic of the all-polarization antireflection
quency above which scalar methods cannot adequately
coating. r and r are the radii of the holes in the first
1 2 describe light propagation and reflection in PCs.
two layers, and d and d are the thicknesses of the
1 2 ThestrongertheexcitationofaBlochmode,themore
extra silicon background layers between the first few
accuratelyourmethodcalculatesitsproperties. Thusthe
rows of holes. For this coating, r =0.13 a, d =0.89 a,
1 1 method is well-suited to calculating reflection and trans-
r =0.17 a, and d =0.9 a.
2 2 missionthrougharbitraryPCstacks,wherethemostim-
portant modes are those that are strongly excited. Since
On a single core of a 16×2.4 GHz Intel Xeon-Quad PC impedances make it so easy to calculate the reflec-
workstation, it took a total of 15 minutes to find the tion and transmission properties of many combinations
modes of the 36 PCs in the two polarizations. For E of PCs in a stack, it is feasible to search large parameter
z
polarization,w2 (cid:39)10−5 formostradii,andforH polar- spaces of PC stacks for particular reflective properties
z
izationw2rangedroughlyfrom3×10−3forthinunitcells over a range of frequencies, incident angles and polariza-
to 10−7 for the thicker cells with larger radius. Due to tions. The method can be used to design not only all-
9
polarization antireflection coatings, but also broadband 3N.Stefanou,V.Karathanos,andA.Modinos,J.Phys.Condens.
antireflection coatings12, polarization filters, angular fil- Matter4,7389(1992).
ters, and other devices. 4F.J.Lawrence,L.C.Botten,K.B.Dossou,andC.M.deSterke,
Appl.Phys.Lett.93,1114(2008).
Ha et al. have applied their method to slab PC 5C. M. de Sterke, K. B. Dossou, T. P. White, L. C. Botten, and
waveguides19. We have not yet applied our method to R.C.McPhedran,Opt.Express17,17338(2009).
any 3D structure. As long as the x − z plane mirror 6V.Heine,Surf.Sci.2,1(1964).
symmetry is present, our method for finding the com- 7R. Biswas, Z.-Y. Li, and K. M. Ho, Appl. Phys. Lett. 84, 1254
(2004).
plex band structure remains valid. The field of a slab
8W. Smigaj, P. Lalanne, J. Yang, T. Paul, C. Rockstuhl, and
waveguide might be sampled only over the PC’s surface F.Lederer,Appl.Phys.Lett.98,111107(2011).
(as in a SNOM experiment19) or throughout the entire 9C.R.Simovski,Metamaterials1,62(2007).
volume of the structure (as in a simulation); either case 10M.Miri,A.Khavasi,K.Mehrany,andB.Rashidian,Opt.Lett.
35,115(2010).
provides sufficient information to determine the modal
11T.-T. Kim, S.-G. Lee, M.-W. Kim, H. Y. Park, and J.-E. Kim,
fields within the sampled region and the associated com-
Appl.Phys.Lett.95,011119(2009).
plex band structure. However, the impedance formalism 12F. J. Lawrence, L. C. Botten, K. B. Dossou, C. M. de Sterke,
is yet to be developed for 3D structures. andR.C.McPhedran,Phys.Rev.A80,23826(2009).
13F. J. Lawrence, L. C. Botten, K. B. Dossou, R. C. McPhedran,
OurmethodisalsovalidforfindingmodesofPCwave-
andC.M.deSterke,Phys.Rev.A82,053840(2010).
guides, using supercells. Calculation of reflection and
14B.Gralak,S.Enoch,andG.Tayeb,J.Opt.Soc.Am.A17,1012
transmission matrices between PC waveguides is yet to
(2000).
bedemonstratedusingimpedances, buttheyhaveprevi- 15L. C. Botten, N.-A. P. Nicorovici, R. C. McPhedran, C. M.
ouslybeencalculateddirectlyfromthesupercell’sEand deSterke,andA.A.Asatryan,Phys.Rev.E64,046603(2001).
H matrices5. 16Y.-C.HsueandT.-J.Yang,Phys.Rev.E70,016706(2004).
17S. Ha, A. A. Sukhorukov, K. B. Dossou, L. C. Botten, C. M.
Bloch mode analysis is a valuable tool in understand-
deSterke,andY.S.Kivshar,Opt.Lett.34,3776(2009).
ing light’s interactions with PCs. Using an EM solver 18A. A. Sukhorukov, S. Ha, I. V. Shadrivov, D. A. Powell, and
and our method, for which source code is available20, it Y.S.Kivshar,Opt.Express17,3716(2009).
isstraightforwardtofindaPC’scomplexbandstructure 19S. Ha, M. Spasenovi´c, A. A. Sukhorukov, T. P. White, C. M.
de Sterke, L. K. Kuipers, T. F. Krauss, and Y. S. Kivshar, J.
and its impedance. Respectively, these quantities dic-
Opt.Soc.Am.B28,955(2011).
tate how the Bloch modes travel through the PC, and 20https://launchpad.net/blochcode
whichmodestheycouplewithataPCinterface. Ifthese 21J.Park,S.-G.Lee,H.Park,andM.Lee,Opt.Express18,13083
quantities are known for a set of PCs, then it is fast and (2010).
22L.C.Botten,N.-A.P.Nicorovici,A.A.Asatryan,R.C.McPhe-
efficient to calculate how light travels through arbitrary
dran,C.M.deSterke,andP.A.Robinson,J.Opt.Soc.Am.A
stacks of the PCs.
17,2165(2000).
23R. C. McPhedran, N.-A. P. Nicorovici, L. C. Botten, and K. A.
Grubits,J.Math.Phys.41,7808(2000).
24L.C.Botten,T.P.White,A.A.Asatryan,T.N.Langtry,C.M.
ACKNOWLEDGMENTS deSterke,andR.C.McPhedran,Phys.Rev.E70,56606(2004).
25E. Popov, L. Mashev, and D. Maystre, J. Mod. Opt. 33, 607
(1986).
This research was conducted by the Australian Re- 26K.Sakoda,Phys.Rev.B52,7982(1995).
search Council Centre of Excellence for Ultrahigh band- 27E. Jones, T. Oliphant, P. Peterson, et al., SciPy: Open source
scientific tools for Python (2001–),http://www.scipy.org/.
width Devices for Optical Systems (project number
28M. H. Wright, in Numerical Analysis 1995, Proceedings of the
CE110001018).
16th Dundee Conference, edited by D. F. Griffiths and G. A.
Watson (Addison Wesley Longman, Harlow, Essex, UK, 1996),
vol.344ofPitmanResearchNotesinMathematics,pp.191–208.
1R.W.Wood,Philos.Mag.4,396(1902). 29W. Stein et al., Sage Mathematics Software (Version 4.7), The
2R.Smaˆali,D.Felbacq,andG.Granet,PhysicaE18,443 (2003). SageDevelopmentTeam(2011),http://www.sagemath.org.
