Computer Lab 9: Lecture 26

As with all of the exercises you will complete in this course, you will be expected to turn your *.nb file in to the course TA before the next lecture.  Please do so by e-mailing your completed exercise file to Lan Yu (ece329.uiuc@gmail.com), make sure you name the file using the following notation: “Exn_Last_First.nb”, where n is the Exercise number (in this case, 26), and First and Last would be your first and last name.  So my submission to Lan for this exercise would be titled “Ex26_Wasserman_Dan.nb”.  You will receive full credit for your exercise if it turned in before the next lecture (Lecture n+1).  No late submissions will be accepted.

Objective/Overview:

In class and our exercises last week, we looked at signal/carrier waves, polarization (circular, elliptical, and linear), and last, but certainly not least, reflection of E&M waves at an interface between two materials.  This last topic is hugely important for a range of applications, including optics, telecommunications, RADAR, LIDAR and even transmission lines.  In the exercise 25 you looked at what the waves incident upon a perfect dielectric looked like, and at the end, introduced a conductivity to the second material.  Today, we are going to dive a little deeper into the behavior of waves incident on a conducting medium.

Setting Up the System

You already wrote the code for an E&M wave incident upon a conducting surface, so we are in good shape to start today’s work.  Let’s bring that code into our program, again assuming  that our wave is incident on a surface at z=0, and is traveling in the +(OverHat[z]) direction, with a frequency ω=2π*(10^14)rad/s in vacuum with (Subscript[μ, 1])=(Subscript[μ, o]) and (Subscript[ϵ, 1])=(Subscript[ϵ, o]) and amplitude (Subscript[E, i])=(OverHat[x])1V/m.  We will allow our (Subscript[μ, 2]),(Subscript[ϵ, 2]) and (Subscript[σ, 2]) to be inputs:

etot[mu2_,eps2_,sigma2_]:=(muo=4*Pi*10^-7;epso=8.85*10^-12;w=2*Pi*10^14;
etao=Sqrt[muo/epso];
bo=w*Sqrt[muo*epso];g=Sqrt[I*w*muo*mu2*(sigma2+I*w*epso*eps2)];eta2s=Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)];
gammas=(eta2s-etao)/(etao+eta2s);taus=1+gammas;etxs=If[z>0,Re[taus*Exp[-g*z]*Exp[I*w*t]],0];
erxs=If[z<0,Re[gammas*Exp[I*bo*z]*Exp[I*w*t]],0];
eixs=If[z<0,Re[Exp[-I*bo*z]*Exp[I*w*t]],0];Animate[Plot[(eixs+etxs+erxs)/.{t->to},{z,-10*10^-6,10*10^-6}, PlotRange->{-2,2},PlotLabel->{gammas,taus,etao,eta2s}],{to,0,10^-13}])

You’ll notice that I am also printing to the Plot a PlotLabel that includes the reflection and transmission coefficients, as well as the characteristic impedance of free space and our conducting material.  Let’s write a Manipulate command, and look at what happens to our wave as sigma goes to infinity, holding (Subscript[μ, 2])=(Subscript[ϵ, 2])=1:

Manipulate[etot[mu2,eps2,sigma2],{mu2,1},{eps2,1},{sigma2,0},ControlType->InputField]

