Coupled Mode Equation Modeling for Out-of-Plane Gap Solitons in 2D Photonic Crystals

# Coupled Mode Equation Modeling for Out-of-Plane Gap Solitons in 2D Photonic Crystals

Tomáš Dohnal111Institute for Applied and Numerical Mathematics, Karlsruhe Institute of Technology, D–76128 Karlsruhe, Germany  and Willy Dörfler111Institute for Applied and Numerical Mathematics, Karlsruhe Institute of Technology, D–76128 Karlsruhe, Germany
July 9, 2019
###### Abstract

Out-of-plane gap solitons in 2D photonic crystals are optical beams localized in the plane of periodicity of the medium and delocalized in the orthogonal direction, in which they propagate with a nonzero velocity. We study such gap solitons as described by the Kerr nonlinear Maxwell system. Using a model of the nonlinear polarization, which does not generate higher harmonics, we obtain a closed curl-curl problem for the fundamental harmonic of the gap soliton. For gap solitons with frequencies inside spectral gaps and in an asymptotic vicinity of a gap edge we use a slowly varying envelope approximation based on the linear Bloch waves at the edge and slowly varying envelopes. We carry out a systematic derivation of the coupled mode equations (CMEs) which govern the envelopes. This derivation needs to be carried out in Bloch variables. The CMEs are a system of coupled nonlinear stationary Schrödinger equations with an additional cross derivative term. Examples of gap soliton approximations are numerically computed for a photonic crystal with a hexagonal periodicity cell and an annulus material structure in the cell.

g

ap soliton, photonic crystal, Maxwell’s equations, Kerr nonlinearity, out of plane propagation, coupled mode equations, slowly varying envelope approximation, Bloch transformation

{AMS}

41A60, 35Q61, 35C20, 78M35

## 1 Introduction

Maxwell’s equations for electromagnetic waves in Kerr nonlinear dielectric materials read

 ∂tD =∇×H, (1.1a) μ0∂tH =−∇×E, (1.1b) ∇⋅D =0, (1.1c) ∇⋅H =0 (1.1d)

for the electric field , magnetic field , the electric displacement field with the constitutive relations

 D=ε0(n2E+PNL),PNL,i=3∑j,l,m=1χ(3)ijlmEjElEmfor i=1,2,3. (1.2)

are the electric permittivity and magnetic permeability of vacuum, respectively, is the refractive index of the medium, and is the cubic electric susceptibility of the medium.

We consider a 2D photonic crystal, i.e. we assume that the material coefficients change periodically on a plane and are independent of the orthogonal component on that plane. Let be linearly independent lattice vectors defining the Bravais lattice of the crystal. Then the required periodicity reads

 n(x)=n(x+R)∈R,χ(3)(x)=χ(3)(x+R)∈R3×3×3×3for all x∈R3 and % all R∈Λ. (1.3)

Without loss of generality we assume that the crystal is homogeneous in the -direction, i.e. and for all . We denote by the Wigner–Seitz cell corresponding to the Bravais lattice. We use to denote the pair of vectors satisfying for , and let the reciprocal lattice be . denotes the first Brillouin zone, i.e. the Wigner–Seitz cell of the reciprocal lattice.

Note that from the relations in (1.2) it is clear that we are neglecting losses, material dispersion as well as higher order nonlinearities and assuming that the third order nonlinear response of the medium is instantaneous.

We will consider monochromatic waves propagating in the -direction, i.e. waves propagating out of the plane of periodicity of the 2D crystal, and use the ansatz

 (E,H,D)(x,t)=ei(κx3−ωt)(E,H,D)(x1,x2;ω)+c.c., (1.4)

where and c.c. denotes the complex conjugate of the first term on the right. The ansatz (1.4) contains no higher harmonics, which is valid if the above form of is replaced by a time averaged one, see below. Alternatively, a physical justification of neglecting higher harmonics is based on the lack of phase matching and absorption.

Note that for the field (1.4) the divergence free conditions (1.1c) and (1.1d) are automatically satisfied provided since the spatially dependent parts

 (^E,^H,^D)(x;ω):=eiκx3(E,H,D)(x1,x2;ω)

satisfy

 ^D=iω∇×^Handμ0^H=−iω∇×^E, (1.5)

and thus . Since our analysis below is for gap solitons with close to a band edge, the condition is for us restrictive only when is in a gap and lies near a band edge. Note also that even if higher harmonics are accounted for, the divergence free conditions are still satisfied for as (1.5) then holds for each generated harmonic. Clearly, only odd, i.e., -th, , harmonics are generated.

