Plasmonics basic introduction- rebuilt

Wave equation (Helmholtz's equation)

Loading package

1. Introduction

1.1 Basic discussion

1.1.1 Wave is dielectric - evanescent

A plane wave solution of electromagnetic waves in dielectric
                                               E=(Subscript[E, 0])(e^(i ((Subscript[k, d]). r -ω t)))  
Where                                        (Subscript[k, d])=(Subscript[k, 0]) n    and    n=(Sqrt[ϵ μ]); assuming ϵ μ is real.
Consider the evanescent solution along the surface:
                                                E=(Subscript[E, 0])(e^(i (β z -ω t))) (e^(-q x ))    
                                                (β^2)-(q^2)=((Subscript[k, 0])^2) ϵ  (we take μ ~ 1)
 We define:   β =n((Subscript[k, 0])/(Cos[θ])), then:      q=(Sqrt[(β^2)-((Subscript[k, 0])^2) ϵ ])=(Sqrt[((Subscript[k, 0])^2) ϵ ((1/((Cos[θ])^2))-1) ])=(Subscript[k, d]) Tan[θ] =(Subscript[k, 0]) n Tan[θ]

ClearAll[λ, k0, β, q,n, ω, θ, s, l, Δz] ;
λ = 1; k0=((2 π)/λ); n= 1.4; ω =2 π ;
Manipulate[β =n(k0/(Cos[θ]))  ;
 q=(Sqrt[(β^2)-(n^2) (k0^2) ]) ; 
Animate[ Plot3D[Re[ (E^(I (β x - ω t)-q z))] ,{x,0,4},{z,0,2}
,BoxRatios->{5,2,2},PlotPoints->{40,15},ViewPoint->{-3,1,2}, Mesh->False, PlotRange->{-1,1}
	, ImageSize -> {300,200}] , {t,0,1}, AnimationRunning-> False]
, {θ, 10.*π/180, 30*Pi/180}]

Consider transverse magnetic solution:
                                                    H=(Subscript[E, 0])(OverHat[y]) (e^(i (β z -ω t))) (e^(-q x ))
                                                   (1/c)ϵ ((∂ (Overscript[E, ->]))/(∂ t))  = ∇×(Overscript[H, ->])  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x), (∂_y), (∂_z)}, {0, (Subscript[H, y]), 0}}))  
                                            -(1/c)i ω ϵ E  = ∇×H  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x), (∂_y), (∂_z)}, {0, (Subscript[H, y]), 0}}))=(Subscript[E, 0]) (- (OverHat[x]) i β - q (OverHat[z]) ) (e^(i (β z -ω t))) (e^(-q x ))  
                                          E  = (1/ϵ)(Subscript[E, 0]) ((OverHat[x]) (β/( (Subscript[k, 0]))) -i (q/(Subscript[k, 0])) (OverHat[z]) ) (e^(i (β z -ω t))) (e^(-q x ))        
                                          E  = (1/n)(Subscript[E, 0]) ((OverHat[x]) (β/(n (Subscript[k, 0]))) -i (q/(n (Subscript[k, 0]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^(-q x ))

ClearAll[λ, k0, β, q,n, ω, θ, s, l, Δz] ;
λ = 1; k0=((2 π)/λ); n= 1.4; ω =2 π ;θ=20*Pi/180.  ;
β =n(k0/(Cos[θ]))  ;   q=(Sqrt[(β^2)-(n^2) (k0^2) ]) ; 
Δz=2 π/(β*20.)  ;  s=0.5  ;
Manipulate[
Graphics3D[{Table[{Hue[0,1,1],
Arrow[{{ Δz l,0,0},{Δz l+s Tan[θ] Sin[2 π (Δz l- t)] /n ,0,s* Cos[2 π (Δz l- t)]/(n Cos[θ] )}} ]},{l,0,80}] , 
Table[{Hue[0.7,1,1],
Arrow[{{ Δz l,0,0},{Δz l  ,s Cos[2 π (Δz l- t)],0}} ]},{l,0,80}]},
ViewPoint-> {5,-8,3}, PlotRange-> {{0,2},{-1,1},{-0.5,0.5}}],
{t,0,1,0.025}, AnimationRunning-> False, DisplayAllSteps-> True]

ClearAll[λ, k0, β, q,n, ω, θ, s, l, Δz] ;
λ = 1; k0=((2 π)/λ); n= 1.4; ω =2 π ;θ=20*Pi/180.  ;
β =n(k0/(Cos[θ]))  ;   q=(Sqrt[(β^2)-(n^2) (k0^2) ]) ; 
Manipulate[
VectorPlot[{Tan[θ] Re[- I (E^(I (β x - ω  t)-q z))]/n,Re[ (E^(I (β x - ω  t)-q z))]/(n Cos[θ])}, {x, 0,1}, {z, 0,0.25}, AspectRatio-> 0.5, VectorPoints->Flatten[Table[{i,j},
{i,0,1,0.025},{j,0,0.25,0.025}],1] ]

,{t, 0,1}]

Field lines:                               E  =(1/n) (Subscript[E, 0]) ((OverHat[x])  (β/(n (Subscript[k, 0]))) -i (q/(n (Subscript[k, 0]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^(-q x ))  
                                   dx/du=β Cos[(β z -ω t)] (e^(-q x ))  dz/du=q Sin[(β z -ω t)] (e^(-q x ))
                                   {x[u], y[u]}
                                   ((d x)/(β Cos[(β z -ω t)] (e^(-q x ))))=((d z)/(q Sin[(β z -ω t)] (e^(-q x ))))
                             
                                   q Sin[(β z -ω t)] d x=β Cos[(β z -ω t)] d z
Calculate contour:       q x Sin[(β z -ω t)]   +C=Sin[(β z -ω t)]

This is not correct because of the alias integration

β= .;  ω= .;  q= .;
Integrate[(β/q) Cot[β z -ω t],z]

λ = 1; k0=((2 π)/λ); n= 1.4; ω =2 π ;θ=20*Pi/180.  ;
β =n(k0/(Cos[θ]))  ;   q=(Sqrt[(β^2)-(n^2) (k0^2) ]) ; 
Manipulate[Plot[Table[i+((log(Abs[Sin[z β-t ω]]))/q),{i,0,1,0.2}],{z,0,1 }, PlotRange-> {0,1}, AspectRatio->1],{t,0,1}]

1.1.2 Wave in metal

A solution of electromagnetic waves in metal is:
                                               E=(Subscript[E, 0])(e^(i ((OverHat[k]). r -ω t)))  
Where                                       (OverHat[k])=(Subscript[k, 0]) (OverHat[n])    and    
                                     ((OverHat[n])^2)=(n^2)  (((1+i κ))^2) =μ (OverHat[ϵ]) =μ (1+4 (Subscript[πχ, 0]) -((((Subscript[ω, p])^2) )/(ω (ω +i β ))))    
                                                                    =1+4 (Subscript[πχ, 0]) -((((Subscript[ω, p])^2) )/(ω (ω +i β )))  
since we approximate μ=1. 
Wave in metal is highly attenuated because of strong damping. Strictly speaking, there are both real and imaginary components that make the real wave propagation complicated, but if one can approximate:
                                                            (OverHat[ϵ])  ~1 -((((Subscript[ω, p])^2) )/(ω^2)) 
Then one can treat (OverHat[ϵ]) as real negative. In general, 
                                                            (OverHat[k])=(Subscript[k, 0]) (OverHat[n]) =(Subscript[k, 0])(Sqrt[(OverHat[ϵ])])=i (Subscript[k, 0]) (Sqrt[|(OverHat[ϵ])|]) 
means that there is propagating solution. However, if we write:
                                                          ((OverHat[k])^2)=((Subscript[k, 1])^2)-((Subscript[k, 2])^2)=((Subscript[k, 0])^2) (OverHat[ϵ])      
such as from an evanescent wave:          (e^(i (Subscript[k, 1]) . x )) (e^( -(Subscript[k, 2]) y )) 
then we can see the possibility:         ((Subscript[k, 1])^2)=((Subscript[k, 2])^2)+((Subscript[k, 0])^2) (OverHat[ϵ]) >0     if    ((Subscript[k, 2])^2)>-((Subscript[k, 0])^2) (OverHat[ϵ])   
Or                                                        ((Subscript[k, 2])^2)>((Subscript[k, 0])^2)  (((((Subscript[ω, p])^2) )/(ω^2))-1)   
This can happen near the plasmon frequency so that (((((Subscript[ω, p])^2) )/(ω^2))-1) is not too large.

Consider a TM wave at a metal surface, such as when there is a wave impinging on it. 
                                                         H=(OverHat[y]) (Subscript[E, 0]) (e^(i (β. z -ω t))) (e^( (Subscript[q, m]) x )) (x <0 )
The condition is:                                (β^2)-((Subscript[q, m])^2)=((Subscript[k, 0])^2) (OverHat[ϵ])

We can solve for E field
                                                          ∇×H-(ϵ/c)((∂ E)/(∂ t)) =((4 π)/c) σ E   
                                                          ∇×H+(ϵ/c)i ω E=((4 π)/c) σ E   
                                                          ∇×H=-(ϵ/c)i ω E+((4 π)/c) σ E   
                                                          ∇×H=-(ω/c)i  (ϵ +i((4 π σ)/( ω))  )E  
                                                          ∇×H=-(Subscript[k, 0]) i  (OverHat[ϵ]) E  
                                      -(Subscript[k, 0]) i  (OverHat[ϵ]) E=(({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x), (∂_y), (∂_z)}, {0, (Subscript[H, y]), 0}}))=(Subscript[E, 0]) (- (OverHat[x]) i β + (Subscript[q, m]) (OverHat[z]) ) (e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))
                                     E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (OverHat[ϵ]))) +i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))

ClearAll[λ, k0, β, q,n, ω,ωp,  θ, s, l, Δz] ;
λ = 1; k0=((2 π)/λ); ωp = 2 π  ; 
Manipulate[ ϵm=1-((ωp)^2)/(ω^2) ;
β=(Sqrt[(q^2)+ϵm (k0^2) ]) ; z=0. ;
Plot[{q  Re[I (E^(I (β x - ω  t)+q z))]/(k0 ϵm),β Re[ (E^(I (β x - ω  t)+q z))]/(k0 ϵm) }, {x, 0,1},  AspectRatio-> 0.5 ]

,{t, 0,1}, {q, -ϵm (k0^2), 100},{{ω,0.9 ωp},0.3 ωp, 0.99  ωp}]

ClearAll[λ, k0, β, q,n, ω,ωp,  θ, s, l, Δz] ;
λ = 1; k0=((2 π)/λ); ωp = 2 π  ; 
Manipulate[ ϵm=1-((ωp)^2)/(ω^2) ;
β=(Sqrt[(q^2)+ϵm (k0^2) ]) ;
VectorPlot[{q  Re[I (E^(I (β x - ω  t)+q z))]/(k0 ϵm),β Re[ (E^(I (β x - ω  t)+q z))]/(k0 ϵm) }, {x, 0,1}, {z,-0.25,0}, AspectRatio-> 0.5, VectorPoints->Flatten[Table[{i,j},
{i,0,1,0.025},{j,0,-0.25,-0.025}],1] ]
,{t, 0,1}, {q, -ϵm (k0^2), 100},{{ω,0.9 ωp},0.3 ωp, 0.99  ωp}]

1.1.3 Plasmon: Matching boundary condition

Comparing two solutions of E filed
                                    ( (OverHat[x])  (β/((Subscript[k, 0]) (OverHat[ϵ]))) +i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ]))) (OverHat[z]) )    and  ( (OverHat[x])  (β/((Subscript[k, 0]) (Subscript[ϵ, d]))) -i ((Subscript[q, d])/((Subscript[k, 0]) (Subscript[ϵ, d]))) (OverHat[z]) )