At (Subscript[σ, 2])~1000 [S/m] we start to see our wave decaying upon incidence on the interface.  However, because (Subscript[μ, 2])=(Subscript[ϵ, 2])=1, we don’t see much ‘bouncing’ of the wave at z<0.  This is because, if you look at our impedances (Subscript[η, o])=(Subscript[η, 1]) and (Subscript[η, 2]), you see that they are very close to equal, with only a small change in the real part of η as we cross the interface.  As we begin to further increase (Subscript[σ, 2]), our wave starts to die out much fast in z>0, and we see our ‘bouncing’ wave at z<0.  Let’s keep pushing (Subscript[σ, 2])...By the time we get to (Subscript[σ, 2])~(10^7), our z<0 wave looks like a stationary Sin wave oscillating in space.  You can also see that our (Subscript[η, 2]) is getting closer and closer to 0, and because of this, our τ~0 and our Γ~-1, as you would expect.
What we are doing here is approximating the behavior of a perfect conductor ((Subscript[σ, 2])=∞), and the resulting E&M wave is called a “standing wave”. Analytically we see that for the case of (Subscript[σ, 2])=∞, (Subscript[η, 2])=0, and therefore  Γ=-1 and τ=1+Γ=0.  In phasor notation, we can write the waves in z<0 as:
(OverTilde[E])=(OverHat[x])(Subscript[E, o])[(e^(-(Subscript[jβ, o])z))-(e^(+(Subscript[jβ, o])z))]=-(OverHat[x])2(Subscript[E, o])sin((Subscript[β, o])z)
If we add the time dependence (e^(jωt))  and take the Re part of the product of the time dependence and the phasor, we get
(OverVector[E])(z,t)=-(OverHat[x])2(Subscript[E, o])Sin((Subscript[β, o])z)Sin(ωt), in other words, a Sin wave in space that oscillates with a radial frequency ω, which is just what we see in our Animation.

Let’s look at the H field for the same system:

htot[mu2_,eps2_,sigma2_]:=(muo=4*Pi*10^-7;epso=8.85*10^-12;w=2*Pi*10^14;
etao=Sqrt[muo/epso];
bo=w*Sqrt[muo*epso];g=Sqrt[I*w*muo*mu2*(sigma2+I*w*epso*eps2)];eta2s=Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)];
gammas=(eta2s-etao)/(etao+eta2s);taus=1+gammas;htxs=If[z>0,Re[(taus/eta2s)*Exp[-g*z]*Exp[I*w*t]],0];
hrxs=If[z<0,Re[-(gammas/etao)*Exp[I*bo*z]*Exp[I*w*t]],0];
hixs=If[z<0,Re[(1/etao)*Exp[-I*bo*z]*Exp[I*w*t]],0];Animate[Plot[(hixs+htxs+hrxs)/.{t->to},{z,-10*10^-6,10*10^-6}, PlotRange->{-.01,.01},PlotLabel->{gammas,taus,etao,eta2s}],{to,0,10^-13}])

Manipulate[htot[mu2,eps2,sigma2],{mu2,1},{eps2,1},{sigma2,0},ControlType->InputField]

Hmmmm....in some ways this looks very similar to our E field, but if we look closely, we see that the H field is not an oscillating Sin wave, but an oscillating Cos wave!  This is a very important distinction, because our E field, as (Subscript[σ, 2])->∞ goes to 0 at z=0, but our H field actually has a maximum at z=0!!  We can perform the same analytical calculations for H, and we see:
(OverTilde[H])=(OverHat[y])((Subscript[E, o])/(Subscript[η, o]))[(e^(-(Subscript[jβ, o])z))+(e^(+(Subscript[jβ, o])z))]=-(OverHat[y])2((Subscript[E, o])/(Subscript[η, o]))Cos((Subscript[β, o])z)
If we add the time dependence (e^(jωt))  and take the Re part of the product of the time dependence and the phasor, we get
(OverVector[H])(z,t)=(OverHat[y])2((Subscript[E, o])/(Subscript[η, o]))Cos((Subscript[β, o])z)Cos(ωt), in other words, a Cos wave in space that oscillates with a radial frequency ω, which is just what we see in our Animation.

Surface Currents on Conducting Interfaces

