Grating and Bragg Grating - An example of optical field in periodic structure

(Advanced - 6358 or above)

Utility Modules

CPPlot[z_List,range_,sizeparm_,hueparm_,label_]:=Module[{R,n,t,r,h, huelist, m},
 n=Length[z] ;
If[TrueQ[range== Null],R=Max[Abs[z]] ,R=range]  ;
If[TrueQ[sizeparm== Null],{t,h,r}={0.005,0.05,0.025} ,{t,h,r}=sizeparm]  ;
If[TrueQ[hueparm== Null] ,
huelist=Table[(i-1) 0.95/n,{i,1, n}] 
 ,huelist=hueparm]  ;  m=Length[huelist]+1;
If[TrueQ[label== Null],
Graphics[{Arrow[{{-1.2R,0},{1.2 R,0}}],Arrow[{{0,-1.2 R},{0,1.2 R}}]
,Inset["Re",{1.1 R,-0.1}]
,Inset["Im",{-0.1,1.1 R} ]
,Table[{Hue[ huelist[[ Mod[i,m]  ]],1,1]
,Opacity[0.75],Disk[{Re[z[[i]] ],Im[z[[i]] ]},r]
,Thickness[t],Arrowheads[h],Arrow[{{0,0},{Re[z[[i]] ],Im[z[[i]] ]}}]},{i,1, n}]},
PlotRange-> {1.2*{-R,R},1.2*{-R,R}},
Frame-> True, GridLines-> Automatic]
,
Graphics[{Arrow[{{-1.2R,0},{1.2 R,0}}],Arrow[{{0,-1.2 R},{0,1.2 R}}]
,Inset["Re",{1.1 R,-0.1}]
,Inset["Im",{-0.1,1.1 R} ]
,Table[{Hue[ huelist[[ Mod[i,m]  ]],1,1]
,Inset[label[[i]],1.15 {Re[z[[i]] ],Im[z[[i]] ]} ]
,Opacity[0.75],Disk[{Re[z[[i]] ],Im[z[[i]] ]},r]
,Thickness[t],Arrowheads[h],Arrow[{{0,0},{Re[z[[i]] ],Im[z[[i]] ]}}]},{i,1, n}]},
PlotRange-> {1.2*{-R,R},1.2*{-R,R}},
Frame-> True, GridLines-> Automatic]
]
]

GratingDiffG[λ_,θi_,d_,ϵf2_,mEvan_] - current - useable