Tangential continuity requirement: 
                                                          ((Subscript[q, d])/(Subscript[ϵ, d]))=-((Subscript[q, m])/(OverHat[ϵ]))
This is possible because                     (OverHat[ϵ])<0
Thus, as a function of                   (1/(Subscript[ϵ, d]))(Sqrt[(β^2)-((Subscript[k, 0])^2) (Subscript[ϵ, d])])=-(1/(OverHat[ϵ])) (Sqrt[(β^2)-((Subscript[k, 0])^2) (OverHat[ϵ])])  
                                                    ((β^2)-((Subscript[k, 0])^2) (Subscript[ϵ, d]))((OverHat[ϵ])^2)= ((Subscript[ϵ, d])^2)  ((β^2)-((Subscript[k, 0])^2) (OverHat[ϵ]))

1.2 Plasmon dispersion

1.2.1 Basic formula

β=. ; (Subscript[k, 0])= .;
FullSimplify[Solve[ ((β^2)-((Subscript[k, 0])^2) (Subscript[ϵ, d])) ((OverHat[ϵ])^2)==  ((Subscript[ϵ, d])^2)  ((β^2)-((Subscript[k, 0])^2) (OverHat[ϵ])) ,β]]

FullSimplify[(((Subscript[k, 0]^2) (OverHat[ϵ]) ((OverHat[ϵ])-(Subscript[ϵ, d])) (Subscript[ϵ, d]))/(((OverHat[ϵ])^2)-(Subscript[ϵ, d]^2)))]

Thus, we have the dispersion relation:
                                                            (β^2)=(Subscript[k, 0]^2)(( (OverHat[ϵ]) (Subscript[ϵ, d]))/((Subscript[ϵ, d])+(OverHat[ϵ])))
Hence:                                  ((Subscript[q, d])^2)=(β^2)-((Subscript[k, 0])^2)  (Subscript[ϵ, d])=((Subscript[k, 0])^2) ((((OverHat[ϵ]) (Subscript[ϵ, d]))/((Subscript[ϵ, d])+(OverHat[ϵ])))- (Subscript[ϵ, d])) =-((Subscript[k, 0])^2) (((Subscript[ϵ, d])^2)/((Subscript[ϵ, d])+(OverHat[ϵ])))   
                                              ((Subscript[q, m])^2)=(β^2)-((Subscript[k, 0])^2)  (OverHat[ϵ])=-((Subscript[k, 0])^2) (((OverHat[ϵ])^2)/((Subscript[ϵ, d])+(OverHat[ϵ])))   
                           (OverHat[ϵ])  (Subscript[q, d])=(Subscript[k, 0]) (OverHat[ϵ])(Sqrt[(((Subscript[ϵ, d])^2)/(-((Subscript[ϵ, d])+(OverHat[ϵ]))))]) =(Subscript[k, 0])(( (OverHat[ϵ])  (Subscript[ϵ, d]))/(Sqrt[-((Subscript[ϵ, d])+(OverHat[ϵ]))]))  ;   (Subscript[ϵ, d]) (Subscript[q, m])=(Subscript[k, 0]) (Subscript[ϵ, d])(Sqrt[- (((OverHat[ϵ])^2)/((Subscript[ϵ, d])+(OverHat[ϵ])))]) =(Subscript[k, 0])(( (OverHat[ϵ])  (Subscript[ϵ, d]))/(Sqrt[-((Subscript[ϵ, d])+(OverHat[ϵ]))]))

EpsMetal[λ_,Parm_List]:=Module[{ω, wp, β, χ},
{wp, β, χ}=Parm ; ω=1.23984187542/λ  ;
1+4 Pi χ -(((wp^2) )/(ω  ( ω +I β )))  ]  ;
B2[ϵd_,ϵm_]:=(( ϵm ϵd)/(ϵm+ ϵd));
Qd[ϵd_,ϵm_]:=Sqrt[-((  ((ϵd)^2))/(ϵm+ ϵd))];
Qm[ϵd_,ϵm_]:= Sqrt[-((  ((ϵm)^2))/(ϵm+ ϵd))] ;

wp=6. ; 
Manipulate[
GraphicsRow[{
Plot[Re[EpsMetal[λ,{wp, γ, χ}] ],
{λ, 0.2,0.5}, Frame-> True, GridLines-> Automatic]
,Plot[Im[EpsMetal[λ,{wp, γ, χ}] ],
{λ, 0.2,0.5}, Frame-> True, GridLines-> Automatic] }]
, {γ,0.01,0.1},{χ,0.,0.1}]

plasmonparm={wp, γ, χ} ;
wp=6. ;γ=0.02; χ=0.001  ;
GraphicsRow[{Plot[Re[Sqrt[B2[1.,EpsMetal[λ,plasmonparm]]]] ,{λ,0.4, 1}, PlotRange-> All],
Plot[Im[Sqrt[B2[1.,EpsMetal[λ,plasmonparm] ]]],{λ,0.4, 1}, PlotRange-> All]}]
GraphicsRow[{Plot[Re[Sqrt[B2[1,EpsMetal[1.23984187542/ω ,plasmonparm]]]] ,{ω,2, 4}, PlotRange-> All],
Plot[Im[Sqrt[B2[1.,EpsMetal[1.23984187542/ω ,plasmonparm]  ]]],{ω,2, 4}, PlotRange-> All]}]

We notice that plasmon has quite large dispersion

plasmonparm={wp, γ, χ} ;
wp=6. ; γ=0.02; χ=0.001  ;
GraphicsGrid[{{Plot[Re[Qm[1.4^2,EpsMetal[λ,plasmonparm]]] ,{λ,0.25, 0.5}, PlotRange-> All]
,Plot[Im[Qm[1.4^2,EpsMetal[λ,plasmonparm]]],{λ,0.25, 0.5}, PlotRange-> All]}
, {Plot[Re[Qd[1.4^2,EpsMetal[λ,plasmonparm]]] ,{λ,0.25, 0.5}, PlotRange-> All]
,Plot[Im[Qd[1.4^2,EpsMetal[λ,plasmonparm] ]],{λ,0.25, 0.5}, PlotRange-> All]}}]

1.2.2 Group velocity

ClearAll[wp, γ, χ, ω, c, (OverHat[ϵ])] ;
B2[ϵd_,ϵm_]:=(( ϵm ϵd)/(ϵm+ ϵd));
FullSimplify[(∂_ω)((Sqrt[B2[(Subscript[ϵ, d])[ω],(OverHat[ϵ])[ω]]]))]

(((((OverHat[ϵ])[ω])^2) (((Subscript[ϵ, d])))'(ω)+(((Subscript[ϵ, d])(ω))^2) (OverHat[ϵ])'(ω))/(2 ((Subscript[ϵ, d])[ω]  (OverHat[ϵ])[ω])^2  FullSimplify[(Sqrt[(((Subscript[ϵ, d])(ω) (OverHat[ϵ])[ω])/((OverHat[ϵ])[ω]+(Subscript[ϵ, d])(ω)))]) (((((OverHat[ϵ])[ω]+(Subscript[ϵ, d])(ω)))^2)/(((Subscript[ϵ, d])[ω]  (OverHat[ϵ])[ω])^2))]))

((((((Subscript[ϵ, d])(ω) (OverHat[ϵ])[ω])/((OverHat[ϵ])[ω]+(Subscript[ϵ, d])(ω)))))^(3/2))FullSimplify[((((((OverHat[ϵ])[ω])^2) (((Subscript[ϵ, d])))'(ω)+(((Subscript[ϵ, d])(ω))^2) (OverHat[ϵ])'(ω)))/(2 (((Subscript[ϵ, d])(ω))^2) (((OverHat[ϵ])[ω])^2)))]

Hence:  (∂_ω)β[ω]=(Subscript[k, 0]) (1/2) ((((((Subscript[ϵ, d])(ω) (OverHat[ϵ])[ω])/((OverHat[ϵ])[ω]+(Subscript[ϵ, d])(ω)))))^(3/2)) ((((OverHat[ϵ])'(ω))/(((OverHat[ϵ])[ω])^2))+(((((Subscript[ϵ, d])))'(ω))/(((Subscript[ϵ, d])(ω))^2)))  
              (1/(Subscript[v, g]))= (∂_ω)β[ω]=(Subscript[k, 0]) (1/2) (((((β[ω])/(Subscript[k, 0]))))^(3/2)) ((((OverHat[ϵ])'(ω))/(((OverHat[ϵ])[ω])^2))+(((((Subscript[ϵ, d])))'(ω))/(((Subscript[ϵ, d])(ω))^2)))

(OverHat[ϵ])'(ω)= (∂_ω)(1+4 π χ(ω)-((wp^2)/(ω (ω+I γ))))

Simplify[((wp^2)/((ω^2) (ω+I γ)))+((wp^2)/(ω (((ω+I γ))^2)))]+4 π χ'(ω)

Hence:
            (1/(Subscript[v, g]))= (∂_ω)β[ω]=(Subscript[k, 0]) (1/2) (((((β[ω])/(Subscript[k, 0]))))^(3/2)) ((((OverHat[ϵ])'(ω))/(((OverHat[ϵ])[ω])^2))+(((((Subscript[ϵ, d])))'(ω))/(((Subscript[ϵ, d])(ω))^2)))  
            (1/(Subscript[v, g]))= (∂_ω)β[ω]=(Subscript[k, 0]) (1/2) (((((β[ω])/(Subscript[k, 0]))))^(3/2)) (((4 π χ'(ω)+(((wp^2) (2 ω+i γ))/((ω^2) (((ω+i γ))^2))))/(((OverHat[ϵ])[ω])^2))+(((((Subscript[ϵ, d])))'(ω))/(((Subscript[ϵ, d])(ω))^2)))  
            (c/(Subscript[v, g]))=c (∂_ω)β[ω]=ω (1/2) (((((β[ω])/(Subscript[k, 0]))))^(3/2)) (((4 π χ'(ω)+(((wp^2) (2 ω+i γ))/((ω^2) (((ω+i γ))^2))))/(((OverHat[ϵ])[ω])^2))+(((((Subscript[ϵ, d])))'(ω))/(((Subscript[ϵ, d])(ω))^2)))  
            (c/(Subscript[v, g]))=c (∂_ω)β[ω]= (1/2) (((((β[ω])/(Subscript[k, 0]))))^(3/2)) (((4 π χ'(ω) ω+(((wp^2) (2 ω+i γ))/(ω (((ω+i γ))^2))))/(((OverHat[ϵ])[ω])^2))+((ω (((Subscript[ϵ, d])))'(ω))/(((Subscript[ϵ, d])(ω))^2)))  
 In the limit ω-> 0,
             (c/(Subscript[v, g]))=c (∂_ω)β[ω]= (1/2) (((((β[ω])/(Subscript[k, 0]))))^(3/2)) (((4 π χ'(ω) ω+(((wp^2) (2 ω+i γ))/(ω (((ω+i γ))^2))))/(((OverHat[ϵ])[ω])^2))+((ω (((Subscript[ϵ, d])))'(ω))/(((Subscript[ϵ, d])(ω))^2)))  
 Example:

Limit[((((wp^2) (2 ω+I γ))/(ω (((ω+I γ))^2)))/(((1+4 π (Subscript[χ, 0])-((wp^2)/(ω (ω+I γ)))))^2)),ω-> 0]

Limit[ (Sqrt[B2[(Subscript[ϵ, d])[0],(1+4 π (Subscript[χ, 0])-((wp^2)/(ω (ω+I γ))))]]) , ω-> 0]