We will assume a centrosymmetric and isotropic -tensor, which leads to the simplification

 PNL=χ(3)ci(E⋅E)E,

where for , see [21, Sec. 2d]. Inserting the ansatz (1.4) in the nonlinearity clearly generates the harmonics . These are, however, typically neglected based on the physical arguments that the fundamental harmonics and the higher harmonics are not phase matched and that at the higher values of frequency (i.e. at ) material absorption is usually large preventing the generation of significant fields at these frequencies, see e.g. [9]. Considering only the fundamental harmonics, the nonlinear polarization for the ansatz (1.4) becomes

 PNL=χ(3)ci(2|E|2E+E⋅E¯E)ei(κx3−ωt)+c.c., (1.6)

i.e.

 PNL=χ(3)ci((3|E1|2+2|E2|2+2|E3|2)E1+(E22+E23)¯E1(2|E1|2+3|E2|2+2|E3|2)E2+(E21+E23)¯E2(2|E1|2+2|E2|2+3|E3|2)E3+(E21+E22)¯E3)ei(κx3−ωt)+c.c.. (1.7)

Another widely used model for the nonlinear polarization is

 PNL=χ(3)ci[E⋅E]avE,

where denotes the time average of over the period of , i.e. over for , cf. [26, 27]. The averaging generates no higher harmonics so that in this model (1.6) is exact. Note that the Kerr nonlinear problem including all higher harmonics has been recently considered for a 1D periodic structure in [25].

In the following we rescale the frequency by defining

 ˜ω:=ωc

but drop the tilde again for better readability. For convenience we will denote the square of the refractive index by

 η(x):=n2(x)for all x∈R3.

With the ansatz (1.4) equations (1.1a) and (1.1b) become

 −icωD =∇×H+i(00κ)×H, (1.8a) icωμ0H =∇×E+i(00κ)×E. (1.8b)

Since all our functions are independent of , we let from now on . Using the fact that depends only on and , a second order formulation of (1.8) reads

 (L−ω2η)E=ω2PNL, (1.9)

where

 LE:=∇×∇×E+iκ(∂x1E3∂x2E3∂x1E1+∂x2E2)+κ2(E1E20), (1.10)

and

Having determined , the magnetic field can be recovered by

Based on the analogy with the periodic nonlinear Schrödinger equation [23], equation (1.9) is expected to have localized -solutions for any in a spectral gap of the linear problem . Such solutions are called gap solitons. The aim of this paper is to provide an approximation of gap solitons of (1.9) for in an -vicinity () of a gap edge using a slowly varying envelope approximation. As we show, envelopes of such gap solitons satisfy a system of nonlinear constant coefficient equations, so called coupled mode equations (CMEs) posed in the slow variables . The CMEs can be numerically solved with less effort than the nonlinear Maxwell system (1.9) in the variable . An asymptotic approximation of a gap soliton of (1.9) near a gap edge is then the sum of linear Bloch waves at the edge, modulated by the corresponding envelopes.

Asymptotic approximations via CMEs have been analyzed for gap solitons of the stationary periodic nonlinear Schrödinger equation in 1D [24] as well as in 2D [12, 13, 14]. In these works the approximation via CMEs was also rigorously justified using Lyapunov–Schmidt reductions. Gap solitons of the nonlinear Maxwell’s equations have been approximated by CMEs in the case of 1D photonic crystals with a small (infinitesimal) contrast in the periodicity [16, 24, 25], where [16] considers gap solitons modulated also in time. To our knowledge the problem of a systematic CME approximation of gap solitons of nonlinear Maxwell’s equations describing 2D or 3D photonic crystals does not appear in the literature. Although CMEs have been formally derived for pulses in Maxwell’s equations with a 2D periodic medium of small contrast [2, 1, 11], these pulses cannot be true gap solitons because in 2D and 3D a large enough contrast is necessary for the opening of spectral gaps. In this paper we consider a 2D photonic crystal with a finite contrast in the periodicity. For our examples we use a photonic crystal which has several spectral gaps [4].

Besides the above cited works on coupled mode modeling of gap solitons there are a number of papers on the slowly varying envelope approximation of nonlinear pulses in periodic structures with the pulse frequency lying within the spectral bands. The envelope in this case can be typically modeled by the time dependent nonlinear Schrödinger equation and the approximation holds on large but finite time intervals [9, 6, 10].