GratingDiffG[λ_,θi_,d_,ϵf2_,mEvan_]:=Module[{ϵ2,nindex, nI,γ,k0,K1,K2, mminR,mmaxR
(* This version generate k vectors for plotting
with horizontal x, and y vertical is positive, so that downward k vectors are negative *) 
,mER,mET, U,U0,UEvan,Uk,UkEvan,β1,κ1
, mminT,mmaxT,V,V0,VEvan,Vk,VkEvan,k2,k2Evan},
(*  ef2 is a function input for complex dielectric with absorption  *)
If[NumberQ[ϵf2]== True ,ϵ2=ϵf2,ϵ2=ϵf2[λ] ] ;
{nindex, nI}={Re[(Sqrt[ϵ2])],Im[(Sqrt[ϵ2])]} ;
If[TrueQ[mEvan==Null]== True  ,{mER,mET}={2,2}, 
If[ListQ[mEvan] && Length[mEvan]≥  2,{mER,mET}=Take[mEvan,2],If[ListQ[mEvan] && Length[mEvan]==  1,{mER,mET}={mEvan[[1]],mEvan[[1]]} 
,{mER,mET}={mEvan,mEvan}] ]] ;
If[mER≤ 0, mER=1, ]; If[mET≤ 0, mET=1, ] ;
k0=2 π/λ ;K1=k0; 
(*   Reflected diffraction  *)
mminR=Min[{Ceiling[((d )/λ) (Sin[θi]-1)  ],Floor[ ((d )/λ)  (Sin[θi]+1)  ] }];
mmaxR=Max[{Ceiling[((d )/λ) (Sin[θi]-1)  ],Floor[ ((d )/λ)  (Sin[θi]+1)  ] }];
U=Table[Sin[θi]-m (λ/d),{m,mminR-mER,mmaxR+mER,1}] ;
U0=Select[U,Abs[#]≤ 1 & ] ; UEvan=Select[U,Abs[#]> 1 & ] ;
Uk=Table[{{K1 U0[[i]],K1(Sqrt[1 -((U0[[i]])^2) ])},mminR+i-1} ,{i,1,Length[U0]}] ;
UEvan=Sort[UEvan,Abs[#1]<Abs[#2]  & ] ;
(*   Print["UEvan=    ",UEvan];     *)
β1=K1*UEvan ;κ1=-(Sqrt[((β1)^2)-(K1^2)]) ;
UkEvan=Table[{β1[[i]],κ1[[i]]} ,{i,1,mER}] ;
(*   Transmitted diffraction  *)
mminT=Min[{Ceiling[((d )/λ) (Sin[θi]-nindex)  ],Floor[ ((d )/λ)  (Sin[θi]+nindex)  ] }];
mmaxT=Max[{Ceiling[((d )/λ) (Sin[θi]-nindex)  ],Floor[ ((d )/λ)  (Sin[θi]+nindex)  ] }];
V=Table[Sin[θi]-m (( λ)/d),{m,mminT-mET, mmaxT+mET,1}] ;
V0=Select[V,((Abs[#])/nindex)≤ 1 & ] ; VEvan=Select[V,((Abs[#])/nindex)> 1 & ] ;
k2=-k0 (Sqrt[ ϵ2-(V0^2)]); k2Evan=k0(Sqrt[ ϵ2-(VEvan^2)]) ;
Vk=Table[{{k0 V0[[i]],k2[[i]]},mminT+i-1} ,{i,1,Length[V0]}] ;
VEvan=Sort[VEvan,Abs[#1]<Abs[#2]  & ] ;
VkEvan=Table[{k0 VEvan[[i]],k2Evan[[i]]} ,{i,1,mET}] ;
{{Uk,UkEvan ,ArcSin[U0]},{Vk,VkEvan, ArcSin[(V0/(Sqrt[(V0^2)+((Re[k2/k0])^2)]))]}}
]

FieldArray0[fieldparm_,xyarr_]:=Module[{A,β,k,x,y,U},
{A,β,k}=fieldparm ; {x,y}=xyarr ;
U=Table[{x[[i]],y[[j]],Re[A Exp[I (β x[[i]] +k y[[j]] )]]},{j,1,Length[y]},{i,1,Length[x]}] ;
{U, Transpose[Flatten[U,1]][[3]]}] ;

DiffractDiagram:=Manipulate[

DiffractDiagram:=Manipulate[  dw=d/2.5 ;Ng=Ceiling[3/dw];
mEvan={mE,mE}; 
(*   Graphics for grating and diffraction disk  *)

Grating=Table[Rectangle[{(i-0.25) dw,-0.2*h},{(i+0.25) dw,0.2 h}],{i,-Ng,Ng}];
PrependTo[Grating,Rectangle[nindex*{-1.2,-1.2},{nindex*1.2,-0.2*h} ] ];
diffdisk={Hue[0.16,1,1,0.4],Disk[{0,0},1,{0,π}],Hue[0.3,1,1,0.3],Disk[{0,0},nindex,{-π,0}]  } ;
pi={-Sin[θi],Cos[θi]}; k0=2 π/λ ; K1=k0 ; ϵ2=(((nindex+I nI))^2) ;
ϵ2R=Re[ϵ2] ; ϵ2I=Im[ϵ2] ;
(*   Output from GratingDiffG module   *)
{{Uk,UkEvan ,θdR},{Vk,VkEvan,θdT}} = GratingDiffG[λ,θi,d, ϵ2,mEvan] ;
nR=Length[Uk] ;
mmin=Vk[[1,2]] ; mmax=Last[Vk][[2]] ;
(*   β vertical grid graphics    *)
BETApoint=Table[{Red,Disk[{Sin[θi]-m λ/d,0},0.035],White,Disk[{Sin[θi]-m λ/d,0},0.02]},{m,mmin-2,mmax+2}]; 
BETAgrid=Table[Line[{{Sin[θi]-m λ/d,-Cos[θdT[[m-mmin+1]] ](Sqrt[((Vk[[m-mmin+1,1,1]])^2)+((Re[Vk[[m-mmin+1,1,2]]])^2)])/k0},{Sin[θi]-m λ/d
,If[(((Sin[θi]-m λ/d))^2)≤ 1,(Sqrt[1-(((Sin[θi]-m λ/d))^2)]), 0]}}]
,{m,mmin,mmax}]; 

 {β1,κ1}=Transpose[UkEvan]/K1  ;
PD=Table[{Hue[(i-0.5)/(nR+1),1,0.9,0.8],Arrow[{{0,0},{Sin[ θdR[[i]] ],Cos[θdR[[i]] ]}}] },{i,1,nR}]  ;
Mlabel=Table[Inset[Uk[[i,2]],1.1*{Sin[ θdR[[i]] ],Cos[θdR[[i]] ]}
,BaseStyle->Directive[14,FontFamily->"Helvetica"]],{i,1,nR}]  ; 

PDEvan=Table[{Hue[(i-1)/2,1,0.75,0.3],Rectangle[{β1[[i]],0 }, {0,-κ1[[i]]  }] 
,Hue[(i-1)/2,1,0.9,0.8],Thickness[0.005],Arrowheads[0.035]
,Arrow[{{0,0},{β1[[i]],0 }}],Arrow[{{0,0},{0,-κ1[[i]] }  }]
 },{i,1,mEvan[[1]] }]  ; 

(*   Reflected diffraction - evanescent field  *)

(*   Transmitted diffraction  *)

nT=Length[Vk] ; 

QD=Table[{Hue[(i-0.5)/(nT+1),1,0.9,0.8],Arrow[{{0,0},((Sqrt[((Vk[[i,1,1]])^2)+((Re[Vk[[i,1,2]]])^2)])/k0){Sin[θdT[[i]] ],-Cos[θdT[[i]] ]}}]
},{i,1,nT}]  ;

Nlabel=Table[Inset[Vk[[i,2]] ,1.1*((Sqrt[((Vk[[i,1,1]])^2)+((Re[Vk[[i,1,2]]])^2)])/k0){Sin[θdT[[i]]],-Cos[θdT[[i]] ] }
,BaseStyle->Directive[14,FontFamily->"Helvetica"]],{i,1,nT}]  ; 
QDEvan=Table[{Hue[(i-1)/2,1,0.75,0.3]
,Polygon[{{0,0},{VkEvan[[i,1]],-Re[VkEvan[[i,2]]  ]}/k0 , {VkEvan[[i,1]],-Re[VkEvan[[i,2]]  ]  -Im[VkEvan[[i,2]]   ] }/k0, {0,-Im[VkEvan[[i,2]]   ]  }/k0} ] 
,Hue[(i-1)/(1+1),1,0.9,0.8],Thickness[0.005],Arrowheads[0.035]
,Arrow[{{0,0},{VkEvan[[i,1]],-Re[VkEvan[[i,2]]  ] }/k0}]
,Arrow[{{0,0},{0,-Im[VkEvan[[i,2]] ] }/k0}]
},{i,1,mEvan[[2]]}]  ;
{Graphics[{Hue[0.6,0.5,1,1],Grating, diffdisk
,Hue[0,0,0,1],Dashed,Thickness[0.0001], Line[{{0,0},{0,1}}]
,{Red,BETAgrid},BETApoint
,Dashing[1],Thickness[0.0075],Arrowheads[.075],Arrow[{pi,{0,0}}]
,PD,Mlabel ,QD,Nlabel}, PlotRange->nindex* {{-1.1,1.1},{-1.1,1.1}}
,ImageSize->{500,500}]
 ,Graphics[{Hue[0.6,0.5,1,1],Grating, diffdisk
,Hue[0,0,0,1],Dashed,Thickness[0.0001], Line[{{0,0},{0,1}}]
,Dashing[1],Thickness[0.0075],Arrowheads[.075],Arrow[{pi,{0,0}}]
,BETApoint
,PDEvan,QDEvan}
,PlotRange->Max[Abs[β1[[1]] ], Abs[VkEvan[[1,1]] ] ]/k0 * {{-1.05,1.05},{-1.05,1.05}}
,ImageSize->{500,500}]  }
,{{λ,0.5},0.4,1.5},{{d,1,"Grating period"},0.25,4}
,{{h,0.1,"Tooth depth"},0.1,2},{{nindex,1.5,"refractive index"},1.,2},{{nI,0,"Im index"},0,0.5},{{θi,π/4.,"incident θ"},0,π/2}
,{{mE,1,"Evanescence"},1,5,1}]

BraggStackGr[Λ_,η_,n_]

BraggStackGr[Λ_,η_,n_,hueparm_]:=Module[{bg,hue},
If[TrueQ[hueparm==Null]== True,hue={0.35,0.13} ,hue=hueparm];
Grating=Table[{Hue[hue[[1]],0.5,1,0.8],Rectangle[{(i+η) Λ,-0.5},{(i+1) Λ,0.5 }]},{i,0,n}];
bg={Hue[hue[[2]],1,1,0.5],Rectangle[{0,-0.5},{(n+1) Λ,0.5 }]} ;
PrependTo[Grating,bg]  ; 
Graphics[Grating]
] ;

1. Introduction

1.1 General discussion

Grating is one of the most important class of structurely simple optical device. 
(see link to Photonics Device in Optical Communications)
It is a wavelength-sensitive device, with chief applications in wavelength filtering, combining or generally wavelength-division multiplexing and de-mux.

The essence is spatially periodic structure. We will see that the concept of periodic structure can be extended to photonic bandgap materials or photonic crystal: essentially 3 D grating.

Do not think “grating” as some specific device structure. Here we want to think more generally “periodic structure”. There is some general principle applied to all of these structures, and hence, while they can look very different and work in different configuration, we want to study the principle behind periodic structure.

1.2 Bragg grating and BG waveguide

Bragg grating is a periodic structure applied on a waveguide. The periodic parameter can be spatial structure (modulation), index-refraction modulation. Effectively, all BG structures can be approximated as a periodic dielectric function, to the first order.
Why Bragg grating is interesting?
- A very useful structure for many waveguide device:
- Function as filter, coupler, spectral shaper with applications in telecom (certainly), lasers, sensors... 
In short, like any periodic structure, it has application on wavelength-based devices. 
Please see presentations on Bragg grating application for what it can do.
This is a simplified introduction to the Bragg grating - it's about how it works.

1.3 Surface diffraction grating

Review from Diffraction Chapter

Below is an array of slit.

Graphics[{Rectangle[{-3,0},{43,10}],Table[{White,Rectangle[{i,0},{i+0.4,10}]},{i,1,40}] }]

Manipulate[ a=afraction*d ; k=2 π/λ ;
Plot[(((2 Sin[a k Sin[θ]])/(k Sin[θ]))((Sin[(((Ns+1)/2)) k d Sin[θ]])/(Sin[(1/2)k d Sin[θ]])))^2,{θ,-π/2,π/2}, PlotRange-> All],{d,1.2,3},{λ,0.8,1.1},{Ns,1,30,1},{afraction,0.05,0.4}]

It can be designed for transmission or reflection applications.

DiffractDiagram

2. Theoretical background: Floquet-Bloch theorem (1-D)

2.1 Summary of Floquet -Bloch theorem for 1-D (can be extended to 3-D)

Suppose we have a homogeneous linear differential equation (DE) of (n^th) order with a "potential energy" term V[x] that is periodic:
                   (Subscript[a, n]) ((∂_x)^n) ψ[x] +(Subscript[a, n-1]) ((∂_x)^(n-1)) ψ[x] + ... + V[x] ψ [x] =0            (2.1.1)
Does the solution have some thing "periodic" in its property?
The answer is “yes”, based on Floquet-Bloch theorem, which, as we known, crucial for the band theory in crystalline condensed matter physics.
The theorem states that any solution can be written as:
                            ψ[x]=(e^(i k x)) (Subscript[u, k])[x]                 (2.1.2)
where k is a certain characteristic value of the DE (which doesn’t have to be real, it can be a general complex number), and (Subscript[u, k])[x] is a periodic function.

Example: of 1st order DE

DSolve[ψ'[x] +((Cos[x])^2)  ψ[x]== 0,ψ[x],x]

We see that the solution can be written in the form:
               ψ(x)=(Subscript[c, 1]) (e^(-(x/2)))(e^(-(1/4) sin(2 x)))=(Subscript[c, 1]) (e^(-(x/2)))u[x]       
where (e^(-(x/2))) is the part “(e^(i k x))” with characteristic value  k=i(1/2) and the periodic part is: u[x]=(e^(-(1/4) sin(2 x))). They can be plotted:

Plot[{(E^(-(x/2))), (E^(- (1/4)Sin[2 x])) ,(E^(-(x/2))) * (E^(- (1/4)Sin[2 x]))},{x,0,10}, PlotRange-> All]

2.2 Discussion

We don’t need any formal proof to understand and use FB theorem (there are various proofs; and for Bloch theorem, which is a case of Floquet’s theorem with Schroedinger equation, the proof is given in most condensed matter theory texts). But it is useful to gain some insight to see how it comes about. (Skip if it appears too “taxing” for now)

Let  ψ[x] be a solution of (2.1.1). If we define a new function:
                               φ[x]≡ ψ[x+T]                                       (2.1.2)
where T is the period, it can be shown that it is also a solution. First: 
    (Subscript[a, n]) ((∂_x)^n) φ[x] +(Subscript[a, n-1]) ((∂_x)^(n-1)) φ[x] + ... + V[x] φ [x] =
	  (Subscript[a, n]) ((∂_x)^n) ψ[x+T] +(Subscript[a, n-1]) ((∂_x)^(n-1)) ψ[x+T] + ... + V[x] ψ[x+T]    (2.1.3)
Since:                                 V[x]=V[x+T]                  (2.1.4)
    (Subscript[a, n]) ((∂_x)^n) ψ[x+T] +(Subscript[a, n-1]) ((∂_x)^(n-1)) ψ[x+T] + ... + V[x] ψ[x+T] 
=  (Subscript[a, n]) ((∂_x)^n) ψ[x+T] +(Subscript[a, n-1]) ((∂_x)^(n-1)) ψ[x+T] + ... + V[x+T] ψ[x+T]=0   (2.1.5)
and thus:
                  (Subscript[a, n]) ((∂_x)^n) φ[x] +(Subscript[a, n-1]) ((∂_x)^(n-1)) φ[x] + ... + V[x] φ[x] =0.              (2.1.6)
which shows that φ[x]≡ ψ[x+T]  is indeed a solution (but not necessarily the same as ψ[x]).

In the example above, it is first order DE, hence there can be ONLY ONE independent solution. Thus:
                        φ[x]≡ψ[x+T]=α  ψ[x]          (2.1.7)
In general:             ψ[x+m T]=(α^m)  ψ[x]        (2.1.8)   (m: arbitrary integer).
Define a new function:       u[x]=(e^(-i λ x))  ψ[x] , which means:      ψ[x]=(e^(i λ x)) u[x]        (2.1.9)
where λ is some value to be determined later.
Then:                   ψ[x+m T]=(e^(i m  λ T))(e^( i λ x))  u[x+m T]           (2.1.10)
But:                     ψ[x+m T]=(α^m)  ψ[x]=(α^m)  (e^(i λ x)) u[x]          (2.1.11)
If we choose:               (e^(i λ T))=α                                             (2.1.12)
then, comparing  (2.1.10) and (2.1.11), using  (2.1.12):
                 (e^( i λ x))  u[x+m T]=(e^(i λ x)) u[x]                      (2.1.13)
or:                       u[x+m T]= u[x]                         (2.1.14)
which means that u[x] is a periodic function. QED

Note that the above proof can be extended for arbitrary nth-order DE by just using a vector of basis functions instead of a single function:
Let  Ψ[x]=(({{(Subscript[ψ, 1])[x ]}, {(Subscript[ψ, 2])[x ]}, {⋮}, {(Subscript[ψ, n])[x ]}})) be a set of linearly indep. solutions. Since   Ψ[x+T]=(({{(Subscript[ψ, 1])[x +T]}, {(Subscript[ψ, 2])[x+T ]}, {⋮}, {(Subscript[ψ, n])[x+T ]}}))  is also a set of solutions, we have:
         Ψ[x+T]=(({{(Subscript[ψ, 1])[x +T]}, {(Subscript[ψ, 2])[x+T ]}, {⋮}, {(Subscript[ψ, n])[x+T ]}}))= (({{(Subscript[c, 11]), (Subscript[c, 12]), ⋯, (Subscript[c, 1n])}, {(Subscript[c, 21]), (Subscript[c, 22]), ⋯, (Subscript[c, 2n])}, {⋮, ⋮, ⋱, ⋮}, {(Subscript[c, n 1]), (Subscript[c, n 2]), ⋯, (Subscript[c, n n])}})) .  (({{(Subscript[ψ, 1])[x ]}, {(Subscript[ψ, 2])[x ]}, {⋮}, {(Subscript[ψ, n])[x ]}}))     (2.1.15a)
         Ψ[x+T]=M .  Ψ[x]                 (2.1.15b)

Let A be the linear transformation matrix that diagonalize M:
                            A. M .(A^(-1)) =(({{(e^(Subscript[λ, 1])), 0, ⋯, 0}, {0, (e^(Subscript[λ, 2])), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(Subscript[λ, n]))}})) ≡ Λ         (2.1.16a)    
                           (A^(-1)) Λ  A = M           (2.1.16b)    
(we can express the eigenvalues in exponential form because no eigenvalue can vanish for a complete basis set of linearly independent solutions). Then, 
         Ψ[x+T]=M .  Ψ[x]=(A^(-1)) Λ  A  Ψ[x]                 (2.1.17a)
         A Φ[x+T]=Λ  A  Φ[x]                    (2.1.17b)
Define:
         A . (({{(e^(-i (Subscript[k, 1])x)), 0, ⋯, 0}, {0, (e^(-i (Subscript[k, 2])x)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(-i (Subscript[k, n])x))}})) . (({{(Subscript[ψ, 1])[x ]}, {(Subscript[ψ, 2])[x ]}, {⋮}, {(Subscript[ψ, n])[x ]}}))≡U[x]≡ (({{(Subscript[u, 1])[x ]}, {(Subscript[u, 2])[x ]}, {⋮}, {(Subscript[u, n])[x ]}}))     (2.1.18a)
or:           A (({{(Subscript[ψ, 1])[x ]}, {(Subscript[ψ, 2])[x ]}, {⋮}, {(Subscript[ψ, n])[x ]}}))=  (({{(e^(i (Subscript[k, 1])x)), 0, ⋯, 0}, {0, (e^(i (Subscript[k, 2])x)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(i (Subscript[k, n])x))}}))  (({{(Subscript[u, 1])[x ]}, {(Subscript[u, 2])[x ]}, {⋮}, {(Subscript[u, n])[x ]}}))     (2.1.18b)
where (Subscript[k, i]), i=1, n are to be determined.
           A (({{(Subscript[ψ, 1])[x +T]}, {(Subscript[ψ, 2])[x+T ]}, {⋮}, {(Subscript[ψ, n])[x +T]}}))=  (({{(e^(i (Subscript[k, 1])T)), 0, ⋯, 0}, {0, (e^(i (Subscript[k, 2])T)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(i (Subscript[k, n])T))}})) (({{(e^(i (Subscript[k, 1])x)), 0, ⋯, 0}, {0, (e^(i (Subscript[k, 2])x)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(i (Subscript[k, n])x))}}))  (({{(Subscript[u, 1])[x+T ]}, {(Subscript[u, 2])[x +T]}, {⋮}, {(Subscript[u, n])[x+T ]}}))     (2.1.19)
Substitute (2.1.19) and (2.1.18b) in (2.1.17b), we obtain:
       (({{(e^(i (Subscript[k, 1])T)), 0, ⋯, 0}, {0, (e^(i (Subscript[k, 2])T)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(i (Subscript[k, n])T))}}))(({{(e^(i (Subscript[k, 1])x)), 0, ⋯, 0}, {0, (e^(i (Subscript[k, 2])x)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(i (Subscript[k, n])x))}})) (({{(Subscript[u, 1])[x +T]}, {(Subscript[u, 2])[x +T]}, {⋮}, {(Subscript[u, n])[x +T]}})) 
=(({{(e^(Subscript[λ, 1])), 0, ⋯, 0}, {0, (e^(Subscript[λ, 2])), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(Subscript[λ, n]))}})) (({{(e^(i (Subscript[k, 1])x)), 0, ⋯, 0}, {0, (e^(i (Subscript[k, 2])x)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(i (Subscript[k, n])x))}}))  (({{(Subscript[u, 1])[x ]}, {(Subscript[u, 2])[x ]}, {⋮}, {(Subscript[u, n])[x ]}}))     (2.1.20)
  Clearly, we want to define the k’s values such that:
       (({{(e^(i (Subscript[k, 1])T)), 0, ⋯, 0}, {0, (e^(i (Subscript[k, 2])T)), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(i (Subscript[k, n])T))}}))=(({{(e^(Subscript[λ, 1])), 0, ⋯, 0}, {0, (e^(Subscript[λ, 2])), ⋯, 0}, {⋮, ⋮, ⋱, ⋮}, {0, 0, ⋯, (e^(Subscript[λ, n]))}}))      (2.1.21)  
  which leads to: 
        (({{(Subscript[u, 1])[x +T]}, {(Subscript[u, 2])[x +T]}, {⋮}, {(Subscript[u, n])[x +T]}})) = (({{(Subscript[u, 1])[x ]}, {(Subscript[u, 2])[x ]}, {⋮}, {(Subscript[u, n])[x ]}}))       (2.1.22a)
Or:                 U[x+T]=U[x]    (2.1.22b)   QED

2.2.1 Application

There are 2 things about this theorem. First:

Note that the theorem doesn’t tell us what U[x] is, and of course we still have to solve for it. Same with the characteristic coefficient k’s; They are specific to the DE and the “potential energy” term V[x] itself, and we have to solve the DE to find out. We can take k as unknown, solve for U[k;x] with the requirement that it is periodic, then it becomes a condition on k. If U[k;x] is a physical entity, such as the E field, then additional physical conditions can be imposed to completely determine physically valid k’s and U[k;x].

Second, one may ask, what is the use of this theorem, if we already know the solution by solving the DE directly which can be more convenient than solving for U[x]? For the example 2.1.1 above:

(∂_x)( (E^(I k x)) u[x] )+((Cos[x])^2)  (E^(I k x)) u[x]

DSolve [u'(x)+I k  u(x)+ u(x) (cos^2)(x) ==0, u[x],x]

We note that in the end, u(x)-> (Subscript[c, 1]) (e^(-i k x-(x/2)-(1/4) sin(2 x))) contains the term (e^(-i k x)) that cancels out the exponent term (e^(i k x )): 
                    ψ[x]=(e^(i k x ))u(x)= (Subscript[c, 1]) (e^(i k x )) ((e^(-i k x-(x/2)-(1/4) sin(2 x))))=(Subscript[c, 1]) (e^(-(x/2)-(1/4) sin(2 x)))   
which can be obtained by solving a simpler DE than that of u[x]. Infact, it is obvious that if we happen to know the solution ψ[x], then we don’t have to bother to solve for u[x] because:
                   u[x]=(e^(-i k x ))ψ[x]     
and we only have to impose the condition in order for u[x] to be periodic:   
                u[x]=(e^(-i k x-(x/2)-(1/4) sin(2 x)))=u[x+ π]=(e^(-(i k +(1/2))  π))(e^(-i k x-(x/2)-(1/4) sin(2 x)))  
The requirement is:                       (e^(-(i k +(1/2)) π))=1
This equation can only be satisfied with: i k +(1/2)=2 m i

Simplify[Solve[I k +(1/2)==2 m I,k] ]

Hence: u(x)=(e^(-i (2 m+(i/2)) x-(x/2)-(1/4) sin(2 x)))=(e^(-2 i m x-(1/4) sin(2 x)))

Formally:
         ψ[x]=(e^(i k x ))u(x)= (Subscript[c, 1]) (e^(i (2 m+(i/2)) x)) ((e^(-2 i m x-(1/4) sin(2 x))))

Again, we see that the term (e^(i 2 m x))  appears irrelevant as it cancels out.

3. Floquet-Bloch theorem for Bragg structure in 1-D Helmholtz

3.1 Bragg waveguide grating

We will find Bragg structure such as grating or dielectric thin-film stack involves the case of 1-D periodic “potential” where FB theorem can be applied. Strictly speaking, we have 3D equation:  

where δϵ(r) is a very small function compared with (Subscript[ϵ, 0]). The variation of ϵ[x,y,z] vs. z is so mall compared with {x,y}, which is the waveguide dielectric profile, hence we can use the perturbation approach: With this approach, we simply assume that the solution is the x-y dimension is simply that of the waveguide. But using these solutions, we find what how the different modes are “mixed” because of this perturbation. This is the concept of couple-mode approach

3.2 Coupled-mode approach for waveguides

Neglecting the small z-dependence, we use the wave solution E(x,y) (e^(i β z)) with (e^(i β z)) being an eigensolution of the equation. But with z-dependence ϵ, we obviously cannot use a single β. Intuitively, the wave solution must contain several β. The perturbation approach is to expand the wave solution into eigencomponents. So, we can use the function

where:

The meaning of this function is that the z-dependence dielectric term can result in many modes mixed together, including both bound and radiative states. For simplicity, let's first consider only the bound state.
The electric field is a linear combination:
                       E[x,y,z] = (∑_(m=-N)^N)(Subscript[a, m])[z]  (Subscript[E, m])[x,y] (e^(i (Subscript[β, m]) z))  
We know  (Subscript[E, m])[x,y] (e^(i (Subscript[β, m]) z)), we only have to solve for (Subscript[a, m])[z]

After extensive simplication by substitution, we obtain the equation for (Subscript[a, n])[z]:
(Subscript[a, n]^′)(z) = ((i  (ω^(2 )))/(2 (c^2) (Subscript[β, n]))) (∑_(m=1)^N)(e^(i z ((Subscript[β, m])-(Subscript[β, n])))) (∫∫δϵ(x,y,z) (Subscript[Ef, m])(x,y) (SuperStar[((Subscript[Ef, n]))])(x,y)dxdy) (Subscript[a, m])(z)

This equation shows that the perturbation δϵ couples the modes together.
It means that the amplitude of one mode depends on the others. That even if we launch some light in a pure mode, it won't remain in that mode, but gradually gains property of other modes.
This is similar to a quantum mechanical system in which a perturbation can couple the various eigenmodes together.  That' s the reason this equation is called coupled - mode equation. Here the variable is z instead of t as in QM system

3.3 Single-mode waveguide:

If there is only one  β, in other words, single-mode waveguide, how do we “couple” single mode together?
Actually, there are always 2 modes: +β and -β: In other words: forward and backward modes.

Substituting, we have this set of equations:

where (Subscript[G, K]), (Subscript[G, -K]) are the coupling coefficient

3.3.1 Solution to coupled mode equations

The general solution for electric field is a linear combination of a forward and backward component: 
                                         E=(Subscript[a, F]) (e^(i β z)) + (Subscript[a, B]) (e^(-i β z))
Or:
                                          E=(({{(e^(i β z)),  (e^(-i β z))}})) (({{(Subscript[a, F])}, {(Subscript[a, B])}})) 
where:                                   (({{(Subscript[a, F])}, {(Subscript[a, B])}}))= M[z] (({{(Subscript[C, 1 ])}, {(Subscript[C, 2])}})) 
the matrix M is
                  M[z]=(({{(e^(i z (Sqrt[((Δβ)^2)-G (OverBar[G])])-i z Δβ)), -(((e^(-i Δβ z-i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/(OverBar[G]))}, {-(((e^(i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/G), (e^(i z Δβ-i z (Sqrt[((Δβ)^2)-G (OverBar[G])])))}}))
where                             Δβ=β-K/2, K is 2 π/Λ, the grating period
and  (({{(Subscript[C, 1 ])}, {(Subscript[C, 2])}}))  is an arbitrary constant.
Example let (({{(Subscript[C, 1 ])}, {(Subscript[C, 2])}}))=(({{1}, {0}}))

M[z_]:=(({{(E^(I z (Sqrt[((Δβ)^2)-G (OverBar[G])])-I z Δβ)), -(((E^(-I Δβ z-I (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/(OverBar[G]))}, {-(((E^(I Δβ z+I (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/G), (E^(I z Δβ-I z (Sqrt[((Δβ)^2)-G (OverBar[G])])))}}))

E1 = (({{(E^(I β z)), (E^(-I β z))}}))  .  M[z]  . (({{1}, {0}}))

Hence the electric field is:
               (Subscript[E, 1])=(e^(i β z-i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z))  -(e^(-i β z+i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (( (Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/G) 
What happens if there is no grating? then G (OverBar[G])-> 0 and
                                                 (Subscript[E, 1])=(e^(i β z)) 
which is the forward wave.
Likewise, (({{(Subscript[C, 1 ])}, {(Subscript[C, 2])}}))=(({{0}, {1}})) will lead to a backward moving wave as expected.

E1field[λ_,n_, Λ_,G_ , z_]:=Module[{β, K, Δβ, Q},
         β= 2 n π /λ ;    K= 2 π /Λ  ;
             Δβ=β- K/2  ;
         Q=If[((Δβ)^2)-Abs[G]^2>0, Sign[Δβ]  (Sqrt[((Δβ)^2)-Abs[G]^2]) ,
I  (Sqrt[Abs[((Δβ)^2)-Abs[G]^2 ]])  ]   ;

(E^(I β z-I Δβ z+I Q z)) -(E^(-I β z+I Δβ z+I Q z)) (( (Δβ-Q))/G)  
];
E1Bfield[λ_,n_, Λ_,G_ , z_]:=Module[{β, K, Δβ, Q},
         β= 2 n π /λ ;    K= 2 π /Λ  ;
             Δβ=β- K/2  ;
         Q=If[((Δβ)^2)-Abs[G]^2>0, Sign[Δβ]  (Sqrt[((Δβ)^2)-Abs[G]^2]) ,
I  (Sqrt[Abs[((Δβ)^2)-Abs[G]^2 ]])  ]   ;
 -(E^(-I β z+I Δβ z+I Q z)) (( (Δβ-Q))/G)
]  ;
E1Ffield[λ_,n_, Λ_,G_ , z_]:=Module[{β, K, Δβ, Q},
         β= 2 n π /λ ;    K= 2 π /Λ  ;
             Δβ=β- K/2  ;
         Q=If[((Δβ)^2)-Abs[G]^2>0, Sign[Δβ]  (Sqrt[((Δβ)^2)-Abs[G]^2]) ,
I  (Sqrt[Abs[((Δβ)^2)-Abs[G]^2 ]])  ]   ;
 (E^(I β z-I Δβ z+I Q z))  
]  ;

The example below will show what happens to the wave when it hits "resonance"

λ=1.55 ; n=1.47 ; Λ= 1.55 /(2* n)  ; G=0.05  ;
Manipulate[
Animate[Plot[{Re[ E1Ffield[λ,n, Λ,G , z] * Exp[-I  2 Pi  t] ]
,Re[ E1Bfield[λ,n, Λ,G , z] * Exp[-I  2 Pi  t] ]
,Re[ E1field[λ,n, Λ,G , z] * Exp[-I  2 Pi  t]  ]} ,{z, 0,20}, AspectRatio-> 0.2, PlotRange-> {-2,2}
, PlotStyle-> {Red, Blue, Black}] 
,{t,0,1}, AnimationRunning-> False]
,  {λ,1.51,1.58}  ]

The grating "resonance" (a lose term) happens when ((Δβ)^2)-G (OverBar[G])  <0
Or               |Δβ|<(Sqrt[G (OverBar[G])])
which is:       2 π |(n/λ)-(1/(2 Λ ))|<(Sqrt[G (OverBar[G])]) 
In other words, when the wavelength on the waveguide is very close to 2 times the grating period. Then the wave has the exponential form as we see. Since this exponential form cannot be at infinite to one end, this solution is not possible everywhere. This is the "forbidden" band, where no waves can propagate forever. A wave that enters will be reflected.
                                 (Subscript[E, 1])=(e^(i β z-i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z))  -(e^(-i β z+i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (( (Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/G) 
At the two edges of the band, the solution is:
                                 (Subscript[E, 1])=(e^(i β z-i Δβ z))  -Sign[Δβ]  (e^(-i β z+i Δβ  z))  
                                (Subscript[E, 1])=(e^(i K z/2))  -Sign[Δβ]  (e^(-i K z/2))
which shows the propagation constant is just half of that of the grating. This is true everywhere inside the band since (Sqrt[((Δβ)^2)-G (OverBar[G])]) is imaginary. We have the interesting result that the dispersion curve k-w is flat completely!!.

3.3.2 DFB reflection coefficient

Let's look at the solution:
 ({{E1=(e^(i β z-i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z))      -    (e^(-i β z+i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (((Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/G)}}) obtained by choosing (({{1}, {0}})) input; and 
E2=-(((Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/(OverBar[G]))((e )^(i β z-i Δβ z-i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) +  (e^(-i β z+i Δβ z-i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z))  for input (({{0}, {1}})).
The waves can NOT just go in one direction! Each has a forward component and a companion backward component. They co-exist together. What this means physically is that if we have a forward wave (e^(i β z-i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)), then there must be a REFLECTED wave -(e^(-i β z+i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (((Δβ-(Sqrt[((Δβ)^2)-G (OverBar[G])])))/G). Thus, we can define a coefficient:
(Subscript[r, B]) =(((Sqrt[((Δβ)^2)-G (OverBar[G])]) -Δβ   )/G)=  - (( (OverBar[G]))/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))  
as the "effective" reflection coefficient of the grating for the forward wave. Similarly,
we define:
(Subscript[r, F]) =(((Sqrt[((Δβ)^2)-G (OverBar[G])]) -Δβ   )/(OverBar[G]))=-(( G)/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))
as the effective reflection coefficient of the backward wave.
If the grating disappear, (Subscript[r, B]) and (Subscript[r, F]) indeed go to zero as expected.
Notice that:
(Subscript[r, B]) (Subscript[r, F]) = ((Δβ- (Sqrt[((Δβ)^2)-G (OverBar[G])    ]))/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))

How "big" is r? The interesting case is the inband case: ((Δβ)^2)-G (OverBar[G]) < 0. Then |r| = 1. This means that the "reflected" wave is as strong as the incident wave in the band. 
Let's assume for simplicity that Δβ = β-K/2 is real. G (OverBar[G]) is also real and of the dimension 1/length. Thus, we can plot as follow:

rF[Δβ_, ϕ_, G_]:= ((Sign[Re[Δβ]] *(Sqrt[((Δβ)^2)-Abs[G]^2])- Δβ)/(G*(E^(I ϕ)))) ;
(* Notice that we must have the Sign[Re[Δβ]] included to avoid the  *)
(* erroneous condition if we change sign of Δβ, as mentioned above *)
Plot[Abs[rF[x,0.,1.]],{x,-5,5},
	       PlotStyle->{RGBColor[1.,0.,0],Thickness[0.01]}
				, PlotRange->{0.,1}, GridLines->Automatic]

Plot[{ Arg[rF[x, 0.,1.]],Arg[rF[x,0.75,1]]}, {x,-2,2},
	       PlotStyle->{{RGBColor[1.,0.,0],Thickness[0.01]}
	                                  ,{RGBColor[0.,0.,1],Thickness[0.01]}}
				, PlotRange->{-4,4}, GridLines->Automatic]

The effect of a complex-coupled grating is a phase shifted difference between the two reflection coefficient. Notice the discontinuity of π at Δβ=0 is also caused by the arbitrary numerical function plotting. 
When Δβ is complex (gain or loss), the reflection coefficient deviates from 1:

Plot[Release[Table[ Abs[rF[x-I*j,0.,1.]],{j,0.,1.,0.2}]],{x,-5,5},
	       PlotStyle->{{RGBColor[1.,0.,0],Thickness[0.01]}
		                                 }
				, PlotRange->{0.,1}, GridLines->Automatic]

Plot[Release[Table[ Abs[rF[x-I*1,0.,j]],{j,0.4,1.,0.2}]],{x,-5,5},
	       PlotStyle->{{RGBColor[1.,0.,0],Thickness[0.01]}
		                                 }
				, PlotRange->{0.,0.5}, GridLines->Automatic]

We will see that in a DFB laser, these terms (Subscript[r, F]) (Subscript[r, B]) play the role of end mirror reflectivity, except that they are distributed along the length of the grating, not from some point. Hence, that's where the concept of distributed feedback (DFB) comes from.
In terms of (Subscript[r, B]) and (Subscript[r, F]) coefficients, we can express the M matrix as:  
M[z_]:=(({{(e^(i z (Sqrt[((Δβ)^2)-G (OverBar[G])])-i z Δβ)), (Subscript[r, B])  (e^(-i Δβ z-i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z ))}, {(Subscript[r, F])  (e^(i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z )), (e^(i z Δβ-i z (Sqrt[((Δβ)^2)-G (OverBar[G])])))}}))

3.3.3 Propagation constant

We see that the propagation in the BG section is modified from the unperturbed section by a shift  -Δβ  ± (Sqrt[((Δβ)^2)-G (OverBar[G]) ]). The forward propagation constant is:
β-Δβ  ± (Sqrt[((Δβ)^2)-G (OverBar[G])])

β-(Δβ±(Sqrt[((Δβ)^2)-G (OverBar[G])]))  /.  (* substitute *) Δβ->  β-K/2

or:  (K/2) ∓ (Sqrt[(((β-(K/2)))^2)-G (OverBar[G])])

Let G -> 0, i. e. there is no grating:

(K/2) ∓ (Sqrt[(((β-(K/2)))^2)-G (OverBar[G])])                    /.                   G->0

We see two propagation constants for forward mode:         β  and       -β +K
The second one does not make sense! Without a grating, why would the system has anything to do with K? Actually, if we go back to the full solution,
(({{(Subscript[a, F])}, {(Subscript[a, B])}})) = (Subscript[C, 1]) (({{(e^(i z (Sqrt[((Δβ)^2)-G (OverBar[G])])-i z Δβ))}, {-(((e^(i Δβ z+i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (OverBar[G]))/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))}}))+(Subscript[C, 2]) (({{-(((e^(-i Δβ z-i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) G )/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))}, {(e^(i z Δβ-i z (Sqrt[((Δβ)^2)-G (OverBar[G])]) ))}}))
we see that the second term, -(((e^(-i Δβ z-i (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) G )/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])]))), vanishes for G->0, so the propagation constant (-β +K) is irrelevant. In fact,

(({{(Subscript[a, F])}, {(Subscript[a, B])}})) = (Subscript[C, 1]) (({{(E^(I z (Sqrt[((Δβ)^2)-G (OverBar[G])])-I z Δβ))}, {-(((E^(I Δβ z+I (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (OverBar[G]))/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))}}))+(Subscript[C, 2]) (({{-(((E^(-I Δβ z-I (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) G )/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))}, {(E^(I z Δβ-I z (Sqrt[((Δβ)^2)-G (OverBar[G])]) ))}})) /. (OverBar[G])->0 /. G->0

PowerExpand[%]

We see that the solution gets back to their unperturbed form! So we know that this solution has proper limit.

3.3.4 Stop band

What happens if  ((Δβ)^2)-G (OverBar[G]) < 0? Now (Sqrt[((Δβ)^2)-G (OverBar[G])]) is imaginary, which we can write
 (Sqrt[((Δβ)^2)-G (OverBar[G])]) = i Γ, 
 the solution has an exponential component:

(({{(Subscript[a, F])}, {(Subscript[a, B])}})) = (Subscript[C, 1]) (({{(E^(I z (Sqrt[((Δβ)^2)-G (OverBar[G])])-I z Δβ))}, {-(((E^(I Δβ z+I (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) (OverBar[G]))/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))}}))+(Subscript[C, 2]) (({{-(((E^(-I Δβ z-I (Sqrt[((Δβ)^2)-G (OverBar[G])]) z)) G )/(Δβ+(Sqrt[((Δβ)^2)-G (OverBar[G])])))}, {(E^(I z Δβ-I z (Sqrt[((Δβ)^2)-G (OverBar[G])]) ))}})) /. (Sqrt[((Δβ)^2)-G (OverBar[G])])->I Γ

This is interesting, as the wave solution now has an exponential amplitude. The forward wave is  
((e^(-z Γ-i z Δβ)) (Subscript[C, 1])-(((e^(z Γ-i z Δβ)) G (Subscript[C, 2]))/(i Γ+Δβ))) (e^(i β z))

(E^(I z β-I z Δβ)) ((E^(-z Γ)) (Subscript[C, 1])-(((E^(z Γ)) G (Subscript[C, 2]))/(I Γ+Δβ))) /.  Δβ->β-K/2

FullSimplify[%]

And for the backward wave:
(e^(-(1/2) i K z)) (((2 (e^(-z Γ)) (OverBar[G]) (Subscript[C, 1]))/(K-2 β-2 i Γ))+(e^(z Γ)) (Subscript[C, 2])).

The wave now propagates with a propagation constant that is of the grating! The exponential amplitude has a characteristic length ~ 1/Γ ≡ 1/(Sqrt[|G((|)^2) -  ((Δβ)^2)]). The larger the modulation depth of the grating, the stronger the exponential decay or growth. If a wave exponentially decays, where does its loss energy go? if it exponentially gains, where does it energy gain come from? Obviously, there are only two waves we consider here, forward and backward. The energy must simplify be exchanged between them. So the decaying energy of the forward wave must be the one adding to the exponential gain of the backward wave and vice versa. This is the reason why they are called coupled waves. In essence, the grating scatters the forward wave into the backward wave and vice versa. More generally, the grating "diffracts" any wave into many modes. Some of these modes are radiative which we have not included in our equation. Here we consider only the special case that a strong diffraction is the coupling between the forward and backward mode of a grating.

Can the wave propagate to infinity? Obviously not, because the waves: (e^((i K z)/2)) ((e^(-z Γ)) (Subscript[C, 1])+((2 (e^(z Γ)) G )/(K-2 β-2 i Γ))(Subscript[C, 2])) or  (e^(-(1/2) i K z)) (((2 (e^(-z Γ)) (OverBar[G]) (Subscript[C, 1]))/(K-2 β-2 i Γ))+(e^(z Γ)) (Subscript[C, 2]))
become infinite at z -infinity or z -> +infinity. Such a solution can not exist in an infinitely-long Bragg grating. The band  (((Δβ)^2)-G (OverBar[G]))<0  is not a physical solution, therefore they are called stop-band. It is equivalent with an energy band gap (some call it "photonics crystal" ). For finite length, it is of course a proper solution and as we can see later that it is part of the BG reflection and transmission band.

3.3.5 Plot of the dispersion curve

We can plot the propagation constant  to see what its like(K/2) ∓ (Sqrt[(((β-(K/2)))^2)-G (OverBar[G])])

dplus[K_,β_,G_, ω_]:=(K/2) + Re[(Sqrt[(((β[ω]-(K/2)))^2)-Abs[G]^2])  ];
dminus[K_,β_,G_, ω_]:=(K/2) - Re[(Sqrt[(((β[ω]-(K/2)))^2)-Abs[G]^2 ]) ];

Testbeta[x_]:=x-0.1 *(x-1)^2
Plot[Testbeta[x],{x,0.,2}]

Plot[{dplus[2.,Testbeta,0.2, ω],
		          dminus[2.,Testbeta,0.2, ω]},{ ω,0.5,2},
             AxesOrigin->{0.,0}, PlotRange->{{0.,2},{0,2}}]

This is a very interesting result. When the propagation constant β is within the grating modulation band region, the wave propagates with a constant wave number: no dispersion.
We can see that ((dω)/dk) -> 0 from the outside of the zone boundary where the exponential wave happens. Does this mean that the group velocity -> 0 at the zone boundary? On the other hand, inside the zone boundary ((dω)/dk) -> ∞, does this means that group velocity -> ∞? Obviously, the concept of group velocity does not mean the same thing here! The dispersion at this point is too large for the group velocity concept relevant.

3.4 Relationship with transfer S matrix approach

The F-B theorem can also be shown with the transfer matrix approach. Applying the transfer matrix method at the first boundary of any segment, which we can set as x=0:
                (Subscript[(Subscript[(({{F}, {F'}}))|, x=Λ])=S (({{F}, {F'}}))|, x=0])    (3.4.1)
For the  (n^th) segment (which can be left or right of the point x=0, with n>0 or <0):
       (({{(Subscript[F, n])}, {(Subscript[F, n])'}}))=(S^n) (({{(Subscript[F, 0])}, {(Subscript[F, 0])'}}))                (3.4.2)     
because all segments are the same, the transfer S matrix is the same. In fact, within each segment, e. g. the (n^th) segment, the E field is given by:       
      (({{F[x]}, {F'[x]}}))=S[ξ;0] . (({{(Subscript[F, n-1])}, {(Subscript[F, n-1])'}}))  ; where  ξ=x -(n-1) Λ ∈ [0,Λ]   (3.4.3a)
or more generally:
         ξ=Modulo[x ,Λ ]         (3.4.3b)
as all we need to know is (({{(Subscript[F, n-1])}, {(Subscript[F, n-1])'}})). Note that the variable ξ is relative within each segment, and S[ξ;0] is the same for all segment, which will be seen later that it is the periodic part of the F-B theorem. Note that in this notation:
                 S=S[Λ;0]           (3.4.4)
In the previous HW on S-matrix, we already discuss “propagating” solutions where (({{F}, {F'}})) at a point is the eigenvector of S. Let’s express:
      S . (({{((Subscript[c, 1])^((+))), ((Subscript[c, 1])^((-)))}, {((Subscript[c, 2])^((+))), ((Subscript[c, 2])^((-)))}}))=(({{(e^(i (ϕ^((+))))), 0}, {0, (e^(i (ϕ^((-)))))}})) . (({{((Subscript[c, 1])^((+))), ((Subscript[c, 1])^((-)))}, {((Subscript[c, 2])^((+))), ((Subscript[c, 2])^((-)))}}))     (3.4.5)
where (e^(i (ϕ^((+))))), (e^(i (ϕ^((-))))) are the eigenvalues, and (ϕ^((+))), (ϕ^((-))) are just quantities that we define them to be so and they can be real or complex - there is no restriction what they should be in this case. The column vectors (({{((Subscript[c, 1])^((+)))}, {((Subscript[c, 2])^((+)))}})), (({{((Subscript[c, 1])^((-)))}, {((Subscript[c, 2])^((-)))}})) are the corresponding eigenvectors. They are the special solutions with (({{(Subscript[F, 0])}, {(Subscript[F, 0])'}})) such that they remains the same except for a multiplication factor (e^(i (ϕ^((+))))), (e^(i (ϕ^((-))))) after each segment. 
If Det[S]=1, we always have: (ϕ^((-)))=- (ϕ^((+)))
If we choose a solution such that (({{(Subscript[F, 0])}, {(Subscript[F, 0])'}})) is an eigenvector:
 if   (({{(Subscript[F, 0])}, {(Subscript[F, 0])'}}))= (({{((Subscript[c, 1])^((+)))}, {((Subscript[c, 2])^((+)))}}))   then     (({{(Subscript[F, n])}, {(Subscript[F, n])'}}))=(S^n) (({{(Subscript[F, 0])}, {(Subscript[F, 0])'}}))=(e^(n i (ϕ^((+))))) (({{(Subscript[F, 0])}, {(Subscript[F, 0])'}})) (3.4.6a)      
or:   if   (({{(Subscript[F, 0])}, {(Subscript[F, 0])'}}))= (({{((Subscript[c, 1])^((-)))}, {((Subscript[c, 2])^((-)))}})) then   (({{(Subscript[F, n])}, {(Subscript[F, n])'}}))=(e^(n i (ϕ^((-)))))(({{(Subscript[F, 0])}, {(Subscript[F, 0])'}}))                     (3.4.6b)

It follows, by substituting in (3.4.3) that for the (n^th) segment:
      (({{F[x]}, {F'[x]}}))=(e^((n-1) i (ϕ^((±)))))S[ξ;0] .  (({{((Subscript[c, 1])^((±)))}, {((Subscript[c, 2])^((±)))}}))             (3.4.7a)
Or:
     F[x]=(e^((n-1) i (ϕ^((±))))) (((Subscript[c, 1])^((±)))(Subscript[S, 11])[ξ;0] +((Subscript[c, 2])^((±)))(Subscript[S, 12])[ξ;0])        (3.4.7b)
Thus, if we define a function u[x] to be:
  u[x]=(e^(-i q x ))F[x]=(e^(-i q x )) (e^((n-1) i (ϕ^((±))))) (((Subscript[c, 1])^((±)))(Subscript[S, 11])[ξ;0] +((Subscript[c, 2])^((±)))(Subscript[S, 12])[ξ;0])       (3.4.8a)
for x belonging to the (n^th) segment:
              x  ∈ [(n-1) Λ,n Λ]          (3.4.8b)
with the periodic requirement: u[x+Λ]=u[x], 
then:   (e^(-i q (x+Λ) )) (e^(n i (ϕ^((±))))) (((Subscript[c, 1])^((±)))(Subscript[S, 11])[ξ;0] +((Subscript[c, 2])^((±)))(Subscript[S, 12])[ξ;0])=
               (e^(-i q x )) (e^((n-1) i (ϕ^((±))))) (((Subscript[c, 1])^((±)))(Subscript[S, 11])[ξ;0] +((Subscript[c, 2])^((±)))(Subscript[S, 12])[ξ;0])     (3.4.8c)
We see that the part (((Subscript[c, 1])^((±)))(Subscript[S, 11])[ξ;0] +((Subscript[c, 2])^((±)))(Subscript[S, 12])[ξ;0]) is automatically periodic because ξ is defined as in (3.4.3b). Thus, the periodic requirement is that:
              (e^(-i q Λ ))=(e^(-i (ϕ^((±)))))                            (3.4.9a)
or we choose to define:
                 (q^((±)))=((ϕ^((±)))/Λ)                                  (3.4.9b)
What (3.4.9b) means is that we can indeed find the periodic function  u[x]=(e^(-i q x ))F[x] with q given in (3.4.9b).
Note: it is our choice to pick these values for (q^((±))), and here, we have 2 choices.
Thus, if we define 2 functions:
  (u^((±)))[x]=(e^(-i (q^((±))) x )) (e^(i (q^((±))) (n-1)  Λ )) (((Subscript[c, 1])^((±)))(Subscript[S, 11])[ξ;0] +((Subscript[c, 2])^((±)))(Subscript[S, 12])[ξ;0])   (3.4.10a) 
or:
   (u^((±)))[x]=(e^(-i (q^((±))) ξ )) (((Subscript[c, 1])^((±)))(Subscript[S, 11])[ξ;0] +((Subscript[c, 2])^((±)))(Subscript[S, 12])[ξ;0])       (3.4.10b) 
         ξ=Modulo[x ,Λ ]         (3.4.10c)   
then they are clearly periodic, and the solutions can be written:
   (({{(F^((+)))[x]}, {(F^((-)))[x]}}))=(({{(e^(i (q^((+))) x )), 0}, {0, (e^(i (q^((-))) x ))}})) . (({{(u^((+)))[x]}, {(u^((-)))[x]}}))         (3.4.11)
where (u^((±)))[x] are defined in (3.4.10b).

This is a proof of FB theorem for the case of Helmholtz equation using the transfer matrix approach. Of course it can be generalized to any order differential equation.

We will see what this means for a Bragg λ/4 stack and Bragg grating.

4. Transmission and reflection of a Bragg grating segment

In practice, all Bragg grating has finite length. (An infinitely long BG can actually be useless!). In other words, it is only a segment of some length L. We can find the reflectivity and transmission through such a segment. The results are:

4.1 Solving by substitution

To calculate the reflection coefficient:

r=-(((OverBar[G]) sinh(L Γ))/(i Γ cosh(L Γ)+Δβ sinh(L Γ)))

The above expression is convenient for real Γ, but for imaginary Γ, we can rewrite it differently.

We can write r for 2 region:
Region 1: when |Δβ| ≤  |G((|)^2),
(Subscript[r, 1])=  -(((OverBar[G]) sinh(L Γ))/(i Γ cosh(L Γ)+Δβ sinh(L Γ))) ~ tanh(L Γ) as Δβ->0
Region 2: when |Δβ| ≥  |G((|)^2),
(Subscript[r, 2])= ((i (OverBar[G]) sin(L χ) Δβ (χ cos(L χ)/Δβ+i  sin(L χ)))/(((Δβ)^2)-G (cos^2)(L χ) (OverBar[G]))) ~ (1/(Δβ)) as Δβ->∞

4.2 Discussion of the reflectivity

In both above expressions, the key parameters are |G| and the product L Γ or Lχ, which is  L |(Sqrt[|G((|)^2) - ((Δβ)^2)])|.
|G|, which is the grating modulation depth determines the boundary of the 2  regions.
When Δβ is small, and within the range bounded by |G|, (Subscript[r, 1]) is ~ -(((OverBar[G]) sinh(L Γ))/(i |G| cosh(L Γ))) which is ~ ((i  (OverBar[G]))/( |G| ))tanh(L Γ). If L |(Sqrt[|G((|)^2) - ((Δβ)^2)])| is large, the reflectivity approaches 1, i. e. the light will be totally reflected.
In the other range, when Δβ is large and outside the range bounded by |G|, (Subscript[r, 2]) ->   ((i (OverBar[G]) sin(L χ) Δβ (χ cos(L χ)/Δβ+i  sin(L χ)))/(((Δβ)^2)-G (cos^2)(L χ) (OverBar[G]))),  which -> (( (OverBar[G]) (sin^(
))(L χ)( cos(L χ)+i  sin(L χ)))/(Δβ)), it is just a small oscillatory function with amplitude decreasing as 1/Δβ. 
So we can say that the grating has a spatial-frequency bandwidth  ~ |G|  = (((Subscript[β, o])^2)/(2β))∫(SuperStar[E])(x,y) E(x,y) δϵ(x,y,z) dxy.
It is clear now how important it is to have a large control over the grating modulation depth: this is how one can control a grating-based filter bandwidth.
With a given grating modulation depth, one can engineer the reflectivity amplitude with grating length. The reflectivity is proportional to the grating length L and the grating "depth" G. Intuitively, this makes sense. The shorter the grating or the weaker the grating, the less is the reflectivity, i. e. the grating does not have much of an effect on the propagating beam.
Obviously the most interesting situation is when the grating strongly couples the forward and backward going wave. To see this we can plot r for the following situation. Let

Dependence of r on length and δϵ

refl[ β_,K_,L_,deps_ ,β0_]:=-(I ((β0)^2) Conjugate(deps) sin(L (Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))])))/(2 β (I (β-(K/2)) sin(L (Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))]))-(Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))]) cos(L (Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))]))))

First, let's see what happens when we increase the grating length:

Manipulate[ Plot[  Abs[ refl[ β,2,L,0.0025 , 1/1.5]]^2,{β,0.99,1.01}, PlotRange->{0,1}
, PlotStyle->{Hue[L*0.001,1,1],Thickness[0.005]}  , ImageSize->{500.,200},Frame->True , GridLines->Automatic ]         
	,	{L,100,50000}]

Let's see what happens when we increase the grating modulation depth

Manipulate[ deps=I δϵ;
              Plot[  Abs[ refl[ β,2,3000,deps , 1/1.5]]^2,{β,0.99,1.01}, PlotRange->{0,1}, PlotStyle->{Hue[δϵ*10^3,1,1],Thickness[0.005]} ],
{ δϵ,0.0001,0.05}]

It should be clear now why the ability to control index modulation depth and grating length are two most important aspects of Bragg grating fabrication technology

4.2 Application:

4.2.1 Plotting example of a notch filter

We can make a notch filter by making a cavity between two highly reflecting mirror. The mirror in this case is a broadband BG. To plot an example of a notch filter, we first define a function for the transmission:

trans [ β_,K_,L_,deps_ ,β0_]:=(I (E^(-I L (β-(K/2)))) (Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))]))/(I (Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))]) cos(L (Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))]))+(β-(K/2)) sin(L (Sqrt[(((β-(K/2)))^2)-((((β0)^4) ((\[LeftBracketingBar]deps\[RightBracketingBar])^2))/(4 (β^2)))])))

Now we can define a cavity transmission function:

CavityTrans[β_, K_, L_,lcav_,deps_, β0_   ]:=Module[{t, r},
		                       t= trans[β, K, L,deps, β0]; r= refl[β, K, L,deps, β0];
	                                                                                 Abs[(((E^(I β (lcav+L))) t^2)/(1-(E^(2 I lcav  β)) Abs[r]^2))]^2 ]

Manipulate[lcav=Pi*(1+δlcav ) ;
  Plot[ 10* Log[10, CavityTrans[ β,2,2500,lcav,0.01 , 1/1.5]] ,{β,0.995,1.005}, PlotRange->{-40,0}, PlotStyle->{Hue[δlcav*10^4,1,0.5],Thickness[0.005]} , Frame->True, GridLines->Automatic
		                       ]
,{δlcav,-0.00025,0.00025}]

Here, we don't even care exactly what's inside a BG section. We treat it as a black box with r and t coefficient. We then can build any thing just like we can do with optical elements.

4.2.2 Example of multiple segments, apodization

Here, we will look at the case with multiple segments BG or modulated BG.

Smatrix[Z1_,Z2_,β_,K_,deps_,β0_]:= (({{(((E^(L (Γ-I Δβ))) (Γ+(E^(2 I L Δβ)) (Γ-I Δβ)+I Δβ))/(2 Γ)), ((I (E^(2 z Γ+L (Γ-I Δβ))) (-1+(E^(2 I L Δβ))) G)/(2 Γ))}, {-((I (E^(-2 z Γ-L (Γ+I Δβ))) (-1+(E^(2 I L Δβ))) (OverBar[G]))/(2 Γ)), (((E^(-L (Γ+I Δβ))) (Γ+(E^(2 I L Δβ)) (Γ+I Δβ)-I Δβ))/(2 Γ))}}))

SingleM[z_,beta_, K_,deps_, beta0_ ]:=Module[{delbeta, gamma,G (*,m11  *)},
		                  G= beta0^2*deps/(2*beta) ;  
		                  delbeta= (beta -K/2) ; 
		           gamma=I*Sqrt[delbeta^2- Abs[G]^2] ;
		           (*  m11=Exp[-z*gamma-I*delbeta*z]; *)
		    (({{Exp[-z*gamma-I*delbeta*z], ((G*Exp[z*gamma-I*delbeta*z])/((delbeta+I*gamma)))}, {((Conjugate[G]*Exp[-z*gamma+I*delbeta*z])/((delbeta+I*gamma))), Exp[z*gamma+I*delbeta*z]}}))         ] ;
InverseM[z_,beta_, K_,deps_, beta0_ ]:=Module[{delbeta, gamma,G (*,m11  *)},
		                  G= beta0^2*deps/(2*beta );  
		                  delbeta= (beta -K/2) ; 
		           gamma=I*Sqrt[delbeta^2- Abs[G]^2] ;
		           (*  m11=Exp[-z*gamma-I*delbeta*z]; *)
		    (({{Exp[z*(gamma+I*delbeta)]*
	(gamma-I*delbeta), I*G*Exp[z*(gamma-I*delbeta)]}, {I*Conjugate[G]*
	Exp[z*(-gamma+I*delbeta)], Exp[-z*(gamma+I*delbeta)]*
	(gamma-I*delbeta)}}))/(2*gamma)         ] ;

Smatrix[Z1_,Z2_,β_,K_,deps_,β0_]:=             
            	SingleM[Z2,β,K,deps,β0].   Inverse[SingleM[Z1,β,K,deps,β0]];

MultiS[Z_List,β_,K_List,deps_List,β0_]:=Module[{n,Q},
		 n=Min[{Length[Z]-1,Length[K],Length[deps]}  ];Q=(({{1, 0}, {0, 1}}));
		          Do[ Q =
				            SingleM[Z[[i+1]],β,K[[i]],deps[[i]],β0].    InverseM[Z[[i]],β,K[[i]],deps[[i]],β0]  
				 (*  Smatrix[ Z[[i]],Z[[i+1]],β,K[[i]],deps[[i]],β0] *) . Q  ,{i,1,n}];Q] ;
refl[S_]:=Abs[S[[1]]/S[[2]]]^2 ;

Here is an example of Gaussian apodized BG

AllPlot={};depslist=0.008` Table[N[Exp[-((i^2)/(70.`))]],{i,-14,14}];
Zlist=Range[30] 200.`;
aBG=Plot[((Abs[r])^2)/.NSolve[{{t},{0}}==MultiS[Zlist,β,Table[2,{i,1,29}],depslist,(1/(1.5`))].{{1},{r}},{r,t}][[1,1]],{β,0.99,1.01},PlotRange->{0,1},PlotStyle->{Hue[0,1,1],Thickness[0.005]},ImageSize->{500.,200}]
ListPlot[depslist,Joined->True]

Compared with a uniform BG

nonaBG=Plot[  Abs[ refl[ β,2,6000,0.004 , 1/1.5]]^2,{β,0.99,1.01}, PlotRange->{0,1}
, PlotStyle->{Hue[0.7,1,1],Thickness[0.005]}
			                      , ImageSize->{500.,200}
		                       ]

Show[aBG,nonaBG]

HW illustration

Illustration of Floquet-Bloch solution on a 2 layer-Bragg Stack

Consider a Bragg stack

Write a module to calculate its S matrix
 S= ((((Subscript[k, H])+(Subscript[k, L])))/2)(({{(1/(Subscript[k, H])) (cos(Φ)+cos(δ) (Subscript[r, HL])), (1/((Subscript[k, L]) (Subscript[k, H])))(sin(Φ)-sin(δ) (Subscript[r, HL]))}, { -(sin(Φ)+sin(δ) (Subscript[r, HL])), (1/(Subscript[k, L])) (cos(Φ)-cos(δ) (Subscript[r, HL]))}}))

Now we plot the wave:
Let:                   (Subscript[ψ, q])[z]=({{{, ({{ (({{(e^(i (Subscript[k, L])z)), (e^(-i (Subscript[k, L])z))}})) . (({{(Subscript[c, q;1])}, {(Subscript[c, q;2])}})), z ∈ [0,(Subscript[l, 1])]  }, {(({{(e^(i (Subscript[k, H]) (z-(Subscript[l, 1])))), (e^(-i (Subscript[k, H])(z-(Subscript[l, 1]))))}})) . (({{(Subscript[d, q;1])}, {(Subscript[d, q;2])}})), z ∈ [(Subscript[l, 1]),(Subscript[l, 1])+(Subscript[l, 2])=Λ]  }})}})    (5.1.1)                   
We choose such that:
                        (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}})) =(({{1, 1}, {i (Subscript[k, L]), -i (Subscript[k, L])}})) . (({{(Subscript[c, q;1])}, {(Subscript[c, q;2])}}))    
Or:                (({{(Subscript[c, q;1])}, {(Subscript[c, q;2])}}))=(((({{1, 1}, {i (Subscript[k, L]), -i (Subscript[k, L])}})))^(-1)) . (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}}))=(1/2)(({{1, -(i/(Subscript[k, L]))}, {1, (i/(Subscript[k, L]))}})) . (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}}))     
and d such that it is continous with the previous segment, then, after one period:
                    (({{(Subscript[F, q])[Λ]}, {(Subscript[F, q])'[Λ]}}))=(e^(i β Λ)) . (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}}))      
So, we define a function                    
 
     (Subscript[u, q])[z]=(e^(-i β z))({{{, ({{ (({{(e^(i (Subscript[k, L])z)), (e^(-i (Subscript[k, L])z))}})) . (({{(Subscript[c, q;1])}, {(Subscript[c, q;2])}})), z ∈ [0,(Subscript[l, 1])]  }, {(({{(e^(i (Subscript[k, H])(z-(Subscript[l, 1])))), (e^(-i (Subscript[k, H])(z-(Subscript[l, 1]))))}})) . (({{(Subscript[d, q;1])}, {(Subscript[d, q;2])}})), z ∈ [(Subscript[l, 1]),(Subscript[l, 1])+(Subscript[l, 2])=Λ]  }})}})    (5.1.2)

                 (({{(Subscript[c, q;1])}, {(Subscript[c, q;2])}}))=(((({{1, 1}, {i (Subscript[k, L]), -i (Subscript[k, L])}})))^(-1)) . (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}}))=(1/2)(({{1, -(i/(Subscript[k, L]))}, {1, (i/(Subscript[k, L]))}})) . (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}}))     
                                    (({{(e^(i (Subscript[k, H])(Subscript[l, 2]))), (e^(-i (Subscript[k, H])(Subscript[l, 2])))}, {i (Subscript[k, H]) (e^(i (Subscript[k, H])(Subscript[l, 2]))), -i (Subscript[k, H])(e^(-i (Subscript[k, H])(Subscript[l, 2])))}})) (({{(Subscript[d, q;1])}, {(Subscript[d, q;2])}}))= (e^(i β Λ)) (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}}))

Hence                              (({{(Subscript[d, q;1])}, {(Subscript[d, q;2])}}))= (e^(i β Λ)) (({{(1/2) (e^(-i (Subscript[k, H]) (Subscript[l, 2]))), -((i (e^(-i (Subscript[k, H]) (Subscript[l, 2]))))/(2 (Subscript[k, H])))}, {(1/2) (e^(i (Subscript[k, H]) (Subscript[l, 2]))), ((i (e^(i (Subscript[k, H]) (Subscript[l, 2]))))/(2 (Subscript[k, H])))}})) (({{(Subscript[F, q])[0]}, {(Subscript[F, q])'[0]}}))

Also:                  (({{(e^(i (Subscript[k, L])(Subscript[l, 1]))), (e^(-i (Subscript[k, L])(Subscript[l, 1])))}, {i (Subscript[k, L]) (e^(i (Subscript[k, L])(Subscript[l, 1]))), -i (Subscript[k, L]) (e^(-i (Subscript[k, L])(Subscript[l, 1])))}})) . (({{(Subscript[c, q;1])}, {(Subscript[c, q;2])}}))= (({{1, 1}, {i (Subscript[k, H]), -i (Subscript[k, H])}})) .  (({{(Subscript[d, q;1])}, {(Subscript[d, q;2])}}))      
Or:                   Inverse[(({{1, 1}, {i (Subscript[k, H]), -i (Subscript[k, H])}}))] . (({{(e^(i (Subscript[k, L])(Subscript[l, 1]))), (e^(-i (Subscript[k, L])(Subscript[l, 1])))}, {i (Subscript[k, L]) (e^(i (Subscript[k, L])(Subscript[l, 1]))), -i (Subscript[k, L]) (e^(-i (Subscript[k, L])(Subscript[l, 1])))}})) . (({{(Subscript[c, q;1])}, {(Subscript[c, q;2])}}))= (({{(Subscript[d, q;1])}, {(Subscript[d, q;2])}}))

Inverse[(({{1, 1}, {I (Subscript[k, H]), -I (Subscript[k, H])}}))] . (({{(E^(I (Subscript[k, L])(Subscript[l, 1]))), (E^(-I (Subscript[k, L])(Subscript[l, 1])))}, {I (Subscript[k, L]) (E^(I (Subscript[k, L])(Subscript[l, 1]))), -I (Subscript[k, L]) (E^(-I (Subscript[k, L])(Subscript[l, 1])))}}))

SBG[BGParm_]:=Module[{ϵH,ϵL, l1,l2, q, λ,K0, sH, sL,kH,kL,Φ,δ,rHL},
        {{ϵL,ϵH}, {l1,l2},q,λ}=BGParm ;K0=2 π/λ ;
kH=K0 (Sqrt[ϵH-(q^2)]) ; kL=K0 (Sqrt[ϵL-(q^2)]) ;
Φ=kH l2+kL l1 ; δ=kH l2-kL l1  ; rHL=((kH-kL)/(kH+kL)) ;
(((kH+kL))/2)(({{(1/kH) (Cos[Φ]+Cos[δ]  rHL), (1/(kL kH))(Sin[Φ]-Sin[δ] rHL)}, { -(Sin[Φ]+Sin[δ] rHL), (1/kL) (Cos[Φ]-Cos[δ]  rHL)}}))
] ;
uFB[z_,β_,Ccoef_,Dcoef_,BGParm_]:=Module[{ζ,ϵH,ϵL, l1,l2, q, λ,K0, sH, sL,kH,kL,Φ,δ,rHL},
        {{ϵL,ϵH}, {l1,l2},q,λ}=BGParm ;K0=2 π/λ ;
kH=K0 (Sqrt[ϵH-(q^2)]) ; kL=K0 (Sqrt[ϵL-(q^2)]) ;
ζ=If[z==l1+l2,0,Mod[z,l1+l2] ];
{ (E^(-I β ζ)) If[ζ≤ l1, (E^(I kL ζ)) Ccoef[[1,1]]+ (E^(-I kL ζ)) Ccoef[[1,2]]
, (E^(I kH(ζ-l1))) Dcoef[[1,1]]+ (E^(-I kH(ζ-l1))) Dcoef[[1,2]]]
,(E^(I β ζ))  If[ζ≤ l1, (E^(I kL ζ)) Ccoef[[2,1]]+ (E^(-I kL ζ)) Ccoef[[2,2]]
, (E^(I kH(ζ-l1))) Dcoef[[2,1]]+ (E^(-I kH(ζ-l1))) Dcoef[[2,2]]]}
];
ψF[z_,β_,Ccoef_,Dcoef_,BGParm_]:=Module[{ζ,ϵH,ϵL, l1,l2, q, λ,K0, sH, sL,kH,kL,Φ,δ,rHL},
        {{ϵL,ϵH}, {l1,l2},q,λ}=BGParm ;K0=2 π/λ ;
kH=K0 (Sqrt[ϵH-(q^2)]) ; kL=K0 (Sqrt[ϵL-(q^2)]) ;
ζ=Mod[z,l1+l2] ;
{  If[ζ≤ l1, (E^(I kL ζ)) Ccoef[[1,1]]+ (E^(-I kL ζ)) Ccoef[[1,2]]
, (E^(I kH(ζ-l1))) Dcoef[[1,1]]+ (E^(-I kH(ζ-l1))) Dcoef[[1,2]]]
,  If[ζ≤ l1, (E^(I kL ζ)) Ccoef[[2,1]]+ (E^(-I kL ζ)) Ccoef[[2,2]]
, (E^(I kH(ζ-l1))) Dcoef[[2,1]]+ (E^(-I kH(ζ-l1))) Dcoef[[2,2]]]}
];

Manipulate[ λ=1;  K0=2 π/λ;{ϵL,ϵH}={(((nLo+I κLo))^2),(((nHi+I κHi))^2) } ;
 Δϵ=ϵH-ϵL ;   l=Λ*{η,1-η} ;{l1,l2}=l ; q=qR+I qI ;
kH=K0 (Sqrt[ϵH-(q^2)]) ; kL=K0 (Sqrt[ϵL-(q^2)]) ; 
BGParm={{ϵL,ϵH}, l,q,λ} ;
{eigenval,eigenvec}=Eigensystem[SBG[BGParm]] ;
eigenval2={Abs[eigenval],Arg[eigenval]}  ;
Ccoef={0.5(({{1, -(I/kL)}, {1, (I/kL)}})) . eigenvec[[1]]
,0.5(({{1, -(I/kL)}, {1, (I/kL)}})) .eigenvec[[2]] };
Dcoef={0.5*eigenval[[1]]  *(({{ (E^(-I kH l2)), -((I (E^(-I kH l2)))/( kH))}, { (E^(I kH l2)), ((I (E^(I kH l2)))/kH)}})) . eigenvec[[1]], 0.5*eigenval[[2]]*(({{ (E^(-I kH l2)), -((I (E^(-I kH l2)))/( kH))}, { (E^(I kH l2)), ((I (E^(I kH l2)))/kH)}})) . eigenvec[[2]]} ;  

If[Abs[Abs[eigenval[[1]] ]-1.]≤ (10^(-6)) , β=Arg[eigenval[[1]]]/Λ ,  β=-I Log[eigenval[[1]] ]/Λ ];
(* Grid[{{"S matrix=  ",SBG[BGParm]}
,{"   eigenvalues=   ",eigenval2,"   eigenvectors=   ",eigenvec}
 ,{" Ccoef 1=  ",  Ccoef[[1]]," Ccoef 2=  ", Ccoef[[2]]}
,{" Dcoef 1=  ",  Dcoef[[1]]," Dcoef 2=  ", Dcoef[[2]] }
,{Plot[Release[Re[ψF[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}]
,Plot[Release[Im[ψF[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}]}
,{Plot[Release[Re[uFB[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}]
,Plot[Release[Im[uFB[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}]}}]  *)
{CPPlot[eigenval,,,,]
,Manipulate[Show[BraggStackGr[Λ,η,Nperiod,{hue1,hue2} ],
Plot[{Re[(E^(I (β z-2 π t))) uFB[z,β,Ccoef,Dcoef,BGParm][[1]] ]+0.2
,Re[(E^(I (-β z-2 π t))) uFB[z,β,Ccoef,Dcoef,BGParm][[2]] ]-0.2}, {z,0,Nperiod Λ}, PlotStyle-> {Red,Purple},Filling->{1-> 0.2,2-> -0.2}, PlotPoints->50], AspectRatio-> 0.5
,PlotRange-> {-1,1}],
{t,0,1}]}
,Grid[{{Control[{{RUN,"OFF"},{"ON","OFF"}}],Control[{{rootgraph,"OFF"},{"ON","OFF"}}]},
{Control[{{nHi,1.6},1.5,3}],Control[{{κHi,0.},0,0.5}],Control[{{qR,0.},0,5}]},{Control[{{nLo,1.},1.,2}],Control[{{κLo,0},0,0.5}],Control[{{qI,0.},0,5}]},{Control[{{Λ,0.75,"period"},0.1,10}],Control[{{η,0.5},0.01,0.99}],Control[{{Nperiod,20},1,40,1}]},{Control[{{hue1,0.35},0,1}],Control[{{hue2,0.13},0,1}]}},Spacings->{5,1}]]

Manipulate[ λ=1;  K0=2 π/λ;{ϵL,ϵH}={(((nLo+I κLo))^2),(((nHi+I κHi))^2) } ;
 Δϵ=ϵH-ϵL ;   l=Λ*{η,1-η} ;{l1,l2}=l ; q=qR+I qI ;
kH=K0 (Sqrt[ϵH-(q^2)]) ; kL=K0 (Sqrt[ϵL-(q^2)]) ; 
BGParm={{ϵL,ϵH}, l,q,λ} ;
{eigenval,eigenvec}=Eigensystem[SBG[BGParm]] ;
eigenval2={Abs[eigenval],Arg[eigenval]}  ;
Ccoef={0.5(({{1, -(I/kL)}, {1, (I/kL)}})) . eigenvec[[1]]
,0.5(({{1, -(I/kL)}, {1, (I/kL)}})) .eigenvec[[2]] };
Dcoef={0.5*eigenval[[1]]  *(({{ (E^(-I kH l2)), -((I (E^(-I kH l2)))/( kH))}, { (E^(I kH l2)), ((I (E^(I kH l2)))/kH)}})) . eigenvec[[1]], 0.5*eigenval[[2]]*(({{ (E^(-I kH l2)), -((I (E^(-I kH l2)))/( kH))}, { (E^(I kH l2)), ((I (E^(I kH l2)))/kH)}})) . eigenvec[[2]]} ;  

If[Abs[Abs[eigenval[[1]] ]-1.]≤ (10^(-6)) , β=Arg[eigenval[[1]]]/Λ ,  β=-I Log[eigenval[[1]] ]/Λ ];
Grid[{{"S matrix=  ",SBG[BGParm]}
,{"   eigenvalues=   ",eigenval2,"   eigenvectors=   ",eigenvec}
 ,{" Ccoef 1=  ",  Ccoef[[1]]," Ccoef 2=  ", Ccoef[[2]]}
,{" Dcoef 1=  ",  Dcoef[[1]]," Dcoef 2=  ", Dcoef[[2]] }
,{Plot[Release[Re[ψF[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}, PlotRange-> All]
,Plot[Release[Im[ψF[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}, PlotRange-> All]}
,{Show[Plot[Release[Re[uFB[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}, PlotRange-> All],BraggStackGr[Λ,η,3,],Plot[Release[Re[uFB[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}, PlotRange-> All]]
,Show[Plot[Release[Im[uFB[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}, PlotRange-> All],BraggStackGr[Λ,η,3,],Plot[Release[Im[uFB[z,β,Ccoef,Dcoef,BGParm]]], {z,0,3Λ}, PlotStyle-> {Red,Blue}, PlotRange-> All]]}  }]
(* Manipulate[Show[BraggStackGr[Λ,η,20,{hue1,hue2} ],
Plot[{Re[(E^(I (β z-2 π t))) uFB[z,β,Ccoef,Dcoef,BGParm][[1]] ]+0.2
,Re[(E^(I (-β z-2 π t))) uFB[z,β,Ccoef,Dcoef,BGParm][[2]] ]-0.2}, {z,0,20Λ}, PlotStyle-> {Red,Purple},Filling->{1-> 0.2,2-> -0.2}], AspectRatio-> 0.5],
{t,0,1}] *)
,Grid[{{Control[{{RUN,"OFF"},{"ON","OFF"}}],Control[{{rootgraph,"OFF"},{"ON","OFF"}}]},
{Control[{{nHi,1.6},1.5,3}],Control[{{κHi,0.},0,0.5}],Control[{{qR,0.},0,5}]},{Control[{{nLo,1.5},1.,2}],Control[{{κLo,0},0,0.5}],Control[{{qI,0.},0,5}]},{Control[{{Λ,0.75,"period"},0.1,10}],Control[{{η,0.5},0.01,0.99}]},{Control[{{Reρmin,0,"Re s plot min"},-5,40}],Control[{{Reρmax,20,"Re s plot max"},0.,20}]},{Control[{{hue1,0.35},0,1}],Control[{{hue2,0.13},0,1}]}},Spacings->{5,1}]]