We see that the limit of (Subscript[n, g]) is 0 for ω.  The group velocity travels faster than the speed of light?

EpsMetal[λ_,Parm_List]:=Module[{ω, wp, β, χ},
{wp, β, χ}=Parm ; ω=1.23984187542/λ  ;
1+4 Pi χ -(((wp^2) )/(ω  ( ω +I β )))  ]  ;

ϵd= 1.4^2 ;
wp=6. ; γ=0.02; χ=0.001  ;
W=1.23984187542  ;
Plot[(1/2)Re[  (((B2[ϵd,EpsMetal[λ,{wp, γ, χ}]]))^(3/2)) (((((wp^2) (2 W/λ+I γ))/((W/λ) (((W/λ+I γ))^2)))/((EpsMetal[λ,{wp, γ, χ}])^2)))],{λ,1,0.34}, PlotRange-> {-10,10}, Frame-> True, GridLines-> Automatic]

We notice that group velocity drops below 1. WHy?

ϵd= 1.4^2 ;
wp=6. ; γ=0.02; χ=0.001  ;
W=1.23984187542  ;
Plot[Re[ (1/2) (((B2[ϵd,EpsMetal[λ,{wp, γ, χ}]]))^(3/2)) (((((wp^2) (2 W/λ+I γ))/((W/λ) (((W/λ+I γ))^2)))/((EpsMetal[λ,{wp, γ, χ}])^2)))],{λ,1,0.34}, PlotRange-> {0,5.}, Frame-> True, GridLines-> Automatic]

Notice

1.2.3 GVD

This section cannot be recovered.

1.3 Plasmon plot

Amplitude plot

EpsMetal[λ_,Parm_List]:=Module[{ω, wp, β, χ},
{wp, β, χ}=Parm ; ω=1.23984187542/λ  ;
1+4 Pi χ -(((wp^2) )/(ω  ( ω +I β )))  ]  ;
B2[ϵd_,ϵm_]:=(( ϵm ϵd)/(ϵm+ ϵd));
Qd[ϵd_,ϵm_]:=Sqrt[-((  ((ϵd)^2))/(ϵm+ ϵd))];
Qm[ϵd_,ϵm_]:= Sqrt[-((  ((ϵm)^2))/(ϵm+ ϵd))] ;

E field in metal:  E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (OverHat[ϵ]))) +i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))  
E field in dielectric:    E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (Subscript[ϵ, d]))) -i ((Subscript[q, d])/((Subscript[k, 0]) (Subscript[ϵ, d]))) (OverHat[z]) )

λ = 0.532; k0=((2 π)/λ); ϵd= 1.48^2 ; ω =2 π ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ 
plasmonparm={wp, γ, χ} ;
β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm]]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm] ]  ;
Plot3D[If[y≤0., Re[I  qm  (E^(I  ( β  x -ω  t))) (E^( qm  y ))/(k0  EpsMetal[λ,plasmonparm ] )]  
,Re[ - I qd  (E^(I  ( β  x - ω  t ))) (E^(-qd  y ))/(k0  ϵd )]  ]
,{x,0,1},{y,-0.2,0.4 }, PlotRange-> All, BoxRatios-> {10,4,5}, PlotPoints-> {50,20}, PlotLabel->Style["E z (parallel) of surface plasmon, showing continuity",Blue,15]] 
,{t,0,1},{λ,0.532,1.24/(0.6*wp)},{{γ,0.02},0.,0.1}]

E field in metal:  E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (OverHat[ϵ]))) +i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))  
E field in dielectric:    E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (Subscript[ϵ, d]))) -i ((Subscript[q, d])/((Subscript[k, 0]) (Subscript[ϵ, d]))) (OverHat[z]) ) (e^(i (β z - ω t))) (e^(-(Subscript[q, d]) x ))

λ = 0.532; k0=((2 π)/λ); ϵd= 1.48^2 ; ω =2 π ;
plasmonparm={wp, γ, χ} ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm]]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm] ] ;
Plot3D[If[y≤0, Re[β (E^(I (β x -ω t))) (E^(qm y ))/(k0  EpsMetal[λ,plasmonparm ] )]
,Re[ β  (E^(I (β x -ω t))) (E^(-qd y ))/(k0  ϵd )] ]
,{x,0,1},{y,-0.2,0.4 }, PlotRange->All, BoxRatios-> {10,4,5}, PlotPoints-> {50,20}, PlotLabel->Style["E x (perp) of surface plasmon",Blue,24]],{t,0,1},{λ, 0.532,1.24/(0.8*wp)}]

Vector and field line plot

EpsMetal[λ_,Parm_List]:=Module[{ω, wp, β, χ},
{wp, β, χ}=Parm ; ω=1.23984187542/λ  ;
1+4 Pi χ -(((wp^2) )/(ω  ( ω +I β )))  ]  ;
B2[ϵd_,ϵm_]:=(( ϵm ϵd)/(ϵm+ ϵd));
Qd[ϵd_,ϵm_]:=Sqrt[-((  ((ϵd)^2))/(ϵm+ ϵd))];
Qm[ϵd_,ϵm_]:= Sqrt[-((  ((ϵm)^2))/(ϵm+ ϵd))] ;

plasmonparm={wp, γ, χ} ;
wp=6. ; γ=0.02; χ=0.001  ;
GraphicsGrid[{{Plot[Re[Qm[1.4^2,EpsMetal[λ,plasmonparm]]] ,{λ,0.25, 0.5}, PlotRange-> All]
,Plot[Im[Qm[1.4^2,EpsMetal[λ,plasmonparm]]],{λ,0.25, 0.5}, PlotRange-> All]}
, {Plot[Re[Qd[1.4^2,EpsMetal[λ,plasmonparm]]] ,{λ,0.25, 0.5}, PlotRange-> All]
,Plot[Im[Qd[1.4^2,EpsMetal[λ,plasmonparm]]],{λ,0.25, 0.5}, PlotRange-> All]}}]

E field in metal:  E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (OverHat[ϵ]))) +i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))  
E field in dielectric:    E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (Subscript[ϵ, d]))) -i ((Subscript[q, d])/((Subscript[k, 0]) (Subscript[ϵ, d]))) (OverHat[z]) )

λ = 0.532; k0=((2 π)/λ); ϵd= 1.48^2 ; ω =2 π ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ 
plasmonparm={wp, γ, χ} ;
β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm]]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm]]  ;
VectorPlot[If[y<0.,{ Re[I  qm  (E^(I  ( β  x -ω  t))) (E^( qm  y ))/(k0  EpsMetal[λ,plasmonparm ] )]  
,Re[ β (E^(I  ( β  x - ω  t ))) (E^(qm  y ))/(k0  EpsMetal[λ,plasmonparm ] )]  }
,{Re[ - I qd  (E^(I  ( β  x - ω  t ))) (E^(-qd  y ))/(k0  ϵd )]
,Re[ β (E^(I  ( β  x - ω  t ))) (E^(-qd  y ))/(k0  ϵd )]}  ]
,{x,0,0.3},{y,-0.1,0.1 }, AspectRatio-> 0.8, VectorPoints->Flatten[Table[{i,j},
{i,0,0.3,0.01},{j,-0.10625,0.1,0.0125}],1] 
,  PlotLabel->Style["E field of surface plasmon",Blue,24]] 
,{t,0,1},{λ,0.532,1.24/(0.6*wp)},{{γ,0.02},0.,0.1}]

Field line metal

E field in metal:  E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (OverHat[ϵ]))) +i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))  
E field in dielectric:    E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (Subscript[ϵ, d]))) -i ((Subscript[q, d])/((Subscript[k, 0]) (Subscript[ϵ, d]))) (OverHat[z]) ) (e^(i (β z - ω t))) (e^(-(Subscript[q, d]) x ))
To consider field line: 
In metal:           ((dx )/(Re[(β/((Subscript[k, 0]) (OverHat[ϵ])))(e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))]))=((dz )/(Re[i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ])))(e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))]))
                          ((dx )/(Re[(β/( (OverHat[ϵ])))(e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))(e^(- (Subscript[β, J]) z )) (e^((Subscript[q, m;R]) x ))]))=((dz )/(Re[i ((Subscript[q, m])/(OverHat[ϵ]))(e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))(e^(- (Subscript[β, J]) z )) (e^((Subscript[q, m;R]) x ))]))
                          ((dx )/(Re[(β/( (OverHat[ϵ])))(e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))]))=((dz )/(Re[i ((Subscript[q, m])/(OverHat[ϵ]))(e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))]))
                          ((dx )/(Re[((|β|)/(|(OverHat[ϵ])| ))(e^(i ((Subscript[ϕ, β])-(Subscript[ϕ, ϵ])))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))]))=((dz )/(Re[ ((|(Subscript[q, m])|)/(|(OverHat[ϵ])| ))(e^(i ((Subscript[ϕ, q])-(Subscript[ϕ, ϵ])+π/2))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))])) 
                          ((dx )/(|β| Re[(e^(i ((Subscript[ϕ, β])-(Subscript[ϕ, ϵ])))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))]))=((dz )/(|(Subscript[q, m])| Re[  (e^(i ((Subscript[ϕ, q])-(Subscript[ϕ, ϵ])+π/2))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(i (Subscript[q, m;I]) x ))])) 
                          ((dx )/(|β| Cos[((Subscript[β, R]) z -ω t) +(Subscript[q, m;J]) x +((Subscript[ϕ, β])-(Subscript[ϕ, ϵ]))]))=((dz )/(-|(Subscript[q, m])| Sin[((Subscript[β, R]) z -ω t) +(Subscript[q, m;J]) x +((Subscript[ϕ, q])-(Subscript[ϕ, ϵ]))])) 
Thus, we must solve:
                 |(Subscript[q, m])| Sin[((Subscript[β, R]) z -ω t) +(Subscript[q, m;J]) x +((Subscript[ϕ, q])-(Subscript[ϕ, ϵ]))] dx +|β| Cos[((Subscript[β, R]) z -ω t) +(Subscript[q, m;J]) x +((Subscript[ϕ, β])-(Subscript[ϕ, ϵ]))] dz=0. 
Change variable     (Subscript[β, R])ζ-> (Subscript[β, R]) z +(Subscript[q, m;J]) x-ω t -(Subscript[ϕ, ϵ]) 
                      dz= dζ-((Subscript[q, m;J])/(Subscript[β, R]))dx 
|(Subscript[q, m])| Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] dx +|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ -|β|((Subscript[q, m;J])/(Subscript[β, R]))  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])]dx=0. 
(|(Subscript[q, m])| Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] -|β|((Subscript[q, m;J])/(Subscript[β, R]))  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])])dx  +|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0. 
(|(Subscript[q, m])| Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] -((|(Subscript[q, m])| Sin[(Subscript[ϕ, q])])/(Cos[(Subscript[ϕ, β])]))  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])])dx  +|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0. 
((|(Subscript[q, m])|)/(Cos[(Subscript[ϕ, β])])) ( Cos[(Subscript[ϕ, β])] Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] -Sin[(Subscript[ϕ, q])]  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])])dx  +|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0.

FullSimplify[Cos[(Subscript[ϕ, b])] Sin[(Subscript[b, R]) ζ +(Subscript[ϕ, q])] -Sin[(Subscript[ϕ, q])]  Cos[(Subscript[b, R]) ζ +(Subscript[ϕ, b])]]

((|(Subscript[q, m])|)/(Cos[(Subscript[ϕ, β])]))  Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])] Sin[(Subscript[β, R]) ζ ] dx  +|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0.