The rest of the paper is organized as follows. In Section 2 we study the linear band structure of (1.9) (with ) and obtain thus the linear spectrum of the problem. We also discuss possible symmetries in the band structure and among the corresponding Bloch waves. An example of a photonic crystal from [4] is then provided, for which the band structure is numerically computed and three band gaps are observed on the positive half axis . In Section 3 we present a slowly varying envelope approximation of gap solitons of (1.9) for in the vicinity of a spectral edge and carry out a systematic formal derivation of CMEs describing the envelopes. Next, examples of CMEs are presented for the concrete photonic crystal given in Section 2 as well as for other theoretical situations. Here the symmetries in the band structure and among the Bloch waves play an important role in determining properties of the CME coefficients. In Section 4 we plot the approximation of two gap solitons in the chosen photonic crystal. The approximation requires computing the Bloch waves at the edge and solving the corresponding CMEs.

## 2 Linear Band Structure

### 2.1 The periodic eigenvalue problem

We study first the linear problem

 Lu=ω2ηuon R2 (2.1)

and define the band structure as well as the linear Bloch waves.

By the Bloch–Floquet theory, see [19] or [15, Ch. 3], solution modes of (2.1) are given by the Bloch waves for that satisfy

 Lun(k;.)=ωn(k)2ηun(k;.),un(k;.+R)=un(k;.)eik⋅Rfor all R∈Λ, (2.2)

where sweeps the first Brillouin zone .

It is well-known that is self-adjoint and has a compact inverse and that there thus exists a sequence of eigenvalues with and each eigenspace is of finite dimension. These eigenvalues are nonnegative and we use the natural ordering for . The mapping is called the -th band of the spectral problem (2.2). Of course, (2.2) allows also non-positive bands . These are typically labeled via and will play no role in our analysis. We therefore restrict ourselves to for . The Bloch waves in (2.2) can be written in the form

 un(k;x)=pn(k;x)eik⋅x,

where the are -periodic in , i.e. for all , . These satisfy the eigenvalue problem

 (˜L(k)−ω2n(k)η(x))pn(k;x)=0for all x∈U,pn(k;x+R)=pn(k;x)for all x∈∂U and all R∈Λ, (2.3)

with

 ˜L(k)pn(k;x)=(∇+ik′)×(∇+ik′)×pn(k;x),

where , . Since is -independent, can be written as

 ˜L(k)=(κ2−(∂x2+ik2)2(∂x1+ik1)(∂x2+ik2)iκ(∂x1+ik1)(∂x1+ik1)(∂x2+ik2)κ2−(∂x1+ik1)2iκ(∂x2+ik2)iκ(∂x1+ik1)iκ(∂x2+ik2)−(∂x1+ik1)2−(∂x2+ik2)2).

In the variable the Bloch waves and the eigenvalues are easily proved to fulfill

 ωn(k)=ωn(k+K),pn(k+K;x)=pn(k;x)e−iK⋅xfor all x∈U,K∈Λ∗. (2.4)

Due to the self-adjoint nature of we can normalize the Bloch functions via

 ⟨pn(k;.),ηpm(k;.)⟩=δn,m, (2.5)

where for .

For purposes of the later asymptotic analysis of gap solitons we also present calculations of first and second order derivatives of the bands at extremal points. Suppose the band has an extremum at and denote . By direct differentiation of (2.3) we see that the “generalized Bloch functions” , for , are solutions of the system

 (˜L(k∗)−ω2∗η)∂kjpn∗(k∗;.)=−∂kj˜L(k∗)pn∗(k∗;.). (2.6)

Applying the differentiation , for , to (2.3) and evaluation at , yields

 (˜L(k∗)−ω2∗η(x))∂2ki,kjpn∗(k∗;x) =2ω∗η(x)∂2ki,kjωn∗(k∗)pn∗(k∗;x)−∂2ki,kj˜L(k∗)pn∗(k∗;x) −∂ki˜L(k∗)∂kjpn∗(k∗;x)−∂kj˜L(k∗)∂kipn∗(k∗;x).