Let’s think about boundary conditions a little bit.  Because our fields are all parallel to the interface, we only have to worry about the tangential boundary contitions (though this would change for wave incident at an angle). We know that 
(OverHat[n])×((SuperPlus[E])-(SuperMinus[E]))=0, which holds because we know that our (OverVector[E]) is zero at the interface for a perfect conductor, and we can see from our animation that regardless of our choice of (Subscript[μ, 2]), (Subscript[ϵ, 2]) and (Subscript[σ, 2]), the E field is always constant across the interface.
If we look at our H field what happens?  If (Subscript[σ, 2])=0, our H field is constant at z=0.  If we add a (Subscript[σ, 2]), it remains constant across z=0, just like the E field, though what we see is that the H field drops VERY quickly to a value of 0 as we move past z=0.  If we go to   (Subscript[σ, 2])->∞, then the H field drops over an infinitesimal region of z (essentially a thickness of 0)!  If we look at our boundary conditions for H:
(OverHat[n])×((SuperPlus[H])-(SuperMinus[H]))=(Subscript[J, s]), we see that the only way that this sudden drop can happen is if we have a surface current!  Because we know what the H field is for the z<0 side of z=0: (OverVector[H])(z->(SuperMinus[0]),t)=(OverHat[y])2((Subscript[E, o])/(Subscript[η, o]))Cos((Subscript[β, o])z->(SuperMinus[0]))Cos(ωt)=(OverHat[y])2((Subscript[E, o])/(Subscript[η, o]))Cos(ωt) and using
(Subscript[J, s])=(OverHat[n])×((SuperPlus[H])-(SuperMinus[H]))=(OverHat[z])×(0-(SuperMinus[H]))=(OverHat[x])2((Subscript[E, o])/(Subscript[η, o]))Cos(ωt)[A/m]

This is pretty interesting!  And can help us to understand reflection from a perfect conductor in a physical manner:
An incident E field sloshes electrons in our material around, the electrons move to cancel out the incident electric field.
These electrons oscillate, creating a surface current (Subscript[J, s])~Cos(ωt)
We know that a surface current creates an H field propagating away from the surface (as well as an E field).
The E field generated by (Subscript[J, s]) cancels out the incident E field at z=0, and the H field enhances the incident H field at z=0.

What about an “imperfect” conductor? Well, there will be a current here too, and we actually know how to find it, because we know that (OverVector[J])=σ·(OverVector[E]).  So let’s plot our current density in a conducting material.  We can basically just cut and paste from our previous work:

Clear[jz,j]

jz[mu2_,eps2_,sigma2_]:=(g=Sqrt[I*w*muo*mu2*(sigma2+I*w*epso*eps2)];eta2s=Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)];gammas=(eta2s-etao)/(etao+eta2s);taus=1+gammas;etxs=If[z>0,Re[taus*Exp[-g*z]*Exp[I*w*t]],0];
j=sigma2*etxs;Animate[Plot[j/.{t->to},{z,-1*10^-6,10*10^-6}, PlotRange->{-sigma2,sigma2},PlotLabel->{gammas,taus,etao,eta2s}],{to,0,10^-13}])

Manipulate[jz[mu2,eps2,sigma2],{mu2,1},{eps2,1},{sigma2,0},ControlType->InputField]

OK, we can see that we have a decaying, oscillating current.  We can write the current phasor as the product of our (Subscript[(OverTilde[E]), t]) and our (Subscript[σ, 2]):

ClearAll[jphas,jphasf,jphastot]

jphas[mu2_,eps2_,sigma2_]:=(sigma2*(2*Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)])/(etao+Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)]))*Exp[-Sqrt[I*w*muo*mu2*(sigma2+I*w*epso*eps2)]*z]

Let’s integrate the current density phasor over the region z>0:

jphastot[mu2_,eps2_,sigma2_]:=Integrate[jphas[mu2,eps2,sigma2],{z,0,Infinity}]

jphastot[1,1,10000]

If we have a perfect conductor, we know that (Subscript[J, s])=2((Subscript[E, o])/(Subscript[η, o]))Cos(ωt), if (Subscript[E, o])=1, then our surface current has an amplitude:

2/etao

Let’s go back to our integral of the phasor and find the amplitude of the phasor:

Abs[jphastot[1,1,10000000]]

We can see that for a good conductor (but not perfect) if we integrate over all of the current density, we end up with an equivalent surface current of about the same magnitude as we get for the surface current of a perfect conductor (the analytical derivation of this result is in the class notes).

Energy

So what about energy in our system?  We know how much energy is incident upon the surface:

siavg=0.5*Re[Exp[-I*bo*z]*(1/etao)*Exp[I*bo*z]]