DSolve[ G'[z]==- ((aβ  Cos[(Subscript[ϕ, b])] Cos[(Subscript[b, R]) z +(Subscript[ϕ, b])])/(aq  Cos[(Subscript[ϕ, b])-(Subscript[ϕ, q])]  Sin[(Subscript[b, R]) z ])),G[z],z]

Hence: the field lines are given by:
            x=-((|β|)/(|q|))  cos((Subscript[ϕ, β])) sec((Subscript[ϕ, β])-(Subscript[ϕ, q])) (((cos((Subscript[ϕ, β])) log(sin(z (Subscript[β, R]))))/(|β| cos((Subscript[ϕ, β]))))-z sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=-(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(sin(z (Subscript[β, R])))-z |β|cos((Subscript[ϕ, β])) sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=-(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(sin(z (Subscript[β, R])))-(Subscript[β, R]) z  sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=-(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(sin(z (Subscript[β, R]) + (Subscript[q, J]) x))-(z (Subscript[β, R]) + (Subscript[q, J]) x)  sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=-(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(sin(z (Subscript[β, R]) + (Subscript[q, J]) x))-(z (Subscript[β, R]) + (Subscript[q, J]) x)  sin((Subscript[ϕ, β]))) + (Subscript[c, 1])

λ = 1; k0=((2 π)/λ); ϵd= 1.4^2 ; ω =2 π ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ 
plasmonparm={wp, γ, χ} ;
β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm] ]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm] ]  ;
bR=Re[β]  ; ϕb=Arg[β]  ;
 qmI=Im[qm]  ;ϕqm=Arg[qm]  ; aqm=Abs[qm]  ;
 qdI=Im[qd]  ;ϕqd=Arg[qd]  ; aqd=Abs[qd]  ;
ContourPlot[z  aqm  Cos[ϕb-ϕqm]+  ( Cos[ϕb] Log[Abs[Sin[x bR+ qmI z-ω t]]]-(x bR+ qmI  z-ω t)  Sin[ϕb]),
{x,0,0.5},{z,-0.2,0}, AspectRatio-> 0.4, Contours-> {-0.3,-0.7,-1.5,-2.5,-4}, ContourShading-> None, AccuracyGoal-> 8, PlotPoints-> {80,10}] 
,{t,0,1},{λ,0.4,1.24/(0.6*wp)},{{γ,0.02},0.,0.1}]

Field line dielectric

E field in dielectric:    E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (Subscript[ϵ, d]))) -i ((Subscript[q, d])/((Subscript[k, 0]) (Subscript[ϵ, d]))) (OverHat[z]) ) (e^(i (β z - ω t))) (e^(-(Subscript[q, d]) x ))
To consider field line: 
In die:           ((dx )/(Re[(β/((Subscript[k, 0]) (OverHat[ϵ])))(e^(i (β z -ω t))) (e^(-(Subscript[q, d]) x ))]))=((dz )/(Re[-i ((Subscript[q, d])/((Subscript[k, 0]) (OverHat[ϵ])))(e^(i (β z -ω t))) (e^(-(Subscript[q, d]) x ))]))
                          ((dx )/(Re[(β/( (OverHat[ϵ])))(e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, d;I]) x ))(e^(- (Subscript[β, J]) z )) (e^(-(Subscript[q, d;R]) x ))]))=((dz )/(Re[-i ((Subscript[q, m])/(OverHat[ϵ]))(e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, m;I]) x ))(e^(- (Subscript[β, J]) z )) (e^(-(Subscript[q, d;R]) x ))]))
                          ((dx )/(Re[(β/( (OverHat[ϵ])))(e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, m;I]) x ))]))=((dz )/(Re[i ((Subscript[q, m])/(OverHat[ϵ]))(e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, m;I]) x ))]))
                          ((dx )/(Re[((|β|)/(|(OverHat[ϵ])| ))(e^(i ((Subscript[ϕ, β])-(Subscript[ϕ, ϵ])))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, m;I]) x ))]))=((dz )/(Re[ ((|(Subscript[q, d])|)/(|(OverHat[ϵ])| ))(e^(i ((Subscript[ϕ, q])-(Subscript[ϕ, ϵ])-π/2))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, m;I]) x ))])) 
                          ((dx )/(|β| Re[(e^(i ((Subscript[ϕ, β])-(Subscript[ϕ, ϵ])))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, m;I]) x ))]))=((dz )/(|(Subscript[q, d])| Re[  (e^(i ((Subscript[ϕ, q])-(Subscript[ϕ, ϵ])-π/2))) (e^(i ((Subscript[β, R]) z -ω t))) (e^(-i (Subscript[q, m;I]) x ))])) 
                          ((dx )/(|β| Cos[((Subscript[β, R]) z -ω t) -(Subscript[q, m;J]) x +((Subscript[ϕ, β])-(Subscript[ϕ, ϵ]))]))=((dz )/(|(Subscript[q, m])| Sin[((Subscript[β, R]) z -ω t) -(Subscript[q, m;J]) x +((Subscript[ϕ, q])-(Subscript[ϕ, ϵ]))])) 
Thus, we must solve:
                 |(Subscript[q, m])| Sin[((Subscript[β, R]) z -ω t) -(Subscript[q, d;J]) x +((Subscript[ϕ, q])-(Subscript[ϕ, ϵ]))] dx -|β| Cos[((Subscript[β, R]) z -ω t) -(Subscript[q, d;J]) x +((Subscript[ϕ, β])-(Subscript[ϕ, ϵ]))] dz=0. 
Change variable     (Subscript[β, R])ζ-> (Subscript[β, R]) z -(Subscript[q, d;J]) x-ω t -(Subscript[ϕ, ϵ]) 
                      dz= dζ+((Subscript[q, d;J])/(Subscript[β, R]))dx 
|(Subscript[q, m])| Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] dx -|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ -|β|((Subscript[q, m;J])/(Subscript[β, R]))  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])]dx=0. 
(|(Subscript[q, m])| Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] -|β|((Subscript[q, m;J])/(Subscript[β, R]))  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])])dx  -|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0. 
(|(Subscript[q, m])| Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] -((|(Subscript[q, m])| Sin[(Subscript[ϕ, q])])/(Cos[(Subscript[ϕ, β])]))  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])])dx  -|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0. 
((|(Subscript[q, m])|)/(Cos[(Subscript[ϕ, β])])) ( Cos[(Subscript[ϕ, β])] Sin[(Subscript[β, R]) ζ +(Subscript[ϕ, q])] -Sin[(Subscript[ϕ, q])]  Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])])dx  -|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0.

FullSimplify[Cos[(Subscript[ϕ, b])] Sin[(Subscript[b, R]) ζ +(Subscript[ϕ, q])] -Sin[(Subscript[ϕ, q])]  Cos[(Subscript[b, R]) ζ +(Subscript[ϕ, b])]]

((|(Subscript[q, m])|)/(Cos[(Subscript[ϕ, β])]))  Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])] Sin[(Subscript[β, R]) ζ ] dx  -|β| Cos[(Subscript[β, R]) ζ +(Subscript[ϕ, β])] dζ =0.

DSolve[ G'[z]== ((aβ  Cos[(Subscript[ϕ, b])] Cos[(Subscript[b, R]) z +(Subscript[ϕ, b])])/(aq  Cos[(Subscript[ϕ, b])-(Subscript[ϕ, q])]  Sin[(Subscript[b, R]) z ])),G[z],z]

Hence: the field lines are given by:
            x=((|β|)/(|q|))  cos((Subscript[ϕ, β])) sec((Subscript[ϕ, β])-(Subscript[ϕ, q])) (((cos((Subscript[ϕ, β])) log(|sin(z (Subscript[β, R]))|))/(|β| cos((Subscript[ϕ, β]))))-z sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(|sin(z (Subscript[β, R]))|)-z |β|cos((Subscript[ϕ, β])) sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(|sin(z (Subscript[β, R]))|)-(Subscript[β, R]) z  sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(|sin(z (Subscript[β, R]) - (Subscript[q, J]) x)|)-(z (Subscript[β, R]) - (Subscript[q, J]) x)  sin((Subscript[ϕ, β]))) + (Subscript[c, 1]) 
            x=(1/(|q| Cos[(Subscript[ϕ, β])-(Subscript[ϕ, q])]))   ( cos((Subscript[ϕ, β])) log(sin(z (Subscript[β, R]) -+ (Subscript[q, J]) x))-(z (Subscript[β, R]) - (Subscript[q, J]) x)  sin((Subscript[ϕ, β]))) + (Subscript[c, 1])

λ = 1; k0=((2 π)/λ); ϵd= 1.4^2 ; ω =2 π ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ 
ClearAll[β, qm, qd, bR, ϕb,  qmI, ϕqm,aqm,  qdI, ϕqd, aqd]  ;
plasmonparm={wp, γ, χ} ;
β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm]]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm]]  ;
bR=Re[β]  ; ϕb=Arg[β]  ;
 qmI=Im[qm]  ;ϕqm=Arg[qm]  ; aqm=Abs[qm]  ;
 qdI=Im[qd]  ;ϕqd=Arg[qd]  ; aqd=Abs[qd]  ;
ContourPlot[z  aqd  Cos[ϕb-ϕqd]-  ( Cos[ϕb] Log[Abs[Sin[x bR- qdI z-ω t]]]-(x bR- qdI  z-ω t)  Sin[ϕb]),
{x,0,0.5},{z,0,0.3}, AspectRatio-> 0.6, Contours-> {0.3,0.7,1.5,2.5,4}, ContourShading-> None, AccuracyGoal-> 8, PlotPoints-> {80,10}] 
,{t,0,1},{λ,0.4,1.24/(0.6*wp)},{{γ,0.02},0.,0.1}]

λ = 0.532; k0=((2 π)/λ); ϵd= 1.48^2 ; ω =2 π ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ 
plasmonparm={wp, γ, χ} ;
ClearAll[β, qm, qd, bR, ϕb,  qmI, ϕqm,aqm,  qdI, ϕqd, aqd]  ;
β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm] ]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm] ]  ;
bR=Re[β]  ; ϕb=Arg[β]  ;
 qmI=Im[qm]  ;ϕqm=Arg[qm]  ; aqm=Abs[qm]  ;
 qdI=Im[qd]  ;ϕqd=Arg[qd]  ; aqd=Abs[qd]  ;
Show[ContourPlot[z  aqd  Cos[ϕb-ϕqd]-  ( Cos[ϕb] Log[Abs[Sin[x bR- qdI z-ω t]]]-(x bR- qdI  z-ω t)  Sin[ϕb]),
{x,0,0.5},{z,0,0.3}, AspectRatio-> 0.6, Contours-> {0.3,0.7,1.5,2.5,4}, ContourShading-> None, AccuracyGoal-> 8,ContourStyle->Thickness[0.01], PlotPoints-> {80,10}] 
, ContourPlot[z  aqm  Cos[ϕb-ϕqm]+  ( Cos[ϕb] Log[Abs[Sin[x bR+ qmI z-ω t]]]-(x bR+ qmI  z-ω t)  Sin[ϕb]),
{x,0,0.5},{z,-0.2,0}, AspectRatio-> 0.4, Contours-> {-0.3,-0.7,-1.5,-2.5,-4}, ContourShading-> None,ContourStyle->Thickness[0.01],AccuracyGoal-> 8, PlotPoints-> {80,10}]
, PlotRange->{{0,0.5},{-0.2,0.3}}, AspectRatio-> 1 ]
,{t,0,1},{λ,0.532,1.24/(0.6*wp)},{{γ,0.02},0.,0.1}]