Necessarily, due to the Fredholm alternative, the right hand side is -orthogonal to , which lies in the kernel of with periodic boundary conditions on . This yields the formula

 (∂2kωn∗(k∗))i,j=∂2ki,kjωn∗(k∗) =12ω∗⟨∂2ki,kj˜L(k∗)pn∗(k∗;.)+∂ki˜L(k∗)∂kjpn∗(k∗;.)+∂kj˜L(k∗)∂kipn∗(k∗;.),pn∗(k∗,⋅)⟩. (2.7)

A straightforward differentiation of yields

 ∂k1˜L(k∗) =(0i(∂x2+ik∗,2)−κi(∂x2+ik∗,2)−2i(∂x1+ik∗,1)0−κ0−2i(∂x1+ik∗,1)), ∂k2˜L(k∗) =(−2i(∂x2+ik∗,2)i(∂x1+ik∗,1)0i(∂x1+ik∗,1)0−κ0−κ−2i(∂x2+ik∗,2)),
 ∂2k1˜L≡(000020002),∂2k2˜L≡(200000002),and∂2k1,k2˜L≡(−0−1−0−1−0−0−0−0−0), (2.8)

where , for , is the -th component of . With these the explicit forms of (2.1) read

 ∂2k1ωn∗(k∗)=1ω∗⟨(i(∂x2+ik∗,2)∂k1pn∗,2(k∗;.)−κ∂k1pn∗,3(k∗;.)i(∂x2+ik∗,2)∂k1pn∗,1(k∗;.)−2i(∂x1+ik∗,1)∂k1pn∗,2(k∗;.)+pn∗,2(k∗;.)−2i(∂x1+ik∗,1)∂k1pn∗,3(k∗;.)−κ∂k1pn∗,1(k∗;.)+pn∗,3(k∗;.)),pn∗(k∗;.)⟩, (2.9)
 ∂2k2ωn∗(k∗)=1ω∗⟨(−2i(∂x2+ik∗,2)∂k2pn∗,1(k∗;.)+i(∂x1+ik∗,1)∂k2pn∗,2(k∗;.)+pn∗,1(k∗;.)i(∂x1+ik∗,1)∂k2pn∗,1(k∗;.)−κ∂k2pn∗,3(k∗;.)−κ∂k2pn∗,2(k∗;.)−2i(∂x2+ik∗,2)∂k2pn∗,3(k∗;.)+pn∗,3(k∗;.)),pn∗(k∗;.)⟩, (2.10)

and

 ∂2k1,k2ωn∗(k∗)=12ω∗⟨(−2i(∂x2+ik∗,2)∂k1pn∗,1(k∗;.)+i(∂x1+ik∗,1)∂k1pn∗,2(k∗;.)+i(∂x2+ik∗,2)∂k2pn∗,2(k∗;.)i(∂x1+ik∗,1)∂k1pn∗,1(k∗;.)+i(∂x2+ik∗,2)∂k2pn∗,1(k∗;.)−2i(∂x1+ik∗,1)∂k2pn∗,2(k∗;.)−2i(∂x2+ik∗,2)∂k1pn∗,3(k∗;.)−2i(∂x1+ik∗,1)∂k2pn∗,3(k∗;.))+(−κ∂k2pn∗,3(k∗;.)−pn∗,2(k∗;.)−κ∂k1pn∗,3(k∗;.)−pn∗,1(k∗;.)−κ(∂k1pn∗,2(k∗;.)+∂k2pn∗,1(k∗;.))),pn∗(k∗;.)⟩. (2.11)

### 2.2 Symmetries of the Band Structure and the Bloch waves

Symmetries in the refractive index function yield symmetries in the band structure and among Bloch waves. We restrict our attention to the cases of discrete rotational and axial reflection symmetry, which are relevant for the example we present below. The results of this section will be important when determining properties of the coefficients of coupled mode equations in Section 3.4.

#### 2.2.1 Rotational symmetry

Assume that the photonic crystal satisfies the rotational symmetry

 η(x)=η(rα(x))for all x∈R2 (2.12)

for some with the rotation defined by

 rα(x)=(cos(α)x1−sin(α)x2sin(α)x1+cos(α)x2)

Below we use the notation if is a two dimensional vector and if is a three dimensional vector .

The symmetry (2.12) implies a symmetry of the Rayleigh quotient corresponding to the eigenvalue problem (2.3) and thus a symmetry of the band structure. In detail, for we have

 ω2n(k)=minV⊂Hcurlper% (U)dimV=nmaxw∈V,w≠0∫U|(∇+ik′)×w(x)|2dx∫Uη(x)|w(x)|2dx,