We can also calculate the reflected energy:

sravgf[mu2_,eps2_,sigma2_]:=(g=Sqrt[I*w*muo*mu2*(sigma2+I*w*epso*eps2)];eta2s=Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)];gammas=(eta2s-etao)/(etao+eta2s);taus=1+gammas;sravg=0.5*Re[gammas*Conjugate[-(gammas/etao)]])

sravgf[1,1,1000000]

Where does the rest of the energy go?  Well, the one thing we haven’t considered is ohmic heating.  Remember, when we wrote out the Poynting Theorem, we were thinking about the (OverVector[J])·(OverVector[E]) term as a source, because (OverVector[J])·(OverVector[E])<0. For z>0 we now also have a (OverVector[J])·(OverVector[E]) term, but it will be positive, it is acting as a sink of energy!
The total power absorbed by a conducting material will be the integral of the (OverVector[J])·(OverVector[E]) from z=0 to z=∞, which we can calculate!

jeavgf[mu2_,eps2_,sigma2_]:=(g=Sqrt[I*w*muo*mu2*(sigma2+I*w*epso*eps2)];eta2s=Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)];gammas=(eta2s-etao)/(etao+eta2s);taus=1+gammas;jez=0.5*Re[sigma2*taus*Exp[-g*z]*Conjugate[taus*Exp[-g*z]]];jeavg=Integrate[jez,{z,0,Infinity}])

jeavgf[1,1,1000000]

So if we add up the reflected and absorbed and incident power (being careful with signs) we see:

powers[mu2_,eps2_,sigma2_]:={siavg,sravgf[mu2,eps2,sigma2],-jeavgf[mu2,eps2,sigma2]}

powers[1,1,10000]

totalpower[mu2_,eps2_,sigma2_]:=siavg+sravgf[mu2,eps2,sigma2]-jeavgf[mu2,eps2,sigma2]

totalpower[1,1,1000000000]

Tada!  We have accounted for all of the energy flow in our system...except...what happens to the total power when (Subscript[σ, 2])=0.  Does this violate conservation of energy?

Radiation Pressure

There’s one more interesting phenomenon we should cover (if we have enough time) today.  If you think about what is going on at the surface of our conductor, you will see that we have a current (which is made up of moving charges) and an H field.  From the Lorentz Force law, this means we should have a force acting on our surface, since (Subscript[(OverVector[F]), B])=q(OverVector[v])×(OverVector[B]).  If we have N electrons moving, then our force is (Subscript[(OverVector[F]), B])=Nq(OverVector[v])×(OverVector[B])=Nqμ(OverVector[v])×(OverVector[H]), and we can recognize that (OverVector[J])=Nq(OverVector[v]), so the incident light exerts a force (Subscript[(OverVector[F]), B])=(OverVector[J])×(OverVector[B])=(OverVector[J])×μ(OverVector[H]) on our surface.  This is actually a Force/Unit Area, in the units we are using, which is Pressure.  This force is called Radiation Pressure and is a real, measurable (though not easily) quantity.
Use the approaches we have covered above to determine the radiation pressure on a conducting surface from an incident wave as a function of (Subscript[μ, 2]), (Subscript[ϵ, 2]), and (Subscript[σ, 2]).
<(Subscript[P, rad])≥∫0.5*Re[(OverTilde[J])(z)×(OverTilde[(SuperStar[B])])(z)]dz

fb[mu2_,eps2_,sigma2_]:=(g=Sqrt[I*w*muo*mu2*(sigma2+I*w*epso*eps2)];eta2s=Sqrt[I*w*muo*mu2/(sigma2+I*w*epso*eps2)];gammas=(eta2s-etao)/(etao+eta2s);taus=1+gammas;etxs=If[z>0,Re[taus*Exp[-g*z]*Exp[I*w*t]],0];
htxs=If[z>0,Re[(taus/eta2s)*Exp[-g*z]*Exp[I*w*t]],0];
jxh=
forceb=Integrate[jxh,{z,0,Infinity}])

fb[1,1,0]