λ = 1; k0=((2 π)/λ); ϵd= 1.4^2 ; ω =2 π ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ 
plasmonparm={wp, γ, χ} ;
ClearAll[β, qm, qd, bR, ϕb,  qmI, ϕqm,aqm,  qdI, ϕqd, aqd]  ;
β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm] ]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm] ]  ;
bR=Re[β]  ; ϕb=Arg[β]  ;
 qmI=Im[qm]  ;ϕqm=Arg[qm]  ; aqm=Abs[qm]  ;
 qdI=Im[qd]  ;ϕqd=Arg[qd]  ; aqd=Abs[qd]  ;
Show[ContourPlot[z  aqd  Cos[ϕb-ϕqd]-  ( Cos[ϕb] Log[Abs[Sin[x bR- qdI z-ω t]]]-(x bR- qdI  z-ω t)  Sin[ϕb]),
{x,0,0.2},{z,0,0.05}, AspectRatio-> 0.6, Contours-> {0.1,0.2,0.3,0.7,1.2,1.8,2.5}, ContourShading-> None, AccuracyGoal-> 8, PlotPoints-> {80,10}] 
, ContourPlot[z  aqm  Cos[ϕb-ϕqm]+  ( Cos[ϕb] Log[Abs[Sin[x bR+ qmI z-ω t]]]-(x bR+ qmI  z-ω t)  Sin[ϕb]),
{x,0,0.2},{z,-0.5,0}, AspectRatio-> 0.4, Contours-> {-0.1,-0.2,-0.3,-0.7,-1.2,-1.8,-2.5}, ContourShading-> None, AccuracyGoal-> 8, PlotPoints-> {80,10}]
, PlotRange->{{0,0.2},{-0.05,0.05}}, AspectRatio-> 0.5 ]
,{t,0,1},{λ,0.36,1.24/(0.6*wp)},{{γ,0.02},0.,0.1}]

Surface charge plot

EpsMetal[λ_,Parm_List]:=Module[{ω, wp, β, χ},
{wp, β, χ}=Parm ; ω=1.23984187542/λ  ;
1+4 Pi χ -(((wp^2) )/(ω  ( ω +I β )))  ]  ;
B2[ϵd_,ϵm_]:=(( ϵm ϵd)/(ϵm+ ϵd));
Qd[ϵd_,ϵm_]:=Sqrt[-((  ((ϵd)^2))/(ϵm+ ϵd))];
Qm[ϵd_,ϵm_]:= Sqrt[-((  ((ϵm)^2))/(ϵm+ ϵd))] ;

plasmonparm={wp, γ, χ} ;
wp=6. ; γ=0.02; χ=0.001  ;
GraphicsGrid[{{Plot[Re[Qm[1.4^2,EpsMetal[λ,plasmonparm]]] ,{λ,0.25, 0.5}, PlotRange-> All]
,Plot[Im[Qm[1.4^2,EpsMetal[λ,plasmonparm] ]],{λ,0.25, 0.5}, PlotRange-> All]}
, {Plot[Re[Qd[1.4^2,EpsMetal[λ,plasmonparm]]] ,{λ,0.25, 0.5}, PlotRange-> All]
,Plot[Im[Qd[1.4^2,EpsMetal[λ,plasmonparm] ]],{λ,0.25, 0.5}, PlotRange-> All]}}]

E field in metal:  E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (OverHat[ϵ]))) +i ((Subscript[q, m])/((Subscript[k, 0]) (OverHat[ϵ]))) (OverHat[z]) ) (e^(i (β z -ω t))) (e^((Subscript[q, m]) x ))  
E field in dielectric:    E=(Subscript[E, 0]) ( (OverHat[x])  (β/((Subscript[k, 0]) (Subscript[ϵ, d]))) -i ((Subscript[q, d])/((Subscript[k, 0]) (Subscript[ϵ, d]))) (OverHat[z]) )
Surface charge:  2 π σ=(Subscript[D, x,metal])-(Subscript[D, x,dielec])=(Subscript[E, 0]) ( ((1+4 π χ)/(Subscript[ϵ, d])) -1)(β/(Subscript[k, 0])) (e^(i (β z -ω t)))

λ = 1; k0=((2 π)/λ); ϵd= 1.4^2 ; ω =2 π ;
wp=6. ; γ=0.02; χ=0.001  ;
Manipulate[ 
plasmonparm={wp, γ, χ} ;
β=k0  Sqrt[B2[ϵd,EpsMetal[λ,plasmonparm]]] ;
qm=k0 Qm[ϵd,EpsMetal[λ,plasmonparm] ]  ;
qd=k0 Qd[ϵd,EpsMetal[λ,plasmonparm]]  ;
Plot[  Re[ (  (1+4 π χ)/EpsMetal[λ,plasmonparm ] - 1 ) ((β )/k0)(E^(I  ( β  x - ω  t )))]
,{x,0,2},Filling-> Axis
,  PlotLabel->Style["Surface charge",Blue,24]] 
,{t,0,1},{λ,0.5,1.24/(0.6*wp)},{{γ,0.02},0.,0.1}]

2. Plasmon coupling

2.1 Introduction

Surface plasmon itself is just a curiosity if it weren't for resonance effects that we will consider later on. Also, there are people who believed that plasmonics would play an important role in many optical device development.
But a relevant question is that if it is just an evanescent wave confined to a surface, how would it interact with other radiative wave?
First, it is a guided wave, so it can be part of any guided wave system. Below is an example.

2.2 Surface plasmon and waveguide.

Characteristic equation

Let                            H=F[x](e^(i(β z-ω t))) (OverHat[y]). 
(1/c)ϵ ((∂ (Overscript[E, ->]))/(∂ t))  = ∇×(Overscript[H, ->])  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}, {(∂_x), (∂_y), (∂_z)}, {0, F[x], 0}})) = (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t)))
                  ((-i ω)/c)ϵ (Overscript[E, ->])  =  (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{-i β F[x]}, {0}, {(∂_x)F[x]}}))(e^(i (β z-ω t))) 
                    (Overscript[E, ->])  = (1/((Subscript[k, 0]) ϵ)) (({{(OverHat[x]), (OverHat[y]), (OverHat[z])}}))(({{β F[x]}, {0}, {i  (∂_x)F[x]}}))(e^(i (β z-ω t)))

Thus, matching boundary condition: (({{(Subscript[H, y])}, {(Subscript[E, z])}})) to be continuous
1- Across a dielectric
  (({{(Subscript[F, L])[(Subscript[x, B])]}, {(1/(Subscript[ϵ, L]))(∂_x)(Subscript[F, L])[(Subscript[x, B])]}}))= (({{(Subscript[F, R])[(Subscript[x, B])]}, {(1/(Subscript[ϵ, R]))(∂_x)(Subscript[F, R])[(Subscript[x, B])]}}))       (({{(Subscript[F, R])[(Subscript[x, B])]}, {(∂_x)(Subscript[F, R])[(Subscript[x, B])]}}))= (({{1, 0}, {0, ((Subscript[ϵ, R])/(Subscript[ϵ, L]))}}))   (({{(Subscript[F, L])[(Subscript[x, B])]}, {(∂_x)(Subscript[F, L])[(Subscript[x, B])]}}))=(({{1, 0}, {0, ((Subscript[ϵ, R])/(Subscript[ϵ, L]))}}))  S[B;A]  (({{(Subscript[F, L])[(Subscript[x, A])]}, {(∂_x)(Subscript[F, L])[(Subscript[x, A])]}}))   
Consider (Subscript[q, 1]),  (Subscript[q, 2]), 
                               (({{F[x]}, {(∂_x)F[x]}}))=(({{(e^(i k x)), (e^(-i k x))}, {i k (e^(i k x)), -i k (e^(-i k x))}})) . (({{(Subscript[c, 1])}, {(Subscript[c, 2])}})) 
                             (((({{(e^(i k (Subscript[x, A]))), (e^(-i k (Subscript[x, A])))}, {i k (e^(i k (Subscript[x, A]))), -i k (e^(-i k (Subscript[x, A])))}})))^(-1)) . (({{F[(Subscript[x, A])]}, {(∂_x)F[(Subscript[x, A])]}}))=(({{(Subscript[c, 1])}, {(Subscript[c, 2])}})) 
  (({{(Subscript[F, R])[(Subscript[x, B])]}, {(∂_x)(Subscript[F, R])[(Subscript[x, B])]}}))=(({{(e^(i k (Subscript[x, B]))), (e^(-i k (Subscript[x, B])))}, {i k (e^(i k (Subscript[x, B]))), -i k (e^(-i k (Subscript[x, B])))}})) . (((({{(e^(i k (Subscript[x, A]))), (e^(-i k (Subscript[x, A])))}, {i k (e^(i k (Subscript[x, A]))), -i k (e^(-i k (Subscript[x, A])))}})))^(-1)) . (({{F[(Subscript[x, A])]}, {(∂_x)F[(Subscript[x, A])]}})) 
  (({{(Subscript[F, R])[(Subscript[x, B])]}, {(∂_x)(Subscript[F, R])[(Subscript[x, B])]}}))=  S  . (({{F[(Subscript[x, A])]}, {(∂_x)F[(Subscript[x, A])]}}))

(Subscript[S, core])=FullSimplify[(({{(E^(I k L)), (E^(-I k L))}, {I k (E^(I k L)), -I k (E^(-I k L))}}))  . Inverse[(({{1, 1}, {I k , -I k}}))] ]  
(Subscript[S, innerclad])=FullSimplify[(({{(E^((Subscript[q, cl]) d)), (E^(-(Subscript[q, cl]) d))}, { (Subscript[q, cl]) (E^( (Subscript[q, cl]) d)), - (Subscript[q, cl]) (E^(- (Subscript[q, cl]) d))}}))  .  Inverse[(({{1, 1}, {(Subscript[q, cl])  , -(Subscript[q, cl]) }}))] ]

1- Plasmon
 (({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(OverHat[ϵ])) }})) . (({{(e^((Subscript[q, m]) x)), (e^(-(Subscript[q, m]) x))}, {(Subscript[q, m]) (e^((Subscript[q, m]) x)), -(Subscript[q, m]) (e^(-(Subscript[q, m]) x))}}))  .  (({{1}, {0}}))

(({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(OverHat[ϵ])) }})) . (({{(E^((Subscript[q, m]) x)), (E^(-(Subscript[q, m]) x))}, {(Subscript[q, m]) (E^((Subscript[q, m]) x)), -(Subscript[q, m]) (E^(-(Subscript[q, m]) x))}}))  .  (({{1}, {0}}))

(({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(OverHat[ϵ])) }})) . ( (({{(E^((Subscript[q, m]) x)), (E^(-(Subscript[q, m]) x))}, {(Subscript[q, m]) (E^((Subscript[q, m]) x)), -(Subscript[q, m]) (E^(-(Subscript[q, m]) x))}})) /. x-> 0 )  .  (({{1}, {0}}))

2-Interface metal Inner clad  
(Subscript[ (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (({{(e^((Subscript[q, cl]) x)), (e^(-(Subscript[q, cl]) x))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) x)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) x))}})) . (({{(Subscript[c, 1])}, {(Subscript[c, 2])}})) = (({{(Subscript[H, y])}, {(Subscript[E, z])}}))|, inner clad-metal]) =  (({{1, 0}, {0 , (1/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}}))   
   (({{(Subscript[c, 1])}, {(Subscript[c, 2])}}))  =Inverse[ (({{(e^((Subscript[q, cl]) x)), (e^(-(Subscript[q, cl]) x))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) x)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) x))}}))]  . (((({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})))^(-1)) .  (({{1, 0}, {0 , (1/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}}))   
3- Core
(Subscript[ (({{(Subscript[H, y])}, {(Subscript[E, z])}}))|, core-inner clad]) = (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (({{(e^((Subscript[q, cl]) d)), (e^(-(Subscript[q, cl]) d))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) d)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) d))}})) . Inverse[ (({{(e^((Subscript[q, cl]) x)), (e^(-(Subscript[q, cl]) x))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) x)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) x))}}))]  . (((({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})))^(-1)) .  (({{1, 0}, {0 , (1/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}}))   
(Subscript[ (({{(Subscript[H, y])}, {(Subscript[E, z])}}))|, core-inner clad]) = (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}}))   
(Subscript[ (({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) . (({{(e^(i k 0)), (e^(-i k 0))}, {i k (e^(i k 0)), -i k (e^(-i k 0))}}))  (({{(Subscript[c, 3])}, {(Subscript[c, 4])}}))=(({{(Subscript[H, y])}, {(Subscript[E, z])}}))|, core-inner clad]) = (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}}))   
    (({{(Subscript[c, 3])}, {(Subscript[c, 4])}})) =(((({{(e^(i k 0)), (e^(-i k 0))}, {i k (e^(i k 0)), -i k (e^(-i k 0))}})))^(-1)) (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) . (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}}))