and the corresponding extremal point is . Due to the relation

 ((∇+irα(k′))×f)(rα(x))=rα[(∇+ik′)×r−α(f(rα(x)))]for all smooth f:R2→R3

we get

 ∫U|(∇+irα(k′))×w(x)|2dx=∫U|(∇+ik′)×r−α(w(rα(x)))|2dx,

and symmetry (2.12) yields

 ∫Uη(x)|w(x)|2dx=∫Uη(x)|r−α(w(rα(x)))|2dx.

As a result we obtain that

 ωn(k)=ωn(rα(k))for all n∈N and all k∈B. (2.13)

If has geometric multiplicity one as an eigenvalue of (2.3), we have also a symmetry of the corresponding Bloch functions, namely

 (2.14)

Note that a renormalization of , in order to obtain in (2.14), is in general impossible when , where reads “ congruent to ” and means for some . This is because in this case and are related by (2.4) and a renormalization of the left hand side of (2.14) would affect the right hand side in the same way. When is not congruent to , e.g. when , then one can set in (2.14).

¿From the symmetry (2.13) we can deduce a symmetry of the second derivatives of . Using the identity , we get by further differentiation

 (2.15)

for all and .

#### 2.2.2 Reflection symmetry

If the photonic crystal satisfies the reflection symmetry

 η(x)=η(S1(x))for all x∈R2, where S1(x)=(−x1,x2)T, (2.16)

then similarly to Section 2.2.1 we have

 ωn(k)=ωn(−k1,k2)for all k∈B and n∈N. (2.17)

Again, if has geometric multiplicity one as an eigenvalue of (2.3), then

 pn(S1(k);x)=eiaS1(pn(k;S1(x)))% for all n∈N and some a=a(n)∈R, (2.18)

where for . Just as above, unless , we can set in (2.18). The symmetry (2.17) implies

 ∂2k1ωn(k)=(∂2k1ωn)(−k1,k2),∂2k2ωn(k)=(∂2k2ωn)(−k1,k2),∂2k1,k2ωn(k)=−(∂2k1,k2ωn)(−k1,k2) (2.19)

for all and .

An analogous discussion, of course, applies for the reflection symmetry for all , where . One the obtains

 ∂2k1ωn(k)=(∂2k1ωn)(k1,−k2),∂2k2ωn(k)=(∂2k2ωn)(k1,−k2),∂2k1,k2ωn(k)=−(∂2k1,k2ωn)(k1,−k2) (2.20)

for all and and if has geometric multiplicity one as an eigenvalue of (2.3), then

 pn(S2(k);x)=eiaS2(pn(k;S2(x)))% for all n∈N and some a=a(n)∈R. (2.21)

#### 2.2.3 Combination of rotational and reflection symmetries

If both the reflection symmetry (2.16) and the rotational symmetry (2.12) for some , , hold, then for along the rays with angles and the mixed derivative can be expressed in terms of and . This is because for along these rays we have or , so that both (2.15) and (2.19) or (2.20) apply. In detail, suppose

 (−k1,k2)=rα(k), i.e.{} k=|k|ei(π/2−α/2)ork=|k|e−i(π/2+α/2)=−|k|ei(π/2−α/2).

Then it follows that

 ∂2k1ωn(k)=(∂2k1ωn)(−k1,k2)=cos2(α)∂2k1ωn(k)−sin(2α)∂2k1,k2ωn(k)+sin2(α)∂2k2ωn(k),

where the first equality holds due to (2.19) and the second due to (2.15). As a result, for , , and we get

 ∂2k1,k2ωn(k)=12tan(α)(∂2k2ωn(k)−∂2k1ωn(k)). (2.22)

Identity (2.22) applies also in the case when the reflection symmetry and the rotational symmetry (2.12) are both present for some , . Then (2.22) holds for that satisfy

 (k1,−k2)=rα(k), i.e.{} k=±|k|e−iα/2.

### 2.3 Example: Hexagonal Lattice with a Circular Material Structure

As an example we consider the hexagonal lattice in the -plane generated by the vectors

 a(1)=a0(cos(π/3)sin(π/3))anda(2)=a0(10)witha0>0.

In the Wigner–Seitz cell the material structure is given by the annulus centered at the lattice point in the origin and having outer and inner radii and respectively. The material properties are given by for and otherwise. This is the same as the crystal used in [4], where the corresponding band structure was also computed. One choice of vectors generating the reciprocal lattice is

where