PDE Constrained Optimization and Design of Frozen ... - Infoscience
COMMUNICATIONS IN COMPUTATIONAL PHYSICS Vol. 1, No. 1, pp. 1-21
Commun. Comput. Phys. December 2006
PDE Constrained Optimization and Design of Frozen Mode Crystals S. Chun1 , J. S. Hesthaven1 1
Division of Applied Mathematics, Brown University, Providence, RI 02912
Abstract. We explore the use of PDE constrained nonlinear optimization techniques to optimize and design electromagnetic crystals which exhibit frozen mode behavior. This is characterized by the association with Van Hove singularities in the dispersion relation, e.g., stationary reflection points and degenerate band edge points. Hence, the optimization process modifies the dispersion relation by adjusting the geometries and material parameters. The resulting algorithm is found to be capable of recovering all known crystal configurations as well as many new configurations, some of which display dramatically improved properties over previously used configuration. We investigate both gyrotropic photonic crystals and degenerate band edge crystals as well as the more complex case of the oblique incidence where we extend the investigation to the three-dimensional case to identify the first three-dimensional crystal exhibiting frozen mode behavior. Key words: PDE constrained optimization, Method of Matched Asymptotics, Frozen mode, Photonic crystal
1
Introduction
Controlling the flow of electromagnetic waves in materials has been the frontier of our technologies in the last decade. Along with prohibiting wave propagation or localizing wave in a certain area, slowing down wave at the significant rate has grabbed broad attentions. Recent studies of periodic structures of several different anisotripic materials shows that when a monochromatic wave with a certain frequency propagates in such structures, it shows abnormal electromagnetic phenomena such as dramatic slow-down of the wave, the significant increase of field amplitude and singularity of the trasmittance rate, which are classified as the frozen mode. These unique features of the frozen mode are attractive for numerous practical applications and recent studies [2] have confirmed
∗
Correspondence to: Division of Applied Mathematics, Brown University, Box F, Providence, RI 02912. Email: [email protected]; [email protected] http://www.global-sci.com/
1
c !2006 Global-Science Press
2
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
the robustness of these phonomena. In the beginning, the frozen mode was regarded as a distinctive phenomenon only related to stationary inflection points of the dispersion relations ω(k) such as, ∂ω = 0; ∂k
∂2ω = 0; ∂k2
∂3ω != 0 ∂k3
at ω = ω0
However, it is worth mentioning that these properties of the dispersion relation are special cases of what is known as Van Hove singularities [7] as points where the first derivative of the dispersion relation ω(k) vanishes, i.e., ∂ω = 0, ∂k
at ω = ω0
This nomenclature is related to density of state(DOS), since it is well-known that DOS is inversely proportional to ∂ω/∂k [8]. However, all previous studies have focused on the analysis of crystal configurations found by a trial-and-error technique which tends to limit the number of parameters one can freely vary. This limits the flexibility in the design of such crystals and leaves open the question of whether it is generally possible to achieve the sought-after properties in the dispersion relation, e.g., when given some specific materials, one can modify the geometric properties of the crystal as needed. To address this problem in a more systematic way, we consider the use of a PDE constrained optimization approach where we use the dispersion relation, given by the transfer matrix method or by solving the Maxwell eigenproblem, in combination with the algorithm of Method of Moving Asymptote(MMA) [9–13] to recreate and optimize known configurations as well as to design entire crystal arrangements. In this study, we consider the optimization and the design of three different types of crystals, all exhibiting the frozen mode phenomena. The first one is a gyrotropic photonic crystal consisting of two misaligned anisotropic layers and one magnetic layer. Its dispersion relation has a stationary inflection point ω0 on the 2nd band. The second class of crystals consists of two misaligned anisotropic layers and vacuum between them. Its dispersion relation has a degenerate band edge ωd in the first band. The final structure consists of only one anisotropic layer with εxz != 0 and vacuum between them. However, the frozen mode phenomenon only emerges when the incident wave propagates into this layer at an oblique angle. This latter case is quite different from the previous two cases. The dispersion relation where a stationary inflection point ω0 lies was originally computed by a one-dimensional transfer matrix method using tangential consistency of wave vector. To verify the existence of a stationary inflection point in general three dimensional setting which was previously unascertained, we first solve the general Maxwell eigenproblem to get the dispersion relation with initial parameter sets and apply the optimization techniques to seek a stationary inflection point in the dispersion relation of three dimension. The remaining part of the paper is organized as follows. In Sec. 2, we briefly recall the physical model, Maxwell’s equations and how one recovers the dispersion relation from
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
3
these equations. This sets the stage for Sec. 2 where we explain the implementation of the optimization method based on Method of Moving Asymptote(MMA). In Secs 3 and 4, we display the methods and results of optimizing the magnetic frozen mode and degenerate band edge for known designs as well as new configurations, respectively. Section 5 is devoted to the design of crystals for the oblique case and its modeling in the time domain. We conclude with a few remarks and future works in Sec. 6.
2
Physical Model
We consider Maxwell’s equations ˜ =− ∇×E
˜ dB , dt
˜ = ∇×H
˜ dD dt
˜ and H ˜ is the electric field and magnetic field, respectively, D ˜ and B ˜ the displacewhere E ment field and the induction field. These fields are related by linear constitutive relations as follows. ˜ ˜ D(z) = εˆ(z)E(z),
˜ ˜ B(z) =µ ˆ(z)H(z)
where we have assumed that (ˆ ε, µ ˆ) are homogeneous in the (x, y)-plane. If we consider harmonic solutions with frequency ω, we recover the Maxwell’s equations in the frequency domain. ˜ = iω µ ˜ . ∇×H ˜ = −iω εˆ(z)E ˜ ∇×E ˆ(z)H
(2.1)
To simplify the above equations, we introduce the normalized quantities using z=
z˜ , ˜ L
t=
t˜ ˜ c˜0 L/
˜ is a reference length and c˜0 = (˜ where, L ε0 µ ˜0 )− 2 is the dimensional speed of light in ˜ and H ˜ are normalized as vacuum. Also, E ! ˜ ˜ ε˜0 E H , H= E= ˜0 ˜0 µ˜0 H H 1
˜ 0 is a dimensional reference magnetic field strength. By substituting these norwhere H malized quantities into Eq.(2.1), we obtain the following non-dimensionalized equations ∇ × E = iω µ ˆ(z)H . ∇ × H = −iω εˆ(z)E , where (ˆ ε, µ ˆ) now represent the relative permittivity and permeability, respectively.
4
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
If we now assume that the solution is periodic and seek solutions of the form # " # " $ x) $k E($ E $ $k = wave number $ k exp(ik · $x) , $ x) = H H($ we recover the modal equation ˆ(z)Hk , (∇ + i$k) × Ek = iω µ
(∇ + i$k) × Hk = −iω εˆ(z)Ek ,
which simplifies as (∇ + i$k) ×
1 (∇ + i$k) × Ek = ω 2 εˆ(z)Ek µ ˆ(z)
For a given value of the wave vector $k, this is recognized as an eigenvalue problem for (ω 2 , Ek ), from which we recover the dispersion relation ($k, ω($k)). Hence, when optimizing the crystals, we must ensure that all considered solutions satisfy the above eigenvalue problem which emerges as the constraints in the optimization process. The goal of the optimization is to identify one type of Van Hove singularities [7] in the dispersion relation, leading to a very strong divergence of the density of states(DOS) near the corresponding frequency. There are several types of singularities, but we are specially interested in finding singularities where the first and second derivative are both zero, i.e. ∂ω = 0, ∂k
∂2ω = 0, ∂k2
at ω = ω0
We state this equivalently as the following minimization problem. $ $ $ 2 $ $ dω $ $d ω $ $ $ $ minimize λ $ ($x)$ + $ 2 ($x)$$ , at $k = $k0 dk dk ≤ xj ≤ xmax , xmin j j
(2.2)
1≤j≤n
where, $x is a n-variable vector and λ is a weight constant. The objective function dω x) dk ($ 2 and ddkω2 ($x) can be obtained directly from the dispersion relation with optimization variable $x by solving the eigenvalue problem discussed in Sec. 2. Once the dispersion relation is obtained, the derivatives of the objective functions can be computed easily. For example, the fourth order approximation of its first and second derivatives are computed by difference approximations such as, & & % % 1 dωj+1 dωj−1 dωj−2 ∂ dω dωj+2 ($x) = +8 −8 + − ∂xj dk 12∆xj dk dk dk dk & & % % 2 dωj+1 dωj dωj−1 dωj−2 dωj+2 1 dω ∂ ($x) = + 16 − 30 + 16 − − dk dk dk dk dk ∂x2j dk 12∆x2j , where
dω dωj = ($x + j∆$x), dk dk
1≤j≤n
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
5
In fact, this fourth order approximation is preferred in the computation of the derivatives to increase the accuracy and to make the fast convergence. Furthermore, it is required to obtain an initial guess $k0 or at least a group of possible candidates of $k0 s. Identifying suitable $k0 is one of most critical and significant factors of convergence of the optimization solution. We have identified the following approach to work well. 1. Divide the dispersion relation into several sections such that each section contains $ only one wave vector with | dω dk (k)| < ε & 1. 2. For each section, find one wave vector $k0 such that D12 ($k0 ) = min | i
i+2 ' dω $ (k" )|, dk
"=i−2
i ∈ section I
Note that the function inside the bars is a way of indicating the amount of the first and the second derivative at the same time 3. Optimize the objective function at $k0 . Since the objective function is generally not convex, we use the algorithm of Method of Moving Asymptote(MMA). The application of MMA to find singularity points of dispersion relation can be stated briefly as follows. First, Eq.(2.2) can be changed into the following equivalent problem. minimize
f0 (x) = z +
4 '
(2.3)
ci yi
i=1
subject to
dω d2 ω ($x) − y1 ≤ 0, f2 (x) = ($x) − y2 ≤ 0 dk dk2 dω d2 ω f3 (x) = −λ ($x) − y3 ≤ 0, f4 (x) = − 2 ($x) − y4 ≤ 0 dk dk z ≥ 0, yi ≥ 0, 1 ≤ i ≤ 4, f1 (x) = λ
≤ xj ≤ xmax , xmin j j
1≤j≤n
where, y1 , · · · , yn and z are artificial optimization variables. Note that at the optimized solution $x, z is zero and yi is given as the maximum of each constraint function. Now, we have the non-linear optimization problem with several constraints. For these fi ($x)s, let’s define gi (x) in the following way, gi (x) = ri +
n ' j=1
(
pij qij + ) Uj − xj xj − L j
pij = max{0, (Uj − xj )2 ri = fi (x) −
n ' ( j=1
∂fi }, ∂xj
(2.4)
qij = max{0, −(xj − Lj )2
pij qij + ), Uj − xj xj − L j
1 ≤ i ≤ 4,
∂fi } ∂xj
1≤j≤n
6
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
The parameters Lj and Uj (known as moving asymptotes) are chosen as Lj < xj < Uj and play a critical role in the convergence of the optimization scheme. In fact, Lj and Uj are adjusted for each iteration according to the convergence of the solution. Note that the second derivatives of gi are given as % & 2pij 2qij ∂ 2 gi = δj" + , 1 ≤ j, ' ≤ n ∂xj x" (Uj − xj )3 (xj − Lj )3 Since pij ≥ 0, qij ≥ 0, gi is a convex function. In fact, as long as Lj and Uj are finite, ∂gi gi is strictly convex in all variables except for ∂x (x) = 0, 1 ≤ j ≤ n. With gi s as in the j Eq.(2.4), we have the following subproblem. n % '
& qij pij + + r0 Uj − xj xj − L j j=1 % & pij qij + + ri ≤ fi , Uj − xj xj − L j
P : minimize subject to
≤ xj ≤ xmax , xmin j j
(2.5) 1≤i≤4
1≤j≤n
Then, the Lagrangian function corresponding to this problem is given by y T $r + W (x, y) = $r0 − $
& n % ' p0j + $y T p$j q0j + $y T $qj + Uj − xj ($y ) xj ($y ) − Lj j=1
In fact, W (x, y) is differentiable, smooth and concave. Finally, the Lagrangian dual problem with this Lagrangian function is as follows. D : maximize
& n % ' q0j + $y T $qj p0j + $y T $pj + W (x, y) = $r0 − $y $r + Uj − xj ($y ) xj ($y ) − Lj T
j=1
subject to
y≥0
This yields a nice convex optimization problem and it can be solved by a primal-dual Newton method. The details for this primal-dual method are given in [13]
3
Optimization and Design of Magnetic Crystals
The frozen mode in the gyrotropic crystal occurs when a monochromatic wave with frequency ω0 propagates into a periodic array consisting of two misaligned anisotropic layers and one magnetic layer(Fig. 1). In the dispersion relation, this ω0 is a stationary inflection point with the following properties ∂ω = 0, ∂k
∂2ω = 0, ∂k2
∂3ω != 0, ∂k3
at ω = ω0
(3.1)
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
7
As shown in [4], one set of parameters parameters resulting in an inflection point at ω0 is
For A layer : εˆ =
(
)
13.61 + 12.40 cos(2φ) 12.40 sin(2φ) 12.40 sin(2φ) 13.61 − 12.40 cos(2φ)
µ ˆ=I
,
(
) ( ) 5.0 0.0 60.0 i37.0 For F layer : εˆ = , µ ˆ= 0.0 5.0 −i37.0 60.0 φ1 = 0, φ2 = −π/4, LF = 0.0047454, LA = 0.5(1.0 − LF ) k0 = 2.6329,
ω0 = 0.607676756(c/L)
where φ1 , φ2 is the misalignment angle of the first and second anisotropic layer. Also LF is the thickness of F layer and LA1 , LA2 are those of the first and second anisotropic layer, respectively. The first test is to recover this set of parameters by the aforementioned optimization technique and, subsequently, seek sets of parameters if possible. Since the periodicity is along z direction only, the dispersion relation of the magnetic frozen mode can be obtained by the transfer matrix method rather than solving the full Maxwell eigenvalue problem. The transfer matrix method is a way of finding the dispersion relation from Bloch solutions of Maxwell’s equations in planar structures. The transfer matrix(TL (ω)) can be obtained with bases of tangential components of electrical and magnetic fields. Then, it is used for Bloch solutions of the Maxwell’s equations in a periodic media such as qk (z + L) = eikL qk (z),
qk = [Ex Eu Hx Hy ]Tk
Using the property of transfer matrix TL (ω)qk (z) = qk (z + L), we get the following charicteristic equation. (TL (ω) − eikL )qk (L) = 0 =⇒
det(TL (ω) − eikL ) = 0
0.8 0.7 w0 0.6 0.5 0.4 0.3 0.2 0.1 0
0
1
2
3
4
5
6
Figure 1: Periodic structure of gyrotropic photonic crystals consisting of two different anisotropic materials with one magnetic layer(left) and its dispersion relation(right).
8
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Dispersion relation when rho0= 0.018924
Dispersion relation when rho0= 0.0095342 0.7
0.6
0.6
0.6
0.5
0.5
0.5
0.4
0.4
0.4
0.3
freqency
0.7
freqency
freqency
Dispersion relation when rho0= 0.05 0.7
0.3
0.3
0.2
0.2
0.2
0.1
0.1
0.1
0
0
1
2
3 4 wave vector
5
6
0
7
0
1
2
3 4 wave vector
5
6
7
0
0
1
2
3 4 wave vector
5
6
7
Figure 2: Scenes frm optimization process with one variable LF (left to right). LF =0.0244(left), 0.009453(cen2 ter), 0.00474(right) after 6 iterations, dω = 3.97d-13, ddxω2 = 0.007907 with λ = 10000. dx
which yields the dispersion relation between ω and the wave vector k. The transfer matrix of an anisotropic layer with misalignment angle φ (TA ) and a magnetic layer (TF ) are given as TA (φ, LA ) = W (φ, LA )W −1 (φ, 0), TF (LF ) = W (LF )W −1 (0) (cos φ)ein1 a (cos φ)e−in1 a −(sin φ)ein2 a −(sin φ)e−in2 a in a −in a in a 1 1 2 (sin φ)e (sin φ)e (cos φ)e (cos φ)e−in2 a W (φ, LA ) = −η1 (sin φ)ein1 a η1 (sin φ)e−in1 a −η2 (cos φ)ein2 a η2 (cos φ)e−in2 a η1 (cos φ)ein1 a −η1 (cos φ)e−in1 a −η2 (sin φ)ein2 a η2 (sin φ)e−in2 a ein1 f e−in1 f −iein2 f −ie−in2 f −iein1 f −ie−in1 f ein2 f e−in2 f W (LF ) = iη1 ein1 f −iηe−in1 f −η2 ein2 f η2 e−in2 f η1 ein1 f −ηe−in1 f −iη2 ein2 f iη2 e−in2 f
where, a =
ω c LA
and f =
ω c LF .
Then, the transfer matrix of three-cell layer is
TL (φ1, φ2, LA1 , LA2 , LF ) = TF (LF )TA (φ2 , LA2 )TA (φ1 , LA1 ) and the dispersion relation of this layer is obtained by solving det(TL − ξI) = 0,
ξ = eikL , L = LA1 + LA2 + LF .
The left hand side is a fourth order polynomial and we can compute the roots of it by simple root-finding. As for choosing wave vector k0 corresponding to ω0 , we first find the wave vector kb corresponding to a band edge ωb (the highest peak such that dω dk |ω=ωb = 0) and seek k less than kb where the minimum of D12 occurs as discussed above. However, this may still miss the true candidates of wave vector k0 for optimized solution ω0 . To address this issue, we
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Input
Solution
n variables φ2 LF LA1 εA δA ω0 k0 dω dk |ω=ω0 d2 ω | dk 2 ω=ω0 TR rate
Set1 1 LF , (LA1 = LA2 ) −π/4 0.00474 0.49763 5.1 1.1 0.60769 2.63894 3.9666d-13 0.007907 0.3388
Set2 2 LF , LA1 −π/4 0.00480 0.44555 5.1 1.1 0.60950 2.63894 2.41590d-16 0.00898 0.4081
Set3 2 LF , LA1 −π/3 0.00490 0.22056 5.1 1.1 0.67177 2.70177 6.5758d-16 0.00917 0.7786
9
Set4 4 LF , LA1 , εA , δA −π/4 0.00470 0.45001 5.0 1.18189 0.62171 2.63894 1.80310d-12 0.00191 0.3109
Table 1: Optimized solutions with different number of variables or different initial values for magnetic frozen mode. φ2 =misalignment angle of second anisotropic layer(we let φ1 = 0). LF = thickness of magnetic layer, LA1 =thickness of the first anisotropic layer. εA and δA are used in the permittivity of anisotropic layer as follows. ( )
εˆ =
εA + δA cos(2φ) δA sin(2φ) , TR rate = Transmittance rate δA sin(2φ) εA − δA cos(2φ)
change the initial values if the converged solution in the sense of *x − xold * < ε0 does not satisfy dω d2 ω < ε1 , < ε2 . dk dk2 In general, we shall allow the thickness of the magnetic layer(LF ), the thickness of the first anisotropic layer (LA1 ) and material constants of the anisotropic layer(εA , δA ) to vary during the optimization, although this could be extended as needed. In Fig. 2, we show the sequence of convergence of the dispersion relation for a stationary inflection point using optimization with just one variable, LF . We note that the optimized solution for LF is accurate up to 5.0e-06 with the analytically derived LF [4] and the difference is negligible in the sense of sensitivity. Other optimized solutions with different number of variables(n) are shown in Table 1. We observe that the thickness of F layer is approximately the same for all sets, but the thickness of A layer varies depending on the misalignment angle or the material constants. The absolute values of the first and second derivative at each ω0 are sufficiently close to zero, so all sets show every characteristics of the frozen mode. However, the closer to zero the first derivative is, the slower the group velocity(velocity of energy transfer) of the wave is. Also, the intensity of the field amplitude inside the layer and the amount of energy density stacked inside the layer indicated by energy density flux difference(EDFD) [2] are closely related to the transmittance rate. If a larger amount of energy enters into layer, the growth of the field amplitude gets higher and energy density
10
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
EDFD 0.45
90
0.4
80
0.35
70
0.3
60
0.25
EDFD
Max. of Field amplitude
Max. of Field amplitude 100
50 40
0.2 0.15
30
0.1
20
0.05
10 0
Set1 Set2 Set3 Set4
0 0
100
200
300
400
500
!0.05
0
100
200
time
300
400
500
time
Figure 3: Maximum of the field amplitude(left) and the Energy Density flux difference(ef of the first cell - ef of the last cell, ef =energy density flux) (right) for set 1(solid line), set 2(dashed line), set 3(dot-dashed line) and set 4(dotted line) of Table 1. Number of cell = 16. E"x polarization.
inside the layer gets stronger (Fig. 3).
4
Optimization and Design of Degenerate Band Edge Crystals
Different from the stationary inflection points in the gyrotropic crystals, the third derivative of the dispersion relation also needs to be zero in the degenerate band edge crystals, i.e. ∂ω = 0, ∂k
∂2ω = 0, ∂k2
∂3ω =0 ∂k3
at ω = ω0
According to the paper [17], the following set of parameters generate a degenerate band
1.2 w0 1
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
Figure 4: Periodic crystal structure allowing a degenerate band edge. The crystal consists of two misaligned anisotropic layers with vacuum between them(left) and its dispersion relation(right)
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
11
edge point at k0 = π. (
13.61 + 12.40 cos(2φ) 12.40 sin(2φ) For A layer : εˆ = 12.40 sin(2φ) 13.61 − 12.40 cos(2φ) φ1 = 0, φ2 = π/4, LV = 0.45891, LA = 0.5(1.0 − LV ) k0 = π,
)
,
µ ˆ=I
ω0 = 1.0094052895(c/L)
,where LA and LV are the thiness of anisotropic layer and vacuum, respectively. One of main advantages of this configuration is the lack of a magnetic components which are both complicated and difficult to manufacture with a sufficiently low loss. As shown in Fig. 4, the dispersion relation is indeed flat at the degenerate band edge point and it increases the speed of convergence when minimizing $ 2 $ $ $ $d ω $ $ dω $ J = λ $$ $$ + $$ 2 $$ dk dk As for the magnetic case, the one-dimensional nature of the problem enables the use of the transfer matrix method, which is expressed as for the three cells. TL (φ1, φ2, LA1 , LA2 , LV ) = TA (φ2 , LA2 )TV (LV )TA (φ1 , LA1 ) where
eiv e−iv −eivf iv e e−iv eiv TV (LV ) = W (LV )W −1 (0), W (LV ) = −eiv e−iv −eiv iv −iv e −e −eiv
−e−iv e−iv e−iv e−iv
and v = ωc LV . The structure of the periodic array shown in Fig. 4 always displays a degenerate band edge at k0 = π independent of input parameters, so we don’t need to find another candidates sets of k0 . As for the parameters of optimization, similar to magnetic frozen mode, we consider the thickness of the vacuum layer (LV ), the thickness of the first anisotropic layer(LA1 ) and the material constants of the anisotropic layer(εA , δA ). Figure 5 shows the sequence of convergence of the dispersion relation for a degenerate band edge point for optimization using one variable LV to recover the known solution, confirming the robustness of the method and the computed solution. In Table 2 we show other solutions found by using different parameters(n) for optimizations. We observe that the thickness of vacuum is approximately the same for all sets, but the thickness of anisotropic layer varies. The first derivative and the second derivative of all sets are less than 1.0e-10, which is sufficiently small to make the phenomena distinctive in the time domain. However, one big difference is that the transmittance rate in set 4 is significantly larger than those for the other sets in some contrast to the widely held belief that total reflectance or negligible amount of transmittance rate at degenerate band edge is inevitable. When
12
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.6
!
1
0.8
!
!
1.2
0.4
0.4
0.4
0.2
0.2
0.2
0
0
1
2
3
4
5
0
6
0
1
2
k
3
4
5
0
6
0
1
2
k
3
4
5
6
k
Figure 5: Scenes from pptimization process with one variable Vw(left to right). Vw=0.3029(left), 0.3994(cen2 ter), 0.4589 (right) after 6 iterations, dω = 2.6258e-015, ddxω2 = 1.3466e-013 with λ = 100. dx
it was believed so, as one of way to increase the transmittance rate, Figotin and Vitebsky suggested to use Fabry-P´erot peak cavity resonance in the paper [17]. In other words, instead of sending a monochromatic wave with the exact degenerate band edge frequency ω0 , we send a monochromatic wave with ω1 which is slightly less than ω0 , but with a transmittance rate which is significantly higher because of the finiteness of the slab. One of optimized solutions shows a transmittance rate of 0.3474 approximately, which is more than 50 times larger than those of the other sets. This result seems to be encouraging, since we don’t need to blur the distinctive features of the frozen mode by using ω1 different from ω0 . Fig. 6 and 7 show the field amplitude at T=300.0 and the increase of the maximum field amplitude with set 1 and set 4, respectively. These figures clearly display the significance of large transmittance rate at the degenerate band edge.
PSI
Max. of Field amplitude
2
1.6
1.8
1.4
1.6 1.2
Max. of Field amplitude
1.4 1.2 1 0.8 0.6
1 0.8 0.6 0.4
0.4 0.2
0.2 0
5
10
15 X
20
0
0
50
100
150
200
250
300
350
time
Figure 6: Field amplitude(left) and Maximum of field amplitude (right) of Degenerate band edge with set 1. Number of unit cell=32. Number of grid points per each domain=12. One unit cell = A1+V+A2 as shown in Figure 4, A1, V and A2 constitute each domain. E"x polarization.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
PSI
Max. of Field amplitude 20
18
18
16
16
14
14
Max. of Field amplitude
20
12 10 8 6 4
12 10 8 6 4
2 0
13
2 6
8
10
12
14
16
18
20
0
0
X
50
100
150
200
250
300
350
time
Figure 7: Field amplitude(left) and maximum of field amplitude (right) of degenerate band edge with parameter set 4. Number of unit cell=32. One unit cell = A1+V+A2 as shown in Figure 4, A1, V and A2 constitute each domain. Number of grid points per each domain=12. E"x polarization.
Input
Solution
n variables φ2 LV LA1 εA δA ω0 k0 dω dk |ω=ω0 d2 ω | dk 2 ω=ω0 TR rate
Set1 1 LV , LA1 = LA2 π/4 0.45890 0.27055 5.1 1.1 1.00940 π 2.9175d-16 6.0364d-14 0.0064
Set2 2 LV , LA1 π/4 0.46000 0.29190 5.1 1.1 1.01331 π 2.9063d-16 1.3414d-13 0.0070
Set3 2 LV , LA1 π/3 0.45005 0.43633 5.1 1.1 1.11726 π 0* 1.5102e-013 0.0440
Set4 4 LV , LA1 , εA , δA π/4 0.45005 0.45050 7.0 3.92 0.88645 π 4.3189d-15 3.8326d-10 0.3474
Table 2: Optimized solutions with different number of variables or different initial values for degenerate band edge. φ2 =misalignment angle of second anisotropic layer(we let φ1 = 0). LV = thickness of vacuum, LA1 =thickness of the first anisotropic layer. εA and δA are used in the permittivity of anisotropic layer as follows. ( * = the value is less than machine accuracy. )
εˆ =
εA + δA cos(2φ) δA sin(2φ) . TR rate = Transmittance rate δA sin(2φ) εA − δA cos(2φ)
14
5
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Optimization and Design of Oblique Frozen Mode Crystals
As proposed by Figotin and Vitebsky [5, 6], stationary inflection points can be also generated without magnetic layers and with structures even simpler as the degenerate band gap structures. To achieve this goal, one needs to break the symmetry and consider incident waves at oblique angles. When an oblique incident waves propagates into a periodic layer consisting of one anisotropic layer with εxz != 0 and vacuum, we can observe similar, but weaker frozen mode inside the layer as compared to the gyrotropic crystal. However, the crystal parameters in [5, 6] were computed by one dimensional transfer matrix method based on the following tangential consistency of the wave to eliminate the derivatives along x, y directions, i.e., 1 % & % &0 ˆ ˆ k k dE d E x,y x,y ˆ= ˆ= = ikx,y E iω E − dx ω ω dt The frozen mode was shown to exist by Chun and Hesthaven [2], after applying the above consistency into Maxwell’s equations. However, the existence of the oblique frozen mode in real three dimension remains in question since the above condition is artificial. To model this general case where the transfer matrix method is not applicable, we must solve the full Maxwell eigenvalue problems for the dispersion relation. The solution of this more complicated problem can be computed with a freely available frequency-domain solver [14]. In addition to choosing a wave vector k0 , for the stationary inflection point ω0 , we need to choose a band containing ω0 . For example, ω0 exists on the second band and the first band for the magnetic frozen mode and degenerate band edge respectively, but for the oblique frozen mode, no such information is known in advance. To identify the band where ω0 exist, we search whole bands for a wave vector k0 such as $ $ $' $ $ i+p dω $ $ (kj )$$ < ε, p ≥ 1, ε = small number D12 (k0 ) = min $ $j=i−p dk $
Figure 8: Array of crystals allowing an oblique frozen mode. The periodic structure consists of one anisotropic layer with εxz #= 0 and vacuum between them. The incident waves need to propagate into the layer at an oblique angle.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
1.38
1.32
1.36
1.31
1.29
V13
1.3
1.32
V13
1.34
15
1.3
1.28
1.28
1.27
1.26
1.26
1.24 -0.3
-0.2
-0.1
0
V4
0.1
0.2
0.3
1.25 -0.15
-0.1
-0.05
0
V4
0.05
0.1
0.15
Figure 9: Dispersion relation for the three-dimensional Maxwell problem with material parameters from paper [6] (left), and with the optimized solution for the 8th band(right). A stationary inflection point ω0 is observed in the dispersion relation only with the optimized solution
To increase the speed of the search process, we choose several bands in advance which seem to have k0 satisfying the above inequality. As for parameters for optimization, we use the oblique angle of incident wave(kx = ky ), the thickness of vacuum(LV ) and material constants of anisotropic layer(ε11 , ε33 ) with which material constant matrix is given as ε11 cos θ 2 + ε33 sin θ 2 0 (ε11 − ε33 ) cos θ sin θ , 0 ε22 0 εˆA = 2 2 (ε11 − ε33 ) cos θ sin θ 0 ε33 cos θ + ε11 sin θ
Among several stationary inflection points found by optimization, one of them with the lowest frequency ω0 is given as follows(also see Fig. 9): θ = π/4 ε11 = ε22 = 4.66812966, ε33 = 3.9656, LV = 0.5711535 $k = (−0.4800309, −0.4800309, −0.008982), ω0 = 1.295439 d2 ω dω = 4.82793e − 4, |ω=ω0 = 0.004393 dk ω=ω0 dk2 The dispersion relation resulting from the new parameters clearly show the correct properties for the three-dimensional crystal.
5.1
Time-Domain Modeling of the 3D Crystal
To observe the phenomena in three dimension, the new set of optimized parameters is solved in a time-domain wave propagation model with periodic boundary condition along x, y direction [20]. The periodic boundaries are x and y sides and absorbing boundary condition(ABC) is added at the end of z-direction to truncate unnecessary reflection of the outgoing waves(Fig. 10). Let only εxz be non-zero off-diagonal element for material
16
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Figure 10: Wave propagation problem with periodic boundary condition along x and y direction with ABC at the end of z direction. The incident wave propagates into the layer with an oblique angle
constant matrix.
˜z ∂E ∂y ˜x ∂E ∂z ˜y ∂E ∂x
−
−
−
˜y ∂E ∂z ˜z ∂E ∂x ˜x ∂E ∂y
˜ Hx ∂ =− H ˜y , ∂t ˜z H
˜z ∂H ∂y ˜x ∂H ∂z ˜y ∂H ∂x
−
−
−
˜y ∂H ∂z ˜z ∂H ∂x ˜x ∂H ∂y
˜ Ex εxx 0 εxz = 0 εyy 0 ∂ E ˜y ∂t ˜z εxz 0 εzz E
˜ H ˜ as two dimensional Fourier Assuming periodicity along x, y direction, we can express E, series as follows. ˜ y, z, t) = { E(x,
Ny /2
'
Nx /2
Ny /2
Nx /2
'
2π
2π imy L y
E(mx , my , z, t)eimx Lx x e
y
my =Ny /2 mx =Nx /2
˜ H(x, y, z, t) = {
'
'
2π
2π imy L y
H(mx , my , z, t)eimx Lx x e
my =Ny /2 mx =Nx /2
}
y
}
where, Nx , Ny are number of points along each x and y direction and Lx , Ly are length of domain along x and y direction, respectively. Note that Lx and Ly depend on the angle of incidence of plane wave. By substituting these expansions into the above Maxwell’s equations after implementing ABC, we have ∂E sz 0 (imy )Ez − ∂zy ∂Ex 0 sz = iω ∂z − (imx )Ez 0 0 (imx )Ey − (imy )Ex ∂Hy εxx (imy )Hz − ∂z ∂Hx 0 = −iω − (im )H x z ∂z ε (imx )Hy − (imy )Hx xz
0 Hx 0 Hy Hz 1/sz 0 εxz sz 0 0 Ex εyy 0 0 sz 0 Ey 0 εzz Ez 0 0 1/sz
where ω stands for the frequency of the incident plane wave.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
1.1
17
45 40
1.05
35 1
TR angle
TR rate
30 0.95
0.9
25 20 15
0.85 10 0.8
0.75
5
0
10
20
30
40 50 Incidence angle
60
70
80
0
0
10
20
30
40 50 Incidence angle
60
70
80
Figure 11: Results of testing the code for 3D time domain modeling with periodicity along x, y direction. Transmittance rate vs. incidence angle(left), transmittance angle vs. incidence angle(right). Solid line indicates the exact value and circles indicate computed values.
Then, by substituting sz = 1 + σz /(iω), we have the following Maxwell equation in time domain( ∀ mx ≤ |Nx /2|, ∀ my ≤ |Ny /2|, Dε = εxx εzz − ε2xz ) ∂Hx ∂t ∂Hy ∂t ∂Bz ∂t ∂Hz ∂t ∂Ex Dε ∂t ∂Ey εyy ∂t ∂Dz Dε ∂t ∂Ez ∂t
∂Ey − (imy )Ez − σz Hx ∂z ∂Ex = − + (imx )Ez − σz Hy ∂z =
= −(imx )Ey + (imy )Ex = = = = =
∂Bz + σz Bz ∂t ∂Hy −εzz { − (imy )Hz } + εxz {(imy )Hx − (imx )Hy } − σz Ex ∂z ∂Hx − (imx )Hz − σz Ey ∂z ∂Hy εxz { − (imy )Hz } − εxx {(imy )Hx − (imx )Hy } ∂z ∂Dz + σz Dz dt
A Discontinuous Galerkin Method is used for differentiation along z-direction as in the paper [18, 19]. As for time marching, we use an implicit scheme for the linear derivatives along z and an explicit scheme for what remains [15, 16]. To validate the code, the transmittance rate and transmittance angle are computed when a plane incident wave propagating from vacuum to a dielectric material with various oblique angles. As shown in Fig. 11, it demonstrates the excellent performance of the scheme In Fig. 12, we show the energy density flux and the field amplitude with the optimized parameter values and observe several features of the frozen mode. In behalf of analysis,
18
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
2.2 2
2
0.4
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
fampz
EDFz
0.3
0.2
0.1
0
2.2
fampz
0.5
1 0.8
0.6
0.4
0.4
0.2
0.2
0 -0.1
0
50
z
0
-0.2 -50
100
1 0.8
0.6
0
50
z
100
150
200
-0.2 -50
0
50
z
100
150
200
Figure 12: Energy density flux(left) and the field amplitude with optimized parameters at T=100(center) and at T=200(right), Nx = Ny = 6, Nz = 12, Nx , Ny , Nz =number of grid points along x,y,z. Number of unit cell=100. Unit cell = A+V
0.05
0.5 0 0.4
0
0.3 -0.05
V4
V3
V2
-0.05 0.2
-0.1 0.1 -0.1 -0.15
-0.2 -100
0
0
100
V1
200
300
-100
0
100
V1
200
300
-0.1 -100
0
100
V1
200
300
Figure 13: Energy Flux Density Difference(EDFD) between the first cell and the last cell along x(left), y(center), z(right) direction. Solid line is with optimized parameters, dashed line is with original values. Nx = Ny = 6, Nz = 12, Number of unit cell=64, Unit cell = A+V
let’s define energy density flux difference(EDFD) as shown in [2] EDFD =
2
S
ef · nds =
2
V
∇ef dV ,
ef = energy density flux
note that this EDFD indicates the amount of energy density stacked inside the layer per unit time. Then, we observe that there are significant amounts of EDFD along the z direction between the first cell and the last cell. This leads to an increased field amplitude, compresses the energy density and the slow-down of waves. However, each of these phenomena are weaker than that of magnetic frozen mode as expected [2]. The maximum of the field amplitude is increasing very slowly and the average velocity of the energy transfer is approximately 0.4c. Fig. 13 compares EDFD for all directions for the crystals based on the original values and the optimized values. With the optimized values, we observe that the EDFD along x, y direction are almost zero, but EDFD along the z direction remains significantly large even after a long time. With original values, EDFD along the x, y direction remains nontrivial.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
6
19
Conclusion
We have explored the use of PDE constrained optimization along with design and optimization of photonic crystals exhibiting frozen mode behavior. This phenomenon is related to the existence of Van Hove singularities in the dispersion relation. We have demonstrated the ability to recover known solutions as well as computed several new sets of parameters by optimizing the geometry and/or the materials of the crystals. For both of the magnetic frozen mode and degenerate band edge, each sets of parameters obtained by optimization with different number of variables show the generally similar phenomena, but some features are different. Above all, the transmittance rate varies significantly, which is a crucial element for energy density and field amplitude inside the layer. In particular one set for the degenerate band edge displays a transmittance rate around 0.3475, which is more than 50 times larger than those of the others. For the case with the frozen mode requiring an oblique illumination, we succeed in finding a stationary inflection point in three dimensions with parameters’ sets derived by nonlinear optimization. The existence of the frozen mode in three dimension is shown for the first time, but the phenomena themselves are not as distinctive as the previous two instances. One of the reason can be the intrinsic property of the oblique array which makes wave vector k0 too close to zero. The current computational design tool allows one to specify materials and/or geometries, then seek to design crystals with the sought-after properties, if they exist. The computations shown here illustrate, however, that the Van Hove singularities are found broadly in the dispersion relations, indicating that many configurations often exist. Future works may include finding stationary inflection points or degenerate band edge points showing strong frozen mode in 2 or 3 dimensional structures. We also believe that using a more complex geometry, but employing simpler materials, can be a candidate for exploring similar behavior.
Acknowledgements This work was supported by the U.S. Air Force Office of Scientific Research under the grant FA9550-04-1-0359. References [1] J.S. Hesthaven and T. Warburton, Discontinuous Galerkin Methods for the Time-Domain Maxwell’s Equations: An Introduction, ACES Newsletter, 19(1), 10-29, 2004. [2] Sehun Chun and Jan S. Hesthaven, Modeling of the frozen mode phenomenon and its Sensitivity using Discontinuous Galerkin methods, Communications in Computational Physics, to appear, 2006. [3] A. Figotin and I. Vitebsky, Nonreciprocal magnetic photonic crystals, Physical Review E, 63, 0666609, May, 2001.
20
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
[4] A. Figotin and I. Vitebsky, Electromanetic unidirectionality in magnetic photonic crystals, Physical Review B, 67, 165210, April, 2003. [5] A. Figotin and I. Vitebsky, Oblique frozen modes in periodic layered media, Physical Review E, 68, 036609, September, 2003. [6] J. Ballato and A. Ballato and A. Figotin and I. Vitebsky, Frozen light in periodic stacks of anisotropic layers, Physical Review E, 71, 1, 2005. [7] L. Van Hove, Van Hove singularities, Phys. Rev, 89, 1189, 1953. [8] M. Ibanescu and S.G. Johnson and D. Roundy and C. Luo and Y.Fink and J.D. Joannopoulos, Anomalous Disperion Relations by Symmetry Breaking in Axially Uniform Waveguides, Physical review letters, 92, 6, 2004. [9] K. Svanberg, The method of moving asymptotes: a new method for structural optimization, Int. J. Nuner. Meth. Engng, 24, 359-373, 1987. [10] C. Zillober, A globally convergent version of the method of moving asymptotes, Structural Optimization, 6, pp. 166-174, 1993. [11] W.H. Zhang and C. Fleury and P. Duysinx and V.H. Nguyen and I. Laschet, A generalized method of moving asymptotes(GMMA) including equality constraints, Structural Optimization, 12, pp. =143-146, 1996. [12] K. Svandberg, A new globally convergent version of MMA, Technical Report TRITA-MAT1999-OS2, Dept. of Mathematics, KTH, Stockholm, 1999. [13] K. Svanberg, Two primal-dual interior-point methods for the MMA subproblems. Technical Report TRITA-MAT-1998-OS12, Dept. of Mathematics, KTH, Stockholm, 1998. [14] Johnson, Steven G. and Joannopoulos, J. D., Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis, Opt. Express, 2001, 8, 3, pp. 173-190, url =http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173. [15] Uri M. Ascher and Steven J. Ruuth and Raymond J. Spiteri, Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics, 1997, 25, pp. 151-167. [16] Christopher A. Kennedy and Mark H. Carpenter, Additive Runge-Kutta schemes for convection-diffusion-reaction equations, Applied Numerical Mathematics, 2003, 44, pp. 139181. [17] Alex Figotin and Ilya Vitebskiy, Gigantic transmission and band-edge resonance in periodic stacks of anisotropic layers, Physical review E, 2005, 72, 036619. [18] J.S. Hesthaven and T. Warburton, Discontinuous Galerkin Methods for the Time-Domain Maxwell’s Equations: An Introdunction, ACES Newsletter, 19(1), pp. 10-29, 2004. [19] J.S. Hesthaven and T. Warburton, Nodal High-Order Methods on Unstructured Grids, Journal of Computational Physics, Journal of Computational Physics, 182, pp. 186-221, 2002. [20] Allen Taflove and Susan C. Hagness, Computational Electrodynamics, the Finite-Difference Time-Domain Method, 3rd edition, Artech House, Boston, 2005. [21] Uri M. Ascher and Linda R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, siam, Philadelpha, 1998. [22] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis 2nd edition, Springer-Verlag, New York, 1979.
Commun. Comput. Phys. December 2006
PDE Constrained Optimization and Design of Frozen Mode Crystals S. Chun1 , J. S. Hesthaven1 1
Division of Applied Mathematics, Brown University, Providence, RI 02912
Abstract. We explore the use of PDE constrained nonlinear optimization techniques to optimize and design electromagnetic crystals which exhibit frozen mode behavior. This is characterized by the association with Van Hove singularities in the dispersion relation, e.g., stationary reflection points and degenerate band edge points. Hence, the optimization process modifies the dispersion relation by adjusting the geometries and material parameters. The resulting algorithm is found to be capable of recovering all known crystal configurations as well as many new configurations, some of which display dramatically improved properties over previously used configuration. We investigate both gyrotropic photonic crystals and degenerate band edge crystals as well as the more complex case of the oblique incidence where we extend the investigation to the three-dimensional case to identify the first three-dimensional crystal exhibiting frozen mode behavior. Key words: PDE constrained optimization, Method of Matched Asymptotics, Frozen mode, Photonic crystal
1
Introduction
Controlling the flow of electromagnetic waves in materials has been the frontier of our technologies in the last decade. Along with prohibiting wave propagation or localizing wave in a certain area, slowing down wave at the significant rate has grabbed broad attentions. Recent studies of periodic structures of several different anisotripic materials shows that when a monochromatic wave with a certain frequency propagates in such structures, it shows abnormal electromagnetic phenomena such as dramatic slow-down of the wave, the significant increase of field amplitude and singularity of the trasmittance rate, which are classified as the frozen mode. These unique features of the frozen mode are attractive for numerous practical applications and recent studies [2] have confirmed
∗
Correspondence to: Division of Applied Mathematics, Brown University, Box F, Providence, RI 02912. Email: [email protected]; [email protected] http://www.global-sci.com/
1
c !2006 Global-Science Press
2
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
the robustness of these phonomena. In the beginning, the frozen mode was regarded as a distinctive phenomenon only related to stationary inflection points of the dispersion relations ω(k) such as, ∂ω = 0; ∂k
∂2ω = 0; ∂k2
∂3ω != 0 ∂k3
at ω = ω0
However, it is worth mentioning that these properties of the dispersion relation are special cases of what is known as Van Hove singularities [7] as points where the first derivative of the dispersion relation ω(k) vanishes, i.e., ∂ω = 0, ∂k
at ω = ω0
This nomenclature is related to density of state(DOS), since it is well-known that DOS is inversely proportional to ∂ω/∂k [8]. However, all previous studies have focused on the analysis of crystal configurations found by a trial-and-error technique which tends to limit the number of parameters one can freely vary. This limits the flexibility in the design of such crystals and leaves open the question of whether it is generally possible to achieve the sought-after properties in the dispersion relation, e.g., when given some specific materials, one can modify the geometric properties of the crystal as needed. To address this problem in a more systematic way, we consider the use of a PDE constrained optimization approach where we use the dispersion relation, given by the transfer matrix method or by solving the Maxwell eigenproblem, in combination with the algorithm of Method of Moving Asymptote(MMA) [9–13] to recreate and optimize known configurations as well as to design entire crystal arrangements. In this study, we consider the optimization and the design of three different types of crystals, all exhibiting the frozen mode phenomena. The first one is a gyrotropic photonic crystal consisting of two misaligned anisotropic layers and one magnetic layer. Its dispersion relation has a stationary inflection point ω0 on the 2nd band. The second class of crystals consists of two misaligned anisotropic layers and vacuum between them. Its dispersion relation has a degenerate band edge ωd in the first band. The final structure consists of only one anisotropic layer with εxz != 0 and vacuum between them. However, the frozen mode phenomenon only emerges when the incident wave propagates into this layer at an oblique angle. This latter case is quite different from the previous two cases. The dispersion relation where a stationary inflection point ω0 lies was originally computed by a one-dimensional transfer matrix method using tangential consistency of wave vector. To verify the existence of a stationary inflection point in general three dimensional setting which was previously unascertained, we first solve the general Maxwell eigenproblem to get the dispersion relation with initial parameter sets and apply the optimization techniques to seek a stationary inflection point in the dispersion relation of three dimension. The remaining part of the paper is organized as follows. In Sec. 2, we briefly recall the physical model, Maxwell’s equations and how one recovers the dispersion relation from
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
3
these equations. This sets the stage for Sec. 2 where we explain the implementation of the optimization method based on Method of Moving Asymptote(MMA). In Secs 3 and 4, we display the methods and results of optimizing the magnetic frozen mode and degenerate band edge for known designs as well as new configurations, respectively. Section 5 is devoted to the design of crystals for the oblique case and its modeling in the time domain. We conclude with a few remarks and future works in Sec. 6.
2
Physical Model
We consider Maxwell’s equations ˜ =− ∇×E
˜ dB , dt
˜ = ∇×H
˜ dD dt
˜ and H ˜ is the electric field and magnetic field, respectively, D ˜ and B ˜ the displacewhere E ment field and the induction field. These fields are related by linear constitutive relations as follows. ˜ ˜ D(z) = εˆ(z)E(z),
˜ ˜ B(z) =µ ˆ(z)H(z)
where we have assumed that (ˆ ε, µ ˆ) are homogeneous in the (x, y)-plane. If we consider harmonic solutions with frequency ω, we recover the Maxwell’s equations in the frequency domain. ˜ = iω µ ˜ . ∇×H ˜ = −iω εˆ(z)E ˜ ∇×E ˆ(z)H
(2.1)
To simplify the above equations, we introduce the normalized quantities using z=
z˜ , ˜ L
t=
t˜ ˜ c˜0 L/
˜ is a reference length and c˜0 = (˜ where, L ε0 µ ˜0 )− 2 is the dimensional speed of light in ˜ and H ˜ are normalized as vacuum. Also, E ! ˜ ˜ ε˜0 E H , H= E= ˜0 ˜0 µ˜0 H H 1
˜ 0 is a dimensional reference magnetic field strength. By substituting these norwhere H malized quantities into Eq.(2.1), we obtain the following non-dimensionalized equations ∇ × E = iω µ ˆ(z)H . ∇ × H = −iω εˆ(z)E , where (ˆ ε, µ ˆ) now represent the relative permittivity and permeability, respectively.
4
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
If we now assume that the solution is periodic and seek solutions of the form # " # " $ x) $k E($ E $ $k = wave number $ k exp(ik · $x) , $ x) = H H($ we recover the modal equation ˆ(z)Hk , (∇ + i$k) × Ek = iω µ
(∇ + i$k) × Hk = −iω εˆ(z)Ek ,
which simplifies as (∇ + i$k) ×
1 (∇ + i$k) × Ek = ω 2 εˆ(z)Ek µ ˆ(z)
For a given value of the wave vector $k, this is recognized as an eigenvalue problem for (ω 2 , Ek ), from which we recover the dispersion relation ($k, ω($k)). Hence, when optimizing the crystals, we must ensure that all considered solutions satisfy the above eigenvalue problem which emerges as the constraints in the optimization process. The goal of the optimization is to identify one type of Van Hove singularities [7] in the dispersion relation, leading to a very strong divergence of the density of states(DOS) near the corresponding frequency. There are several types of singularities, but we are specially interested in finding singularities where the first and second derivative are both zero, i.e. ∂ω = 0, ∂k
∂2ω = 0, ∂k2
at ω = ω0
We state this equivalently as the following minimization problem. $ $ $ 2 $ $ dω $ $d ω $ $ $ $ minimize λ $ ($x)$ + $ 2 ($x)$$ , at $k = $k0 dk dk ≤ xj ≤ xmax , xmin j j
(2.2)
1≤j≤n
where, $x is a n-variable vector and λ is a weight constant. The objective function dω x) dk ($ 2 and ddkω2 ($x) can be obtained directly from the dispersion relation with optimization variable $x by solving the eigenvalue problem discussed in Sec. 2. Once the dispersion relation is obtained, the derivatives of the objective functions can be computed easily. For example, the fourth order approximation of its first and second derivatives are computed by difference approximations such as, & & % % 1 dωj+1 dωj−1 dωj−2 ∂ dω dωj+2 ($x) = +8 −8 + − ∂xj dk 12∆xj dk dk dk dk & & % % 2 dωj+1 dωj dωj−1 dωj−2 dωj+2 1 dω ∂ ($x) = + 16 − 30 + 16 − − dk dk dk dk dk ∂x2j dk 12∆x2j , where
dω dωj = ($x + j∆$x), dk dk
1≤j≤n
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
5
In fact, this fourth order approximation is preferred in the computation of the derivatives to increase the accuracy and to make the fast convergence. Furthermore, it is required to obtain an initial guess $k0 or at least a group of possible candidates of $k0 s. Identifying suitable $k0 is one of most critical and significant factors of convergence of the optimization solution. We have identified the following approach to work well. 1. Divide the dispersion relation into several sections such that each section contains $ only one wave vector with | dω dk (k)| < ε & 1. 2. For each section, find one wave vector $k0 such that D12 ($k0 ) = min | i
i+2 ' dω $ (k" )|, dk
"=i−2
i ∈ section I
Note that the function inside the bars is a way of indicating the amount of the first and the second derivative at the same time 3. Optimize the objective function at $k0 . Since the objective function is generally not convex, we use the algorithm of Method of Moving Asymptote(MMA). The application of MMA to find singularity points of dispersion relation can be stated briefly as follows. First, Eq.(2.2) can be changed into the following equivalent problem. minimize
f0 (x) = z +
4 '
(2.3)
ci yi
i=1
subject to
dω d2 ω ($x) − y1 ≤ 0, f2 (x) = ($x) − y2 ≤ 0 dk dk2 dω d2 ω f3 (x) = −λ ($x) − y3 ≤ 0, f4 (x) = − 2 ($x) − y4 ≤ 0 dk dk z ≥ 0, yi ≥ 0, 1 ≤ i ≤ 4, f1 (x) = λ
≤ xj ≤ xmax , xmin j j
1≤j≤n
where, y1 , · · · , yn and z are artificial optimization variables. Note that at the optimized solution $x, z is zero and yi is given as the maximum of each constraint function. Now, we have the non-linear optimization problem with several constraints. For these fi ($x)s, let’s define gi (x) in the following way, gi (x) = ri +
n ' j=1
(
pij qij + ) Uj − xj xj − L j
pij = max{0, (Uj − xj )2 ri = fi (x) −
n ' ( j=1
∂fi }, ∂xj
(2.4)
qij = max{0, −(xj − Lj )2
pij qij + ), Uj − xj xj − L j
1 ≤ i ≤ 4,
∂fi } ∂xj
1≤j≤n
6
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
The parameters Lj and Uj (known as moving asymptotes) are chosen as Lj < xj < Uj and play a critical role in the convergence of the optimization scheme. In fact, Lj and Uj are adjusted for each iteration according to the convergence of the solution. Note that the second derivatives of gi are given as % & 2pij 2qij ∂ 2 gi = δj" + , 1 ≤ j, ' ≤ n ∂xj x" (Uj − xj )3 (xj − Lj )3 Since pij ≥ 0, qij ≥ 0, gi is a convex function. In fact, as long as Lj and Uj are finite, ∂gi gi is strictly convex in all variables except for ∂x (x) = 0, 1 ≤ j ≤ n. With gi s as in the j Eq.(2.4), we have the following subproblem. n % '
& qij pij + + r0 Uj − xj xj − L j j=1 % & pij qij + + ri ≤ fi , Uj − xj xj − L j
P : minimize subject to
≤ xj ≤ xmax , xmin j j
(2.5) 1≤i≤4
1≤j≤n
Then, the Lagrangian function corresponding to this problem is given by y T $r + W (x, y) = $r0 − $
& n % ' p0j + $y T p$j q0j + $y T $qj + Uj − xj ($y ) xj ($y ) − Lj j=1
In fact, W (x, y) is differentiable, smooth and concave. Finally, the Lagrangian dual problem with this Lagrangian function is as follows. D : maximize
& n % ' q0j + $y T $qj p0j + $y T $pj + W (x, y) = $r0 − $y $r + Uj − xj ($y ) xj ($y ) − Lj T
j=1
subject to
y≥0
This yields a nice convex optimization problem and it can be solved by a primal-dual Newton method. The details for this primal-dual method are given in [13]
3
Optimization and Design of Magnetic Crystals
The frozen mode in the gyrotropic crystal occurs when a monochromatic wave with frequency ω0 propagates into a periodic array consisting of two misaligned anisotropic layers and one magnetic layer(Fig. 1). In the dispersion relation, this ω0 is a stationary inflection point with the following properties ∂ω = 0, ∂k
∂2ω = 0, ∂k2
∂3ω != 0, ∂k3
at ω = ω0
(3.1)
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
7
As shown in [4], one set of parameters parameters resulting in an inflection point at ω0 is
For A layer : εˆ =
(
)
13.61 + 12.40 cos(2φ) 12.40 sin(2φ) 12.40 sin(2φ) 13.61 − 12.40 cos(2φ)
µ ˆ=I
,
(
) ( ) 5.0 0.0 60.0 i37.0 For F layer : εˆ = , µ ˆ= 0.0 5.0 −i37.0 60.0 φ1 = 0, φ2 = −π/4, LF = 0.0047454, LA = 0.5(1.0 − LF ) k0 = 2.6329,
ω0 = 0.607676756(c/L)
where φ1 , φ2 is the misalignment angle of the first and second anisotropic layer. Also LF is the thickness of F layer and LA1 , LA2 are those of the first and second anisotropic layer, respectively. The first test is to recover this set of parameters by the aforementioned optimization technique and, subsequently, seek sets of parameters if possible. Since the periodicity is along z direction only, the dispersion relation of the magnetic frozen mode can be obtained by the transfer matrix method rather than solving the full Maxwell eigenvalue problem. The transfer matrix method is a way of finding the dispersion relation from Bloch solutions of Maxwell’s equations in planar structures. The transfer matrix(TL (ω)) can be obtained with bases of tangential components of electrical and magnetic fields. Then, it is used for Bloch solutions of the Maxwell’s equations in a periodic media such as qk (z + L) = eikL qk (z),
qk = [Ex Eu Hx Hy ]Tk
Using the property of transfer matrix TL (ω)qk (z) = qk (z + L), we get the following charicteristic equation. (TL (ω) − eikL )qk (L) = 0 =⇒
det(TL (ω) − eikL ) = 0
0.8 0.7 w0 0.6 0.5 0.4 0.3 0.2 0.1 0
0
1
2
3
4
5
6
Figure 1: Periodic structure of gyrotropic photonic crystals consisting of two different anisotropic materials with one magnetic layer(left) and its dispersion relation(right).
8
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Dispersion relation when rho0= 0.018924
Dispersion relation when rho0= 0.0095342 0.7
0.6
0.6
0.6
0.5
0.5
0.5
0.4
0.4
0.4
0.3
freqency
0.7
freqency
freqency
Dispersion relation when rho0= 0.05 0.7
0.3
0.3
0.2
0.2
0.2
0.1
0.1
0.1
0
0
1
2
3 4 wave vector
5
6
0
7
0
1
2
3 4 wave vector
5
6
7
0
0
1
2
3 4 wave vector
5
6
7
Figure 2: Scenes frm optimization process with one variable LF (left to right). LF =0.0244(left), 0.009453(cen2 ter), 0.00474(right) after 6 iterations, dω = 3.97d-13, ddxω2 = 0.007907 with λ = 10000. dx
which yields the dispersion relation between ω and the wave vector k. The transfer matrix of an anisotropic layer with misalignment angle φ (TA ) and a magnetic layer (TF ) are given as TA (φ, LA ) = W (φ, LA )W −1 (φ, 0), TF (LF ) = W (LF )W −1 (0) (cos φ)ein1 a (cos φ)e−in1 a −(sin φ)ein2 a −(sin φ)e−in2 a in a −in a in a 1 1 2 (sin φ)e (sin φ)e (cos φ)e (cos φ)e−in2 a W (φ, LA ) = −η1 (sin φ)ein1 a η1 (sin φ)e−in1 a −η2 (cos φ)ein2 a η2 (cos φ)e−in2 a η1 (cos φ)ein1 a −η1 (cos φ)e−in1 a −η2 (sin φ)ein2 a η2 (sin φ)e−in2 a ein1 f e−in1 f −iein2 f −ie−in2 f −iein1 f −ie−in1 f ein2 f e−in2 f W (LF ) = iη1 ein1 f −iηe−in1 f −η2 ein2 f η2 e−in2 f η1 ein1 f −ηe−in1 f −iη2 ein2 f iη2 e−in2 f
where, a =
ω c LA
and f =
ω c LF .
Then, the transfer matrix of three-cell layer is
TL (φ1, φ2, LA1 , LA2 , LF ) = TF (LF )TA (φ2 , LA2 )TA (φ1 , LA1 ) and the dispersion relation of this layer is obtained by solving det(TL − ξI) = 0,
ξ = eikL , L = LA1 + LA2 + LF .
The left hand side is a fourth order polynomial and we can compute the roots of it by simple root-finding. As for choosing wave vector k0 corresponding to ω0 , we first find the wave vector kb corresponding to a band edge ωb (the highest peak such that dω dk |ω=ωb = 0) and seek k less than kb where the minimum of D12 occurs as discussed above. However, this may still miss the true candidates of wave vector k0 for optimized solution ω0 . To address this issue, we
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Input
Solution
n variables φ2 LF LA1 εA δA ω0 k0 dω dk |ω=ω0 d2 ω | dk 2 ω=ω0 TR rate
Set1 1 LF , (LA1 = LA2 ) −π/4 0.00474 0.49763 5.1 1.1 0.60769 2.63894 3.9666d-13 0.007907 0.3388
Set2 2 LF , LA1 −π/4 0.00480 0.44555 5.1 1.1 0.60950 2.63894 2.41590d-16 0.00898 0.4081
Set3 2 LF , LA1 −π/3 0.00490 0.22056 5.1 1.1 0.67177 2.70177 6.5758d-16 0.00917 0.7786
9
Set4 4 LF , LA1 , εA , δA −π/4 0.00470 0.45001 5.0 1.18189 0.62171 2.63894 1.80310d-12 0.00191 0.3109
Table 1: Optimized solutions with different number of variables or different initial values for magnetic frozen mode. φ2 =misalignment angle of second anisotropic layer(we let φ1 = 0). LF = thickness of magnetic layer, LA1 =thickness of the first anisotropic layer. εA and δA are used in the permittivity of anisotropic layer as follows. ( )
εˆ =
εA + δA cos(2φ) δA sin(2φ) , TR rate = Transmittance rate δA sin(2φ) εA − δA cos(2φ)
change the initial values if the converged solution in the sense of *x − xold * < ε0 does not satisfy dω d2 ω < ε1 , < ε2 . dk dk2 In general, we shall allow the thickness of the magnetic layer(LF ), the thickness of the first anisotropic layer (LA1 ) and material constants of the anisotropic layer(εA , δA ) to vary during the optimization, although this could be extended as needed. In Fig. 2, we show the sequence of convergence of the dispersion relation for a stationary inflection point using optimization with just one variable, LF . We note that the optimized solution for LF is accurate up to 5.0e-06 with the analytically derived LF [4] and the difference is negligible in the sense of sensitivity. Other optimized solutions with different number of variables(n) are shown in Table 1. We observe that the thickness of F layer is approximately the same for all sets, but the thickness of A layer varies depending on the misalignment angle or the material constants. The absolute values of the first and second derivative at each ω0 are sufficiently close to zero, so all sets show every characteristics of the frozen mode. However, the closer to zero the first derivative is, the slower the group velocity(velocity of energy transfer) of the wave is. Also, the intensity of the field amplitude inside the layer and the amount of energy density stacked inside the layer indicated by energy density flux difference(EDFD) [2] are closely related to the transmittance rate. If a larger amount of energy enters into layer, the growth of the field amplitude gets higher and energy density
10
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
EDFD 0.45
90
0.4
80
0.35
70
0.3
60
0.25
EDFD
Max. of Field amplitude
Max. of Field amplitude 100
50 40
0.2 0.15
30
0.1
20
0.05
10 0
Set1 Set2 Set3 Set4
0 0
100
200
300
400
500
!0.05
0
100
200
time
300
400
500
time
Figure 3: Maximum of the field amplitude(left) and the Energy Density flux difference(ef of the first cell - ef of the last cell, ef =energy density flux) (right) for set 1(solid line), set 2(dashed line), set 3(dot-dashed line) and set 4(dotted line) of Table 1. Number of cell = 16. E"x polarization.
inside the layer gets stronger (Fig. 3).
4
Optimization and Design of Degenerate Band Edge Crystals
Different from the stationary inflection points in the gyrotropic crystals, the third derivative of the dispersion relation also needs to be zero in the degenerate band edge crystals, i.e. ∂ω = 0, ∂k
∂2ω = 0, ∂k2
∂3ω =0 ∂k3
at ω = ω0
According to the paper [17], the following set of parameters generate a degenerate band
1.2 w0 1
0.8
0.6
0.4
0.2
0
0
1
2
3
4
5
6
Figure 4: Periodic crystal structure allowing a degenerate band edge. The crystal consists of two misaligned anisotropic layers with vacuum between them(left) and its dispersion relation(right)
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
11
edge point at k0 = π. (
13.61 + 12.40 cos(2φ) 12.40 sin(2φ) For A layer : εˆ = 12.40 sin(2φ) 13.61 − 12.40 cos(2φ) φ1 = 0, φ2 = π/4, LV = 0.45891, LA = 0.5(1.0 − LV ) k0 = π,
)
,
µ ˆ=I
ω0 = 1.0094052895(c/L)
,where LA and LV are the thiness of anisotropic layer and vacuum, respectively. One of main advantages of this configuration is the lack of a magnetic components which are both complicated and difficult to manufacture with a sufficiently low loss. As shown in Fig. 4, the dispersion relation is indeed flat at the degenerate band edge point and it increases the speed of convergence when minimizing $ 2 $ $ $ $d ω $ $ dω $ J = λ $$ $$ + $$ 2 $$ dk dk As for the magnetic case, the one-dimensional nature of the problem enables the use of the transfer matrix method, which is expressed as for the three cells. TL (φ1, φ2, LA1 , LA2 , LV ) = TA (φ2 , LA2 )TV (LV )TA (φ1 , LA1 ) where
eiv e−iv −eivf iv e e−iv eiv TV (LV ) = W (LV )W −1 (0), W (LV ) = −eiv e−iv −eiv iv −iv e −e −eiv
−e−iv e−iv e−iv e−iv
and v = ωc LV . The structure of the periodic array shown in Fig. 4 always displays a degenerate band edge at k0 = π independent of input parameters, so we don’t need to find another candidates sets of k0 . As for the parameters of optimization, similar to magnetic frozen mode, we consider the thickness of the vacuum layer (LV ), the thickness of the first anisotropic layer(LA1 ) and the material constants of the anisotropic layer(εA , δA ). Figure 5 shows the sequence of convergence of the dispersion relation for a degenerate band edge point for optimization using one variable LV to recover the known solution, confirming the robustness of the method and the computed solution. In Table 2 we show other solutions found by using different parameters(n) for optimizations. We observe that the thickness of vacuum is approximately the same for all sets, but the thickness of anisotropic layer varies. The first derivative and the second derivative of all sets are less than 1.0e-10, which is sufficiently small to make the phenomena distinctive in the time domain. However, one big difference is that the transmittance rate in set 4 is significantly larger than those for the other sets in some contrast to the widely held belief that total reflectance or negligible amount of transmittance rate at degenerate band edge is inevitable. When
12
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
1.2
1.2
1
1
0.8
0.8
0.6
0.6
0.6
!
1
0.8
!
!
1.2
0.4
0.4
0.4
0.2
0.2
0.2
0
0
1
2
3
4
5
0
6
0
1
2
k
3
4
5
0
6
0
1
2
k
3
4
5
6
k
Figure 5: Scenes from pptimization process with one variable Vw(left to right). Vw=0.3029(left), 0.3994(cen2 ter), 0.4589 (right) after 6 iterations, dω = 2.6258e-015, ddxω2 = 1.3466e-013 with λ = 100. dx
it was believed so, as one of way to increase the transmittance rate, Figotin and Vitebsky suggested to use Fabry-P´erot peak cavity resonance in the paper [17]. In other words, instead of sending a monochromatic wave with the exact degenerate band edge frequency ω0 , we send a monochromatic wave with ω1 which is slightly less than ω0 , but with a transmittance rate which is significantly higher because of the finiteness of the slab. One of optimized solutions shows a transmittance rate of 0.3474 approximately, which is more than 50 times larger than those of the other sets. This result seems to be encouraging, since we don’t need to blur the distinctive features of the frozen mode by using ω1 different from ω0 . Fig. 6 and 7 show the field amplitude at T=300.0 and the increase of the maximum field amplitude with set 1 and set 4, respectively. These figures clearly display the significance of large transmittance rate at the degenerate band edge.
PSI
Max. of Field amplitude
2
1.6
1.8
1.4
1.6 1.2
Max. of Field amplitude
1.4 1.2 1 0.8 0.6
1 0.8 0.6 0.4
0.4 0.2
0.2 0
5
10
15 X
20
0
0
50
100
150
200
250
300
350
time
Figure 6: Field amplitude(left) and Maximum of field amplitude (right) of Degenerate band edge with set 1. Number of unit cell=32. Number of grid points per each domain=12. One unit cell = A1+V+A2 as shown in Figure 4, A1, V and A2 constitute each domain. E"x polarization.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
PSI
Max. of Field amplitude 20
18
18
16
16
14
14
Max. of Field amplitude
20
12 10 8 6 4
12 10 8 6 4
2 0
13
2 6
8
10
12
14
16
18
20
0
0
X
50
100
150
200
250
300
350
time
Figure 7: Field amplitude(left) and maximum of field amplitude (right) of degenerate band edge with parameter set 4. Number of unit cell=32. One unit cell = A1+V+A2 as shown in Figure 4, A1, V and A2 constitute each domain. Number of grid points per each domain=12. E"x polarization.
Input
Solution
n variables φ2 LV LA1 εA δA ω0 k0 dω dk |ω=ω0 d2 ω | dk 2 ω=ω0 TR rate
Set1 1 LV , LA1 = LA2 π/4 0.45890 0.27055 5.1 1.1 1.00940 π 2.9175d-16 6.0364d-14 0.0064
Set2 2 LV , LA1 π/4 0.46000 0.29190 5.1 1.1 1.01331 π 2.9063d-16 1.3414d-13 0.0070
Set3 2 LV , LA1 π/3 0.45005 0.43633 5.1 1.1 1.11726 π 0* 1.5102e-013 0.0440
Set4 4 LV , LA1 , εA , δA π/4 0.45005 0.45050 7.0 3.92 0.88645 π 4.3189d-15 3.8326d-10 0.3474
Table 2: Optimized solutions with different number of variables or different initial values for degenerate band edge. φ2 =misalignment angle of second anisotropic layer(we let φ1 = 0). LV = thickness of vacuum, LA1 =thickness of the first anisotropic layer. εA and δA are used in the permittivity of anisotropic layer as follows. ( * = the value is less than machine accuracy. )
εˆ =
εA + δA cos(2φ) δA sin(2φ) . TR rate = Transmittance rate δA sin(2φ) εA − δA cos(2φ)
14
5
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Optimization and Design of Oblique Frozen Mode Crystals
As proposed by Figotin and Vitebsky [5, 6], stationary inflection points can be also generated without magnetic layers and with structures even simpler as the degenerate band gap structures. To achieve this goal, one needs to break the symmetry and consider incident waves at oblique angles. When an oblique incident waves propagates into a periodic layer consisting of one anisotropic layer with εxz != 0 and vacuum, we can observe similar, but weaker frozen mode inside the layer as compared to the gyrotropic crystal. However, the crystal parameters in [5, 6] were computed by one dimensional transfer matrix method based on the following tangential consistency of the wave to eliminate the derivatives along x, y directions, i.e., 1 % & % &0 ˆ ˆ k k dE d E x,y x,y ˆ= ˆ= = ikx,y E iω E − dx ω ω dt The frozen mode was shown to exist by Chun and Hesthaven [2], after applying the above consistency into Maxwell’s equations. However, the existence of the oblique frozen mode in real three dimension remains in question since the above condition is artificial. To model this general case where the transfer matrix method is not applicable, we must solve the full Maxwell eigenvalue problems for the dispersion relation. The solution of this more complicated problem can be computed with a freely available frequency-domain solver [14]. In addition to choosing a wave vector k0 , for the stationary inflection point ω0 , we need to choose a band containing ω0 . For example, ω0 exists on the second band and the first band for the magnetic frozen mode and degenerate band edge respectively, but for the oblique frozen mode, no such information is known in advance. To identify the band where ω0 exist, we search whole bands for a wave vector k0 such as $ $ $' $ $ i+p dω $ $ (kj )$$ < ε, p ≥ 1, ε = small number D12 (k0 ) = min $ $j=i−p dk $
Figure 8: Array of crystals allowing an oblique frozen mode. The periodic structure consists of one anisotropic layer with εxz #= 0 and vacuum between them. The incident waves need to propagate into the layer at an oblique angle.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
1.38
1.32
1.36
1.31
1.29
V13
1.3
1.32
V13
1.34
15
1.3
1.28
1.28
1.27
1.26
1.26
1.24 -0.3
-0.2
-0.1
0
V4
0.1
0.2
0.3
1.25 -0.15
-0.1
-0.05
0
V4
0.05
0.1
0.15
Figure 9: Dispersion relation for the three-dimensional Maxwell problem with material parameters from paper [6] (left), and with the optimized solution for the 8th band(right). A stationary inflection point ω0 is observed in the dispersion relation only with the optimized solution
To increase the speed of the search process, we choose several bands in advance which seem to have k0 satisfying the above inequality. As for parameters for optimization, we use the oblique angle of incident wave(kx = ky ), the thickness of vacuum(LV ) and material constants of anisotropic layer(ε11 , ε33 ) with which material constant matrix is given as ε11 cos θ 2 + ε33 sin θ 2 0 (ε11 − ε33 ) cos θ sin θ , 0 ε22 0 εˆA = 2 2 (ε11 − ε33 ) cos θ sin θ 0 ε33 cos θ + ε11 sin θ
Among several stationary inflection points found by optimization, one of them with the lowest frequency ω0 is given as follows(also see Fig. 9): θ = π/4 ε11 = ε22 = 4.66812966, ε33 = 3.9656, LV = 0.5711535 $k = (−0.4800309, −0.4800309, −0.008982), ω0 = 1.295439 d2 ω dω = 4.82793e − 4, |ω=ω0 = 0.004393 dk ω=ω0 dk2 The dispersion relation resulting from the new parameters clearly show the correct properties for the three-dimensional crystal.
5.1
Time-Domain Modeling of the 3D Crystal
To observe the phenomena in three dimension, the new set of optimized parameters is solved in a time-domain wave propagation model with periodic boundary condition along x, y direction [20]. The periodic boundaries are x and y sides and absorbing boundary condition(ABC) is added at the end of z-direction to truncate unnecessary reflection of the outgoing waves(Fig. 10). Let only εxz be non-zero off-diagonal element for material
16
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
Figure 10: Wave propagation problem with periodic boundary condition along x and y direction with ABC at the end of z direction. The incident wave propagates into the layer with an oblique angle
constant matrix.
˜z ∂E ∂y ˜x ∂E ∂z ˜y ∂E ∂x
−
−
−
˜y ∂E ∂z ˜z ∂E ∂x ˜x ∂E ∂y
˜ Hx ∂ =− H ˜y , ∂t ˜z H
˜z ∂H ∂y ˜x ∂H ∂z ˜y ∂H ∂x
−
−
−
˜y ∂H ∂z ˜z ∂H ∂x ˜x ∂H ∂y
˜ Ex εxx 0 εxz = 0 εyy 0 ∂ E ˜y ∂t ˜z εxz 0 εzz E
˜ H ˜ as two dimensional Fourier Assuming periodicity along x, y direction, we can express E, series as follows. ˜ y, z, t) = { E(x,
Ny /2
'
Nx /2
Ny /2
Nx /2
'
2π
2π imy L y
E(mx , my , z, t)eimx Lx x e
y
my =Ny /2 mx =Nx /2
˜ H(x, y, z, t) = {
'
'
2π
2π imy L y
H(mx , my , z, t)eimx Lx x e
my =Ny /2 mx =Nx /2
}
y
}
where, Nx , Ny are number of points along each x and y direction and Lx , Ly are length of domain along x and y direction, respectively. Note that Lx and Ly depend on the angle of incidence of plane wave. By substituting these expansions into the above Maxwell’s equations after implementing ABC, we have ∂E sz 0 (imy )Ez − ∂zy ∂Ex 0 sz = iω ∂z − (imx )Ez 0 0 (imx )Ey − (imy )Ex ∂Hy εxx (imy )Hz − ∂z ∂Hx 0 = −iω − (im )H x z ∂z ε (imx )Hy − (imy )Hx xz
0 Hx 0 Hy Hz 1/sz 0 εxz sz 0 0 Ex εyy 0 0 sz 0 Ey 0 εzz Ez 0 0 1/sz
where ω stands for the frequency of the incident plane wave.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
1.1
17
45 40
1.05
35 1
TR angle
TR rate
30 0.95
0.9
25 20 15
0.85 10 0.8
0.75
5
0
10
20
30
40 50 Incidence angle
60
70
80
0
0
10
20
30
40 50 Incidence angle
60
70
80
Figure 11: Results of testing the code for 3D time domain modeling with periodicity along x, y direction. Transmittance rate vs. incidence angle(left), transmittance angle vs. incidence angle(right). Solid line indicates the exact value and circles indicate computed values.
Then, by substituting sz = 1 + σz /(iω), we have the following Maxwell equation in time domain( ∀ mx ≤ |Nx /2|, ∀ my ≤ |Ny /2|, Dε = εxx εzz − ε2xz ) ∂Hx ∂t ∂Hy ∂t ∂Bz ∂t ∂Hz ∂t ∂Ex Dε ∂t ∂Ey εyy ∂t ∂Dz Dε ∂t ∂Ez ∂t
∂Ey − (imy )Ez − σz Hx ∂z ∂Ex = − + (imx )Ez − σz Hy ∂z =
= −(imx )Ey + (imy )Ex = = = = =
∂Bz + σz Bz ∂t ∂Hy −εzz { − (imy )Hz } + εxz {(imy )Hx − (imx )Hy } − σz Ex ∂z ∂Hx − (imx )Hz − σz Ey ∂z ∂Hy εxz { − (imy )Hz } − εxx {(imy )Hx − (imx )Hy } ∂z ∂Dz + σz Dz dt
A Discontinuous Galerkin Method is used for differentiation along z-direction as in the paper [18, 19]. As for time marching, we use an implicit scheme for the linear derivatives along z and an explicit scheme for what remains [15, 16]. To validate the code, the transmittance rate and transmittance angle are computed when a plane incident wave propagating from vacuum to a dielectric material with various oblique angles. As shown in Fig. 11, it demonstrates the excellent performance of the scheme In Fig. 12, we show the energy density flux and the field amplitude with the optimized parameter values and observe several features of the frozen mode. In behalf of analysis,
18
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
2.2 2
2
0.4
1.8
1.8
1.6
1.6
1.4
1.4
1.2
1.2
fampz
EDFz
0.3
0.2
0.1
0
2.2
fampz
0.5
1 0.8
0.6
0.4
0.4
0.2
0.2
0 -0.1
0
50
z
0
-0.2 -50
100
1 0.8
0.6
0
50
z
100
150
200
-0.2 -50
0
50
z
100
150
200
Figure 12: Energy density flux(left) and the field amplitude with optimized parameters at T=100(center) and at T=200(right), Nx = Ny = 6, Nz = 12, Nx , Ny , Nz =number of grid points along x,y,z. Number of unit cell=100. Unit cell = A+V
0.05
0.5 0 0.4
0
0.3 -0.05
V4
V3
V2
-0.05 0.2
-0.1 0.1 -0.1 -0.15
-0.2 -100
0
0
100
V1
200
300
-100
0
100
V1
200
300
-0.1 -100
0
100
V1
200
300
Figure 13: Energy Flux Density Difference(EDFD) between the first cell and the last cell along x(left), y(center), z(right) direction. Solid line is with optimized parameters, dashed line is with original values. Nx = Ny = 6, Nz = 12, Number of unit cell=64, Unit cell = A+V
let’s define energy density flux difference(EDFD) as shown in [2] EDFD =
2
S
ef · nds =
2
V
∇ef dV ,
ef = energy density flux
note that this EDFD indicates the amount of energy density stacked inside the layer per unit time. Then, we observe that there are significant amounts of EDFD along the z direction between the first cell and the last cell. This leads to an increased field amplitude, compresses the energy density and the slow-down of waves. However, each of these phenomena are weaker than that of magnetic frozen mode as expected [2]. The maximum of the field amplitude is increasing very slowly and the average velocity of the energy transfer is approximately 0.4c. Fig. 13 compares EDFD for all directions for the crystals based on the original values and the optimized values. With the optimized values, we observe that the EDFD along x, y direction are almost zero, but EDFD along the z direction remains significantly large even after a long time. With original values, EDFD along the x, y direction remains nontrivial.
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
6
19
Conclusion
We have explored the use of PDE constrained optimization along with design and optimization of photonic crystals exhibiting frozen mode behavior. This phenomenon is related to the existence of Van Hove singularities in the dispersion relation. We have demonstrated the ability to recover known solutions as well as computed several new sets of parameters by optimizing the geometry and/or the materials of the crystals. For both of the magnetic frozen mode and degenerate band edge, each sets of parameters obtained by optimization with different number of variables show the generally similar phenomena, but some features are different. Above all, the transmittance rate varies significantly, which is a crucial element for energy density and field amplitude inside the layer. In particular one set for the degenerate band edge displays a transmittance rate around 0.3475, which is more than 50 times larger than those of the others. For the case with the frozen mode requiring an oblique illumination, we succeed in finding a stationary inflection point in three dimensions with parameters’ sets derived by nonlinear optimization. The existence of the frozen mode in three dimension is shown for the first time, but the phenomena themselves are not as distinctive as the previous two instances. One of the reason can be the intrinsic property of the oblique array which makes wave vector k0 too close to zero. The current computational design tool allows one to specify materials and/or geometries, then seek to design crystals with the sought-after properties, if they exist. The computations shown here illustrate, however, that the Van Hove singularities are found broadly in the dispersion relations, indicating that many configurations often exist. Future works may include finding stationary inflection points or degenerate band edge points showing strong frozen mode in 2 or 3 dimensional structures. We also believe that using a more complex geometry, but employing simpler materials, can be a candidate for exploring similar behavior.
Acknowledgements This work was supported by the U.S. Air Force Office of Scientific Research under the grant FA9550-04-1-0359. References [1] J.S. Hesthaven and T. Warburton, Discontinuous Galerkin Methods for the Time-Domain Maxwell’s Equations: An Introduction, ACES Newsletter, 19(1), 10-29, 2004. [2] Sehun Chun and Jan S. Hesthaven, Modeling of the frozen mode phenomenon and its Sensitivity using Discontinuous Galerkin methods, Communications in Computational Physics, to appear, 2006. [3] A. Figotin and I. Vitebsky, Nonreciprocal magnetic photonic crystals, Physical Review E, 63, 0666609, May, 2001.
20
S. Chun, J. S.Hesthaven / Commun. Comput. Phys., 1 (2006), pp. 1-21
[4] A. Figotin and I. Vitebsky, Electromanetic unidirectionality in magnetic photonic crystals, Physical Review B, 67, 165210, April, 2003. [5] A. Figotin and I. Vitebsky, Oblique frozen modes in periodic layered media, Physical Review E, 68, 036609, September, 2003. [6] J. Ballato and A. Ballato and A. Figotin and I. Vitebsky, Frozen light in periodic stacks of anisotropic layers, Physical Review E, 71, 1, 2005. [7] L. Van Hove, Van Hove singularities, Phys. Rev, 89, 1189, 1953. [8] M. Ibanescu and S.G. Johnson and D. Roundy and C. Luo and Y.Fink and J.D. Joannopoulos, Anomalous Disperion Relations by Symmetry Breaking in Axially Uniform Waveguides, Physical review letters, 92, 6, 2004. [9] K. Svanberg, The method of moving asymptotes: a new method for structural optimization, Int. J. Nuner. Meth. Engng, 24, 359-373, 1987. [10] C. Zillober, A globally convergent version of the method of moving asymptotes, Structural Optimization, 6, pp. 166-174, 1993. [11] W.H. Zhang and C. Fleury and P. Duysinx and V.H. Nguyen and I. Laschet, A generalized method of moving asymptotes(GMMA) including equality constraints, Structural Optimization, 12, pp. =143-146, 1996. [12] K. Svandberg, A new globally convergent version of MMA, Technical Report TRITA-MAT1999-OS2, Dept. of Mathematics, KTH, Stockholm, 1999. [13] K. Svanberg, Two primal-dual interior-point methods for the MMA subproblems. Technical Report TRITA-MAT-1998-OS12, Dept. of Mathematics, KTH, Stockholm, 1998. [14] Johnson, Steven G. and Joannopoulos, J. D., Block-iterative frequency-domain methods for Maxwell’s equations in a planewave basis, Opt. Express, 2001, 8, 3, pp. 173-190, url =http://www.opticsexpress.org/abstract.cfm?URI=OPEX-8-3-173. [15] Uri M. Ascher and Steven J. Ruuth and Raymond J. Spiteri, Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Applied Numerical Mathematics, 1997, 25, pp. 151-167. [16] Christopher A. Kennedy and Mark H. Carpenter, Additive Runge-Kutta schemes for convection-diffusion-reaction equations, Applied Numerical Mathematics, 2003, 44, pp. 139181. [17] Alex Figotin and Ilya Vitebskiy, Gigantic transmission and band-edge resonance in periodic stacks of anisotropic layers, Physical review E, 2005, 72, 036619. [18] J.S. Hesthaven and T. Warburton, Discontinuous Galerkin Methods for the Time-Domain Maxwell’s Equations: An Introdunction, ACES Newsletter, 19(1), pp. 10-29, 2004. [19] J.S. Hesthaven and T. Warburton, Nodal High-Order Methods on Unstructured Grids, Journal of Computational Physics, Journal of Computational Physics, 182, pp. 186-221, 2002. [20] Allen Taflove and Susan C. Hagness, Computational Electrodynamics, the Finite-Difference Time-Domain Method, 3rd edition, Artech House, Boston, 2005. [21] Uri M. Ascher and Linda R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, siam, Philadelpha, 1998. [22] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis 2nd edition, Springer-Verlag, New York, 1979.