3- Core-outer
 (({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) . (({{(e^(i k L)), (e^(-i k L))}, {i k (e^(i k L)), -i k (e^(-i k L))}}))  (({{(Subscript[c, 3])}, {(Subscript[c, 4])}})) 
 (({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) . (({{(e^(i k L)), (e^(-i k L))}, {i k (e^(i k L)), -i k (e^(-i k L))}}))  (((({{(e^(i k 0)), (e^(-i k 0))}, {i k (e^(i k 0)), -i k (e^(-i k 0))}})))^(-1)) (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) . (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}})) 
 (({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) . (Subscript[S, core])  .  (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) .  (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}})) 
    (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (({{(e^((Subscript[q, cl]) 0)), (e^(-(Subscript[q, cl]) 0))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) 0)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) 0))}})) . (({{0}, {(Subscript[c, 6])}}))= (({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) . (Subscript[S, core])  .  (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) .  (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}})) 
  (({{(e^((Subscript[q, cl]) 0)), (e^(-(Subscript[q, cl]) 0))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) 0)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) 0))}})) . (({{0}, {(Subscript[c, 6])}})) =(({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, co])) }})) . (Subscript[S, core])  .  (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) .  (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}})) 
Hence, the characteristic equation is:

(Subscript[S, core])=FullSimplify[(({{(E^(I k L)), (E^(-I k L))}, {I k (E^(I k L)), -I k (E^(-I k L))}}))  . Inverse[(({{1, 1}, {I k , -I k}}))] ]  
(Subscript[S, innerclad])=FullSimplify[(({{(E^((Subscript[q, cl]) d)), (E^(-(Subscript[q, cl]) d))}, { (Subscript[q, cl]) (E^( (Subscript[q, cl]) d)), - (Subscript[q, cl]) (E^(- (Subscript[q, cl]) d))}}))  .  Inverse[(({{1, 1}, {(Subscript[q, cl])  , -(Subscript[q, cl]) }}))] ]

(({{1, 1}, {(Subscript[q, cl])  , -(Subscript[q, cl]) }})) . (({{0}, {(Subscript[C, 6])}}))==(({{1, 0}, {0, ((Subscript[ϵ, cl])/(Subscript[ϵ, co]))}})) .  (Subscript[S, core])  .  (({{1, 0}, {0, ((Subscript[ϵ, co])/(Subscript[ϵ, cl]))}}))  .  (Subscript[S, innerclad])  .  (({{1, 0}, {0, ((Subscript[ϵ, cl])/(OverHat[ϵ]))}}))  . (({{1, 1}, {(Subscript[q, m])  , -(Subscript[q, m]) }})) . (({{1}, {0}}))

FullSimplify[(Subscript[q, cl])  (cos(k L) cosh(d (Subscript[q, cl]))+((sin(k L) sinh(d (Subscript[q, cl])) (Subscript[q, cl]) (Subscript[ϵ, co]))/(k (Subscript[ϵ, cl])))+(((Subscript[q, m]) (Subscript[ϵ, cl]) (((cos(k L) sinh(d (Subscript[q, cl])))/(Subscript[q, cl]))+((cosh(d (Subscript[q, cl])) sin(k L) (Subscript[ϵ, co]))/(k (Subscript[ϵ, cl])))))/(OverHat[ϵ])))] +cos(k L) sinh(d (Subscript[q, cl])) (Subscript[q, cl])+(((Subscript[q, m]) (Subscript[ϵ, cl]) (cos(k L) cosh(d (Subscript[q, cl]))-((k sin(k L) sinh(d (Subscript[q, cl])) (Subscript[ϵ, cl]))/((Subscript[q, cl]) (Subscript[ϵ, co])))))/(OverHat[ϵ]))-((k cosh(d (Subscript[q, cl])) sin(k L) (Subscript[ϵ, cl]))/(Subscript[ϵ, co])) ==0

Expand[(((Subscript[ϵ, cl]) (Subscript[q, m]) (cos(k L) cosh(d (Subscript[q, cl]))-((k (Subscript[ϵ, cl]) sin(k L) sinh(d (Subscript[q, cl])))/((Subscript[q, cl]) (Subscript[ϵ, co])))))/(OverHat[ϵ]))+(1/(k (OverHat[ϵ]) (Subscript[ϵ, cl])))((Subscript[q, cl]) (Subscript[ϵ, co]) sin(k L) ((Subscript[ϵ, cl]) (Subscript[q, m]) cosh(d (Subscript[q, cl]))+(OverHat[ϵ]) (Subscript[q, cl]) sinh(d (Subscript[q, cl])))+k (Subscript[ϵ, cl]) cos(k L) ((Subscript[ϵ, cl]) (Subscript[q, m]) sinh(d (Subscript[q, cl]))+(OverHat[ϵ]) (Subscript[q, cl]) cosh(d (Subscript[q, cl]))))-((k (Subscript[ϵ, cl]) sin(k L) cosh(d (Subscript[q, cl])))/(Subscript[ϵ, co]))+(Subscript[q, cl]) cos(k L) sinh(d (Subscript[q, cl]))]==0

FullSimplify[-((k (Subscript[ϵ, cl]^2) (Subscript[q, m]) sin(k L) sinh(d (Subscript[q, cl])))/((OverHat[ϵ]) (Subscript[q, cl]) (Subscript[ϵ, co])))+(((Subscript[q, cl]) (Subscript[ϵ, co]) (Subscript[q, m]) sin(k L) cosh(d (Subscript[q, cl])))/(k (OverHat[ϵ])))+(((Subscript[q, cl]^2) (Subscript[ϵ, co]) sin(k L) sinh(d (Subscript[q, cl])))/(k (Subscript[ϵ, cl])))-((k (Subscript[ϵ, cl]) sin(k L) cosh(d (Subscript[q, cl])))/(Subscript[ϵ, co]))] +Simplify[(((Subscript[ϵ, cl]) (Subscript[q, m]) cos(k L) sinh(d (Subscript[q, cl])))/(OverHat[ϵ]))+(((Subscript[ϵ, cl]) (Subscript[q, m]) cos(k L) cosh(d (Subscript[q, cl])))/(OverHat[ϵ]))+(Subscript[q, cl]) cos(k L) sinh(d (Subscript[q, cl]))+(Subscript[q, cl]) cos(k L) cosh(d (Subscript[q, cl]))]==0

Expand[((1/(k (OverHat[ϵ]) (Subscript[q, cl]) (Subscript[ϵ, cl]) (Subscript[ϵ, co])))sin(k L) ((Subscript[q, cl]^2) (Subscript[ϵ, co]^2) ((Subscript[ϵ, cl]) (Subscript[q, m]) cosh(d (Subscript[q, cl]))+(OverHat[ϵ]) (Subscript[q, cl]) sinh(d (Subscript[q, cl])))-(k^2) (Subscript[ϵ, cl]^2) ((Subscript[ϵ, cl]) (Subscript[q, m]) sinh(d (Subscript[q, cl]))+(OverHat[ϵ]) (Subscript[q, cl]) cosh(d (Subscript[q, cl]))))+((cos(k L) ((Subscript[ϵ, cl]) (Subscript[q, m])+(OverHat[ϵ]) (Subscript[q, cl])) (sinh(d (Subscript[q, cl]))+cosh(d (Subscript[q, cl]))))/(OverHat[ϵ]))) * k (OverHat[ϵ]) (Subscript[q, cl]) (Subscript[ϵ, cl]) (Subscript[ϵ, co]) ]

Simplify[(Subscript[q, cl]^2) (Subscript[ϵ, cl]) (Subscript[ϵ, co]^2) (Subscript[q, m]) sin(k L) cosh(d (Subscript[q, cl]))+(OverHat[ϵ]) (Subscript[q, cl]^3) (Subscript[ϵ, co]^2) sin(k L) sinh(d (Subscript[q, cl]))-(k^2) (Subscript[ϵ, cl]^3) (Subscript[q, m]) sin(k L) sinh(d (Subscript[q, cl]))-(k^2) (OverHat[ϵ]) (Subscript[q, cl]) (Subscript[ϵ, cl]^2) sin(k L) cosh(d (Subscript[q, cl]))] + Simplify[k (Subscript[q, cl]) (Subscript[ϵ, cl]^2) (Subscript[ϵ, co]) (Subscript[q, m]) cos(k L) sinh(d (Subscript[q, cl]))+k (Subscript[q, cl]) (Subscript[ϵ, cl]^2) (Subscript[ϵ, co]) (Subscript[q, m]) cos(k L) cosh(d (Subscript[q, cl]))+k (OverHat[ϵ]) (Subscript[q, cl]^2) (Subscript[ϵ, cl]) (Subscript[ϵ, co]) cos(k L) sinh(d (Subscript[q, cl]))+k (OverHat[ϵ]) (Subscript[q, cl]^2) (Subscript[ϵ, cl]) (Subscript[ϵ, co]) cos(k L) cosh(d (Subscript[q, cl]))]

sin(k L)  Expand[(Subscript[ϵ, cl]) (Subscript[q, m]) ((Subscript[q, cl]^2) (Subscript[ϵ, co]^2) cosh(d (Subscript[q, cl]))-(k^2) (Subscript[ϵ, cl]^2) sinh(d (Subscript[q, cl])))+(OverHat[ϵ]) ((Subscript[q, cl]^3) (Subscript[ϵ, co]^2) sinh(d (Subscript[q, cl]))-(k^2) (Subscript[q, cl]) (Subscript[ϵ, cl]^2) cosh(d (Subscript[q, cl])))]+k (Subscript[q, cl]) (Subscript[ϵ, cl]) (Subscript[ϵ, co]) cos(k L) ((Subscript[ϵ, cl]) (Subscript[q, m])+(OverHat[ϵ]) (Subscript[q, cl])) (sinh(d (Subscript[q, cl]))+cosh(d (Subscript[q, cl])))==0

sin(k L) (Simplify[(Subscript[q, cl]^2) (Subscript[ϵ, cl]) (Subscript[ϵ, co]^2) (Subscript[q, m]) cosh(d (Subscript[q, cl]))-(k^2) (OverHat[ϵ]) (Subscript[q, cl]) (Subscript[ϵ, cl]^2) cosh(d (Subscript[q, cl]))]+Simplify[(OverHat[ϵ]) (Subscript[q, cl]^3) (Subscript[ϵ, co]^2) sinh(d (Subscript[q, cl]))-(k^2) (Subscript[ϵ, cl]^3) (Subscript[q, m]) sinh(d (Subscript[q, cl]))])+k (Subscript[q, cl]) (Subscript[ϵ, cl]) (Subscript[ϵ, co]) cos(k L) ((Subscript[ϵ, cl]) (Subscript[q, m])+(OverHat[ϵ]) (Subscript[q, cl])) (sinh(d (Subscript[q, cl]))+cosh(d (Subscript[q, cl])))==0

sin(k L) (((sinh(d (Subscript[q, cl])))/(cosh(d (Subscript[q, cl])))) ((OverHat[ϵ]) (Subscript[q, cl]^3) (Subscript[ϵ, co]^2)-(k^2) (Subscript[ϵ, cl]^3) (Subscript[q, m]))+(Subscript[q, cl]) (Subscript[ϵ, cl])  ((Subscript[q, cl]) (Subscript[ϵ, co]^2) (Subscript[q, m])-(k^2) (OverHat[ϵ]) (Subscript[ϵ, cl])))+k (Subscript[q, cl]) (Subscript[ϵ, cl]) (Subscript[ϵ, co]) cos(k L) ((Subscript[ϵ, cl]) (Subscript[q, m])+(OverHat[ϵ]) (Subscript[q, cl])) (((sinh(d (Subscript[q, cl])))/(cosh(d (Subscript[q, cl]))))+1)==0

The characteristic equation is: 
sin(k L) (tanh(d (Subscript[q, cl])) ((OverHat[ϵ]) (Subscript[q, cl]^3) (Subscript[ϵ, co]^2)-(k^2) (Subscript[ϵ, cl]^3) (Subscript[q, m]))+(Subscript[q, cl]) (Subscript[ϵ, cl]) ((Subscript[q, cl]) (Subscript[ϵ, co]^2) (Subscript[q, m])-(k^2) (OverHat[ϵ]) (Subscript[ϵ, cl])))+k (Subscript[q, cl]) (Subscript[ϵ, cl]) (Subscript[ϵ, co]) cos(k L) (tanh(d (Subscript[q, cl]))+1) ((Subscript[ϵ, cl]) (Subscript[q, m])+(OverHat[ϵ]) (Subscript[q, cl]))==0 
neff should contain imaginary component for the characteristic equation, as (OverHat[ϵ]) is complex.

PlasmResolve[neff_,λ_,coreParm_List, cladParm_List, metParm_List]:=
Module[{k0, qcl,k,qm, ϵm, ϵc,L,ϵcl,d} ,
k0=2.  π/λ  ; 
{ϵc,L}=coreParm  ; {ϵcl,d}=cladParm  ; ϵm=EpsMetal[λ,metParm] ;
qcl=Sqrt[(neff^2)-ϵcl]  ; qm=Sqrt[((neff )^2)-ϵm]  ;
k=Sqrt[ϵc-(neff^2)] ;
sin(k L k0) (tanh(d qcl k0) (ϵm  (qcl^3) ((ϵc)^2)-(k^2) ((ϵcl)^3) qm)+qcl ϵcl  ( qcl  ((ϵc)^2) qm-(k^2) ϵm ϵcl))+k qcl ϵcl ϵc  cos(k L k0) (tanh(d qcl k0)+1) (ϵcl  qm+ϵm  qcl) ]

EpsMetal[λ_,Parm_List]:=Module[{ω, wp, β, χ},
{wp, γ, χ}=Parm ; ω=1.23984187542/λ  ;
1+4 Pi χ -(((wp^2) )/(ω  ( ω +I  γ )))  ]  ;
B2[ϵd_,ϵm_]:=(( ϵm ϵd)/(ϵm+ ϵd))  ;
Qd[ϵd_,ϵm_]:=Sqrt[-((  ((ϵd)^2))/(ϵm+ ϵd))]  ;
Qm[ϵd_,ϵm_]:= Sqrt[-((  ((ϵm)^2))/(ϵm+ ϵd))]   ;

plasmonparm={wp, γ, χ} ;
wp=6. ; γ=0.02; χ=0.001  ;
coreParm={1.6,1.} ; cladParm={1.4,1.} ;metParm= {wp, γ, χ}  ;
Plot[{Log[Abs[Re[PlasmResolve[neff,0.7,coreParm, cladParm, metParm]]] ],
,Log[Abs[Im[PlasmResolve[neff,0.7,coreParm, cladParm, metParm]]]]},{neff,1.,1.5}, PlotRange-> All]
Plot[{Log[Re[PlasmResolve[neff,0.7,coreParm, cladParm, metParm]] ^2+
 Im[PlasmResolve[neff,0.7,coreParm, cladParm, metParm]]^2]},{neff,1.,1.5}
, PlotRange-> All]

plasmonparm={wp, γ, χ} ;
wp=6. ; γ=0.02; χ=0.001  ;
coreParm={1.6,1.} ; cladParm={1.4,1.} ;metParm= {wp, γ, χ}  ;
Plot3D[Log[Abs[PlasmResolve[xn + I yn,0.5,coreParm, cladParm, metParm]]] 
,{xn,1.,1.3},{yn,-0.001,0.001} ,PlotRange-> All]

We notice nearly real solution

plasmonparm={wp, γ, χ} ; λ=0.7 ;
wp=6. ; γ=0.02; χ=0.001  ;
core={1.6,1.} ; clad={1.4,1.} ;metParm= {wp, γ, χ}  ;
Manipulate[Plot[{Log[Abs[Re[PlasmResolve[neff,λ,core, clad, metParm]]] ],
,Log[Abs[Im[PlasmResolve[neff,λ,core, clad, metParm]]]]},{neff,nl,nh}, PlotRange-> All],{nl,1.2,1.25},{nh,1.201,1.44}]

plasmonparm={wp, γ, χ} ; λ=0.7 ;
wp=6. ; γ=0.02; χ=0.001  ;
core={2.,3.} ; clad={1.4,2.} ;metParm= {wp, γ, χ}  ;
Manipulate[Plot[{Log[Abs[Re[PlasmResolve[neff,λ,core, clad, metParm]]] ],
,Log[Abs[Im[PlasmResolve[neff,λ,core, clad, metParm]]]]},{neff,nl,nh}, PlotRange-> All],{λ,0.7,0.4},{nl,1.1,1.4},{nh,1.15,1.44}]
 Manipulate[ 
FindRoot[{Re[PlasmResolve[xn+I yn ,λ,core, clad, metParm]]==0, 
Im[PlasmResolve[xn+I yn ,λ,core, clad, metParm]]==0}
, {{xn,x0},{yn,y0}}, AccuracyGoal->15],
{λ,0.7,0.4}, {x0,1.228,1.4},{y0,0.,2}]

Solution set up

(({{1, 1}, {(Subscript[q, m])  , -(Subscript[q, m]) }})) . (({{(Subscript[C, 1])}, {0}}))==(({{1, 0}, {0, ((OverHat[ϵ])/(Subscript[ϵ, cl]))}})) . (Subscript[S, innerclad]) . (({{1, 0}, {0, ((Subscript[ϵ, cl])/(Subscript[ϵ, co]))}}))  .  (Subscript[S, core]) .  (({{1, 0}, {0, ((Subscript[ϵ, co])/(Subscript[ϵ, cl]))}}))  . (({{1, 1}, {(Subscript[q, cl])  , -(Subscript[q, cl]) }})) . (({{0}, {1}}))
1- Metal:
 (({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(Subscript[ϵ, m])) }})) . (({{(e^((Subscript[q, m]) x)), (e^(-(Subscript[q, m]) x))}, {(Subscript[q, m]) (e^((Subscript[q, m]) x)), -(Subscript[q, m]) (e^(-(Subscript[q, m]) x))}})) . (({{1}, {0}}))   
2- At junction:
 (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (({{(e^((Subscript[q, cl]) 0)), (e^(-(Subscript[q, cl]) 0))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) 0)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) 0))}})) . (Subscript[V, in cl])=(({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(Subscript[ϵ, m])) }})) . (({{(e^((Subscript[q, m]) 0)), (e^(-(Subscript[q, m]) 0))}, {(Subscript[q, m]) (e^((Subscript[q, m]) 0)), -(Subscript[q, m]) (e^(-(Subscript[q, m]) 0))}})) . (({{1}, {0}}))   
Hence:
  (Subscript[V, in cl])=Inverse[(({{(e^((Subscript[q, cl]) 0)), (e^(-(Subscript[q, cl]) 0))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) 0)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) 0))}}))]  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, m])) }})) . (({{(e^((Subscript[q, m]) 0)), (e^(-(Subscript[q, m]) 0))}, {(Subscript[q, m]) (e^((Subscript[q, m]) 0)), -(Subscript[q, m]) (e^(-(Subscript[q, m]) 0))}})) . (({{1}, {0}}))

FullSimplify[Inverse[(({{(E^((Subscript[q, cl]) 0)), (E^(-(Subscript[q, cl]) 0))}, {(Subscript[q, cl]) (E^((Subscript[q, cl]) 0)), -(Subscript[q, cl]) (E^(-(Subscript[q, cl]) 0))}}))]  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, m])) }})) . (({{(E^((Subscript[q, m]) 0)), (E^(-(Subscript[q, m]) 0))}, {(Subscript[q, m]) (E^((Subscript[q, m]) 0)), -(Subscript[q, m]) (E^(-(Subscript[q, m]) 0))}}))]

3- At core-inner clad junction  
 (({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) . (({{(e^(i k 0)), (e^(-i k 0))}, {i k (e^(i k 0)), -i k (e^(-i k 0))}})) . (Subscript[V, co])=(({{(Subscript[H, y])}, {(Subscript[E, z])}})) =(({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (({{(e^((Subscript[q, cl]) d)), (e^(-(Subscript[q, cl]) d))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) d)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) d))}}))  . (Subscript[V, in cl])  
  (Subscript[V, co]) =Inverse[(({{(e^(i k 0)), (e^(-i k 0))}, {i k (e^(i k 0)), -i k (e^(-i k 0))}}))]. (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) . (({{(e^((Subscript[q, cl]) d)), (e^(-(Subscript[q, cl]) d))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) d)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) d))}}))  . (Subscript[V, in cl])

FullSimplify[Inverse[(({{(E^(I k 0)), (E^(-I k 0))}, {I k (E^(I k 0)), -I k (E^(-I k 0))}}))]. (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) . (({{(E^((Subscript[q, cl]) d)), (E^(-(Subscript[q, cl]) d))}, {(Subscript[q, cl]) (E^((Subscript[q, cl]) d)), -(Subscript[q, cl]) (E^(-(Subscript[q, cl]) d))}})) ]

4- At core-outer clad junction  
 (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) . (({{1, 1}, {(Subscript[q, cl]) , -(Subscript[q, cl])}}))  . (Subscript[V, out cl])=(({{(Subscript[H, y])}, {(Subscript[E, z])}}))=(({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) . (({{(e^(i k L)), (e^(-i k L))}, {i k (e^(i k L)), -i k (e^(-i k L))}})) . (Subscript[V, co])   
 (Subscript[V, out cl])=Inverse[(({{1, 1}, {(Subscript[q, cl]) , -(Subscript[q, cl])}}))]  .  (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, co])) }})) . (({{(e^(i k L)), (e^(-i k L))}, {i k (e^(i k L)), -i k (e^(-i k L))}})) . (Subscript[V, co])

FullSimplify[Inverse[(({{1, 1}, {(Subscript[q, cl]) , -(Subscript[q, cl])}}))]  .  (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, co])) }})) . (({{(E^(I k L)), (E^(-I k L))}, {I k (E^(I k L)), -I k (E^(-I k L))}}))]

Solution for small core numerical illustration

(({{(e^((Subscript[q, cl]) 0)), (e^(-(Subscript[q, cl]) 0))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) 0)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) 0))}})) . (({{0}, {(Subscript[c, 6])}})) =(({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, co])) }})) . (Subscript[S, core])  .  (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) .  (Subscript[S, clad])  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}})) 
full expansion
  (({{(e^((Subscript[q, cl]) 0)), (e^(-(Subscript[q, cl]) 0))}, {(Subscript[q, cl]) (e^((Subscript[q, cl]) 0)), -(Subscript[q, cl]) (e^(-(Subscript[q, cl]) 0))}})) . (({{0}, {(Subscript[c, 6])}})) =(({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, co])) }})) . (({{(e^(i k L)), (e^(-i k L))}, {i k (e^(i k L)), -i k (e^(-i k L))}}))  . (((({{1, 1}, {i k , -i k}})))^(-1)) .  (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }})) . (({{(e^((Subscript[q, cl]) d)), (e^(-(Subscript[q, cl]) d))}, { (Subscript[q, cl]) (e^( (Subscript[q, cl]) d)), - (Subscript[q, cl]) (e^(- (Subscript[q, cl]) d))}}))
  . (((({{1, 1}, {(Subscript[q, cl])  , -(Subscript[q, cl]) }})))^(-1))  . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  (({{1}, {0}})) 
Hence:    V1=(({{1}, {0}}))
Innerclad:          (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) (({{1, 1}, {(Subscript[q, cl])  , -(Subscript[q, cl]) }})) . V2=(({{1, 0}, {0 , (1/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  V1    
                                V2=(((({{1, 1}, {(Subscript[q, cl])  , -(Subscript[q, cl]) }})))^(-1))  (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(OverHat[ϵ])) }})) . (({{1, 1}, {(Subscript[q, m]), -(Subscript[q, m]) }}))  .  V1    
Core:            (({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }})) .  (({{1, 1}, {i k , -i k}}))  . V3=(({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }}))  . (({{(e^((Subscript[q, cl]) d)), (e^(-(Subscript[q, cl]) d))}, { (Subscript[q, cl]) (e^( (Subscript[q, cl]) d)), - (Subscript[q, cl]) (e^(- (Subscript[q, cl]) d))}}))  .  V2    
                       V3=(((({{1, 1}, {i k , -i k}})))^(-1)) .  (({{1, 0}, {0 , ((Subscript[ϵ, co])/(Subscript[ϵ, cl])) }}))  . (({{(e^((Subscript[q, cl]) d)), (e^(-(Subscript[q, cl]) d))}, { (Subscript[q, cl]) (e^( (Subscript[q, cl]) d)), - (Subscript[q, cl]) (e^(- (Subscript[q, cl]) d))}}))  .  V2    
Outerclad:          (({{1, 0}, {0 , (1/(Subscript[ϵ, cl])) }})) .  (({{1, 1}, {(Subscript[q, cl]) , -(Subscript[q, cl]) }}))  . V4=(({{1, 0}, {0 , (1/(Subscript[ϵ, co])) }}))  . (({{(e^(i k L)), (e^(-i k L))}, {i k (e^(i k L)), -i k (e^(-i k L))}}))  .  V3    
                          V4=(((({{1, 1}, {(Subscript[q, cl]) , -(Subscript[q, cl]) }})))^(-1)) . (({{1, 0}, {0 , ((Subscript[ϵ, cl])/(Subscript[ϵ, co])) }}))  . (({{(e^(i k L)), (e^(-i k L))}, {i k (e^(i k L)), -i k (e^(-i k L))}}))  .  V3

PlasmFieldParmDB[neff_,λ_,coreParm_List, cladParm_List, metParm_List]:=
Module[{k0, qcl,k,qm, ϵm, ϵc,L,ϵcl,d,SM,V, Coef } ,
k0=2.  π/λ  ; 
{ϵc,L}=coreParm  ; {ϵcl,d}=cladParm  ; ϵm=EpsMetal[λ,metParm] ;
qcl=Sqrt[(neff^2)-ϵcl]  ; qm=Sqrt[((neff )^2)-ϵm]  ;
k=Sqrt[ϵc-(neff^2)] ;
SM={Inverse[(({{1, 1}, {qcl k0  , -qcl k0 }}))]  .  (({{1, 0}, {0 , ((ϵcl)/(ϵm)) }})) . (({{1, 1}, {qm k0, -qm k0 }}))
,  Inverse[(({{1, 1}, {I k k0, -I k k0}}))] .  (({{1, 0}, {0 , ((ϵc)/(ϵcl)) }}))  .  (({{(E^(d qcl k0)), (E^(-d qcl k0))}, {qcl k0 (E^(d qcl k0)), - qcl k0 (E^(-d qcl k0))}}))
,Inverse[(({{1, 1}, {qcl k0  , -qcl k0 }}))] . (({{1, 0}, {0 , ((ϵcl)/(ϵc)) }})) .  (({{(E^(I k k0 L)), (E^(-I k k0 L))}, {I k k0 (E^(I k k0 L)), -I k k0 (E^(-I k k0 L))}}))
 } ;
V={1.,0. }  ;   Coef ={V}   ;
For[i=1,i≤Length[SM] , i++, V= SM[[i]]   .  V   ; AppendTo[Coef, V]]   ;
{{ϵc,ϵcl,ϵm,k0,k, qcl, qm , L,d},Coef} ] ;

PlasmFieldParm[neff_,λ_,coreParm_List, cladParm_List, metParm_List]:=
Module[{k0, qcl,k,qm, ϵm, ϵc,L,ϵcl,d,SM,V, Coef } ,
k0=2.  π/λ  ; 
{ϵc,L}=coreParm  ; {ϵcl,d}=cladParm  ; ϵm=EpsMetal[λ,metParm] ;
qcl=Sqrt[(neff^2)-ϵcl]  ; qm=Sqrt[((neff )^2)-ϵm]  ;
k=Sqrt[ϵc-(neff^2)] ;
SM={ (1/2)(({{((qm ϵcl)/( qcl ϵm))+1, 1-((qm ϵcl)/(qcl ϵm))}, {1-((qm ϵcl)/(qcl ϵm)), ((qm ϵcl)/(qcl ϵm))+1}})) 
,(({{(E^(d k0 qcl)) (k ϵcl-I qcl ϵc), (E^(-d k0 qcl)) (k ϵcl+I qcl ϵc)}, {(E^(d k0 qcl)) (k ϵcl+I qcl ϵc), (E^(-d k0 qcl)) (k ϵcl-I qcl ϵc)}}))/(2  k ϵcl) 
,(({{(E^(I k k0 L)) (qcl ϵc+I k ϵcl), (E^(-I k k0 L)) (qcl ϵc-I k ϵcl)}, {(E^(I k k0 L)) (qcl ϵc-I k ϵcl), (E^(-I k k0 L)) (qcl ϵc+I k ϵcl)}}))/(2 qcl ϵc)} ;
V={1.,0. }  ;   Coef ={V}   ;
For[i=1,i≤Length[SM] , i++, V= SM[[i]]   .  V   ; AppendTo[Coef, V]]   ;
{{ϵc,ϵcl,ϵm,k0,k, qcl, qm , L,d},Coef} ] ;

PlasmField[x_,PFieldParm_List]:=  Module[{ϵc,ϵcl,ϵm,k0, k, qcl,qm,L,d, Coef} ,
{{ϵc,ϵcl,ϵm, k0, k, qcl,qm,L,d}, Coef}= PFieldParm  ;
If[x≤   0., (1/(ϵm))   (({{(E^(qm k0 x )), (E^(-qm k0 x))}, {qm  k0 (E^(qm k0 x)), -qm  k0 (E^(-qm k0 x))}})) . Coef[[1]]
, If[x≤ d,  (1/(ϵcl))   (({{(E^(qcl k0 x)), (E^(-qcl k0 x))}, {qcl k0 (E^(qcl k0 x)), -qcl  k0 (E^(-qcl k0 x))}})) . Coef[[2]]
,If[x≤ L+d,  (1/(ϵc))  (({{(E^(I k k0  (x-d))), (E^(-I k k0  (x-d)))}, {I k k0 (E^(I k k0  (x-d))), -I k  k0 (E^(-I k k0  (x-d)))}})) . Coef[[3]]
,  (1/(ϵcl))  (({{(E^(qcl k0 (x-d-L))), (E^(-qcl k0 (x-d-L)))}, {qcl  k0 (E^(qcl k0 (x-d-L))), -qcl  k0 (E^(-qcl k0 (x-d-L)))}})) . Coef[[4]]
]
]
]
]

EpsMetal[λ_,Parm_List]:=Module[{ω, wp, β, χ},
{wp, γ, χ}=Parm ; ω=1.23984187542/λ  ;
1+4 Pi χ -(((wp^2) )/(ω  ( ω +I  γ )))  ]  ;
B2[ϵd_,ϵm_]:=(( ϵm ϵd)/(ϵm+ ϵd))  ;
Qd[ϵd_,ϵm_]:=Sqrt[-((  ((ϵd)^2))/(ϵm+ ϵd))]  ;
Qm[ϵd_,ϵm_]:= Sqrt[-((  ((ϵm)^2))/(ϵm+ ϵd))]   ;

plasmonparm={wp, γ, χ} ; 
wp=6. ; γ=0.02; χ=0.001  ;
coreParm={1.6,1.} ; cladParm={1.4,0.01} ;metParm= {wp, γ, χ}  ;
Manipulate[Plot[Log[Abs[PlasmResolve[neff,λ,coreParm, cladParm, metParm]]] ,{neff,nl,nh}, PlotRange-> All]
,{{λ,0.5},0.7,0.4},{nl,1.1,1.4},{nh,1.15,1.44}]
 Manipulate[ λ=wl ;
{nR,nI} ={xn, yn}   /. FindRoot[{Re[PlasmResolve[xn+I yn ,λ,coreParm, cladParm, metParm]]==0
,Im[PlasmResolve[xn+I yn ,λ,coreParm, cladParm, metParm]]==0}
, {{xn,x0},{yn,y0}}, AccuracyGoal->15],
{{wl,0.5},0.7,0.4}, {x0,1.1,1.4},{y0,0.,2}]

λ

neff=nR+ I nI 
PlasmResolve[neff,λ,coreParm, cladParm, metParm]
PFParm=PlasmFieldParm[neff,λ ,coreParm, cladParm, metParm]

ModeP=Plot[-Re[PlasmField[x,PFParm][[1]]],{x,-0.5,4}, PlotRange-> All]

What happen to this mode if we don't have the metal layer?

plasmonparm={wp, γ, χ} ; 
coreParm={1.6,1.} ; cladParm={1.4,0.01} ;
wp=0. ; γ=0.02; χ=(1.4-1)/(4 π )  ;
metParm= {wp, γ, χ}  ;
Manipulate[Plot[Log[Abs[PlasmResolve[neff,λ,coreParm, cladParm, metParm]]] ,{neff,nl,nh}, PlotRange-> All]
,{{λ,0.5},0.7,0.4},{nl,1.1,1.4},{nh,1.15,1.44}]
 Manipulate[ λ=wl ;
{nR,nI} ={xn, yn}   /. FindRoot[{Re[PlasmResolve[xn+I yn ,λ,coreParm, cladParm, metParm]]==0
,Im[PlasmResolve[xn+I yn ,λ,coreParm, cladParm, metParm]]==0}
, {{xn,x0},{yn,y0}}, AccuracyGoal->15],
{{wl,0.5},0.7,0.4}, {x0,1.1,1.4},{y0,0.,2}]

neff=nR+ I nI 
PlasmResolve[neff,λ,coreParm, cladParm, metParm]
PFParm=PlasmFieldParm[neff,λ ,coreParm, cladParm, metParm]

modeNoP=Plot[Re[PlasmField[x,PFParm][[1]]],{x,-0.5,4}, PlotRange-> All]

Show[modeNoP, ModeP]

3. Plasmon resonance phenomena

This section cannot be recovered.