Light Propagation (part 1)
Spatial Behavior

ECE 5368/6358                  han q le - copyrighted
Use solely for students registered for UH ECE 6358/5368 during courses - DO NOT distributed (copyrighted materials)

Segment 3

A digression: a review of the wave concept

Link to periodic and wave lecture

5. Spatial aspects of light waves

5.1 General wave and harmonic wave

5.1.1 General wave - time-dependence

Let’s simplify equations  (4.1.10)  (4.1.11):
                          (∇^2)U -((μ ϵ)/(c^2))((∂^2)/(∂ (t^2)))U=0,           
to just one-dimension to see what kind of solutions we get: 
                         ((∂^2)/(∂ (x^2)))U -((μ ϵ)/(c^2))((∂^2)/(∂ (t^2)))U=0,          (5.1.1)
This has the general wave solution:
                            f[x±v t] where v =(c/(Sqrt[μ ϵ])).  (5.1.2)
Thus, we can similarly reduce (∇^2)U if we can make coordinate reduction. One such reduction is x.(OverHat[s]) where (OverHat[s]) is a unit vector:  
(∇^2)f[x.(OverHat[s])±v t] -((μ ϵ)/(c^2))((∂^2)/(∂ (t^2)))f[x.(OverHat[s])±v t]= (((Subscript[s, x])^2) +((Subscript[s, y])^2) +((Subscript[s, z])^2)) (((d^2) F)/(d (ξ^2)))-(v^2)((μ ϵ)/(c^2)) (((d^2) F)/(d (ξ^2))) =0      (5.1.3)
                                                  =  (((d^2) F)/(d (ξ^2)))-(v^2)((μ ϵ)/(c^2)) (((d^2) F)/(d (ξ^2))) =(1-(v^2)((μ ϵ)/(c^2)))  (((d^2) F)/(d (ξ^2)))=0
which is satisfied by:                 1-(v^2)((μ ϵ)/(c^2))=0   or:  v=(c/(Sqrt[μ ϵ]))

Looks like we have solved the wave equation! Does that mean we can solve all EM wave problems? 
Well, not really. First, remember that the wave equation (∇^2)U -((μ ϵ)/(c^2))((∂^2)/(∂ (t^2)))U=0 is not strictly correct in a medium. It is correct for vaccuum and it is OK to accept the above as a general wave equation. But a generic form  f[x±v t] is NOT very useful for particular problems, in which we need more details about the function f . But for now, it is useful just for the sake of illustration of the wave concept.

The above equation is called plane wave because every point on the plane x.(OverHat[s]) = constant has the same value at a given time. The entire plane has the same value and it moves at the speed v. Here is an illustration:

f[x_]:=Sin[π x]/(π  x) ;
Plot[f[x], {x,-10,10}, PlotRange-> All
,Frame-> True, GridLines-> Automatic,Filling-> Axis
,FillingStyle-> Hue[0.05,1,1]]

Demonstration:

Demo: ENVELOPE WAVE PULSE

Here is an example of an impulse wave

https://www.youtube.com/watch?v=NXqzO5g_juw

Is the general solution above sufficient to describe everything now?

Let’s look at a laser beam, how do you describe the beam 100 m from here, or a 1 M km from here? how about the beam from a flash light? How about the beam through a slit, a pin hole? a grating? How about the light coming out from a diode laser? why does one direction diverge than the other? how does that divergence depend on the laser structure? Afterall, all those beams are propagating in free space and satisfy the wave equation. Don’t we already have a solution?

Well, not all waves are plane wave! In fact, no real life wave is truly a plane wave! Even if we use plane wave in some problems, it is ONLY an approximation that is good enough for the problem of interest. Plane wave is an idealized model of waves. It is very useful, but clearly NOT enough to really describe so many other light waves that we deal around us and in many optical systems.
We need specific solutions of waves for various problems of interest.

Is this a solution: ∫g[(OverHat[s])]f[x.(OverHat[s])±v t]d(OverHat[s]) ?  Which space does the integration extend?

5.1.2 Harmonic wave or monochromatic wave

The above general wave solution is not the only way to solve the wave equation. Let’s take another common approach: Solving by part.

(∇^2)U -((μ ϵ)/(c^2))((∂^2)/(∂ (t^2)))U=0 

We split U into 2 components:
                                   U=F[x] T[t].                               (5.1.4a)
Then:                       T[t] (∇^2)F= -((μ ϵ)/(c^2))F[x](((∂^2) T)/(∂ (t^2)))                (5.1.4b)
or                            (1/(F[x]))(∇^2)F= -((μ ϵ)/(c^2))(1/(T[t] ))(((∂^2) T)/(∂ (t^2)))= constant C       (5.1.4c)

DSolve[(1/(T[t] ))(∂_(t,t))T[t]==Cn,T[t],t]

What should (Sqrt[Cn]) be?

Why don't we want solution where (Sqrt[Cn]) is real? Does it make physical sense?

Let T[t] = (e^(-i ω t)).  This is a fundamental harmonic wave: F[x] (e^(-i ω t)). (See 3.2)

We are interested in this type of solutions because lightwaves are fundamentally harmonic waves. We will see later on that a photon has a unique frequency, and all lightwaves are simply made up of many photons (which all can have different frequencies).
Thus, we are interested in learning specifically harmonic waves: waves that oscillate periodically.

If we substitute T[t] = (e^(-i ω t)) in the equation above, then:
                                     (∇^2)F+ ((ω^(2 ))/(c^2))μ ϵ F=0:                     (5.1.5)
this is Helmholtz equation. 
As mentioned above, it is the frequency-domain approach to the Maxwell’s theory and it is more relevant to the real problem of optical media. It appears as if (5.1.1) is derived from the wave equations (4.1.10) and (4.1.11) and if the formers are not strictly correct, then how can (5.1.1) be correct? The answer is that (5.1.1) can be derived independently from (4.1.10) and (4.1.11) in frequency domain approach to the whole Maxwell’s equation shown in (3.2)

In the following, we will consider aspects of the spatial component F[x] of the wave, rather than the temporal component, which we will assume to be always (e^(-i ω t)). 
The essence is to find solutions to the Helmholtz equation.
We will study waves with different values of frequency ω in later chapters.

Demo: Visible light wave - frequency and physiological perception of color

5.2 Harmonic plane wave

We see that a very easily recognized solution to the Helmholtz equation is the spatial harmonic functions: sin and cos. Or, we use the complex function which is more convenient than sin or cos:
                                                  F[x]=(e^(i x. k  ))           (5.2.1)
(of course, this is NOT the only solution to the Helmholtz eq.).

Combine what we learn above from the general form of wave with harmonic wave, we have harmonic plane waves:
             F[x]=(e^(±i x. k  ))   and   U=F[x] T[t]=(e^(i (±x. k - ω t) ))     (5.2.2)

Is there any relationship between k and ω?

Substitute in the Helmholtz equation:
                 -(k^2) + ((ω^(2 ))/(c^2))μ ϵ = 0, or k = ((ω )/c)(Sqrt[μ ϵ]) = (Subscript[k, 0])(Sqrt[μ ϵ])    (5.2.3)

Link to review of harmonic function and wave

Embedded

1. Definition of harmonic function

The fundamental functions that exhibit periodicity are sine and cosine functions, which are also called harmonic functions. The basic form is:     X[t]=A Cos[2 π f t+ϕ]≡A Cos[ω t+ϕ]  
where                              ω≡ 2 π f : Angular (or radian) frequency
A: amplitude, f: frequency, ϕ: phase, and t is the function variable

2. Illustration

3. Comparing Sine and Cosine

Manipulate[Plot[  {Cos[2 π t+ϕc π],Sin[2 π t+ϕs  π]}, {t,0,4}, PlotRange-> {-1,1},Frame-> True , GridLines-> Automatic
,PlotLabel-> {"red: cos ; blue: sine"}
, PlotStyle-> {{Thick,Red}, {Thick,Blue}}
,Filling-> Axis
],{{ϕc,0,"Cos phase (π)"},-2,2 },{{ϕs,0,"Sine phase (π)"},-2,2 }]

4. Harmonic wave: both space and time are involved - example 1: sound wave:

Manipulate[Plot[A Sin[2 π f (x-v t)],{x,0, 2}, Filling-> Axis,PlotRange-> {-1,1}]
,{{A,1,"amplitude"},0,1},{{f,1.5,"frequency"},1,3},{{v,1,"velocity"},-1,1} ,{{t,0,"time"},0,1}]

5. More sound wave illustration:

Manipulate[f=(10^(0.1*Lf)) ; {f,
Sound[Play[Cos[2 π f t],{t,0,1},SampleRate->22050]]},
{{Lf,13,"Frequency (dBHz)"},13,40,0.1}]

Sound[{SoundNote["C"],SoundNote["C"],SoundNote["G"],SoundNote["G"],SoundNote["A"],SoundNote["A"],SoundNote["G"]},{0,4}]

For sound wave:  X[t]=A Cos[2 π f (x-v t)+ϕ] :
A: amplitude tells us how loud the sound is (or how intense the sound is). We can hear over a wide range of sound amplitude
- See link: http://hyperphysics.phy-astr.gsu.edu/hbase/sound/earsens.html
 f: frequency (or wavelength) tells us the pitch. We are most sensitive to the variation of this, which is why we humans have music! Just think how much we human love music to recognize how important it is. 
v: the speed of sound which is the effect in which we hear thunder long after lightning. We humans are aware of the finite speed of sound whenever we hear echo from a large structure: e.g. hill, mountain, cave, building, etc. 
ϕ: phase: we humans are not sensitive to sound phase. But this is a very important aspect in active anti-noise technology.

See link on anti-noise technology

6. Light wave:

Manipulate[{Plot[A Sin[2 π f (x-v t)],{x,0, 2}, Filling-> Axis,PlotRange-> {-1,1}]
,Graphics[{Hue[(f-1)*0.4,1,(A^2)],Disk[{0,0},1] }] } 
,{{A,1,"amplitude"},0,1},{{f,1.5,"frequency"},1,3},{{v,1,"velocity"},-1,1} ,{{t,0,"time"},0,1}]

For light wave:  X[t]=A Cos[2 π f (x-v t)+ϕ] :
A: amplitude tells us how bright the light is.
 f: frequency (or wavelength) tells us what color the light is (note: human perception of color based on our retina 3-color receptors is NOT one-to-one to pure single frequency, but each single frequency gives us a definitive color perception if it is within our visible wavelength range). Of course, we still refer to f as “color” in the broad sense of the word, even if we can’t see it such as mm-wave, infrared or UV, X-Ray...
v: the speed of light which is a “fundamental” constant of nature in vacuum, but it is different in matter. It is too fast for us to feel it! But we can see the effects on the astronomical scale, such as Doppler red-shift of receding galaxies that tells us expanding universe (Edwin Hubble)
ϕ: phase: we humans are sensitive only to lightwave amplitude squared and its frequency/wavelength but not phase. We can tell bright light vs. dim light, light of different color, but can’t tell if they have difference phase. But phase is a very important if not crucial aspect of light involving in all optical phenomena around us. In fact, to understand how something as intrinsic our eyes work, we need to understand how the phase changes as light goes through the lens of the eye. The engineering of light phase is the essence of numerous most important modern optical technologies.

5.2.1 Illustration and discussion of harmonic plane wave

Link to plane wave

So now we have harmonic plane waves. Is that enough for us to answer question above? (about laser beam etc...)

Let's look at the above solution again. We write:
(e^(i x. k  ))=(e^(i ((Subscript[k, x]) x + (Subscript[k, y]) y+(Subscript[k, z]) z ) )), we can think of this as the product of three waves: x, y, and z.
The requirement is: (k^2) =((Subscript[k, x])^2) + ((Subscript[k, y])^2) + ((Subscript[k, z])^2)  = ((ω^(2 ))/(c^2))μ ϵ= ((Subscript[k, 0])^2)μ ϵ 
What if we accidentally choose ((Subscript[k, x])^2) + ((Subscript[k, y])^2) >((Subscript[k, 0])^2)μ ϵ ? Can we have a solution?

5.2.2 Evanescent wave

Let's consider only x and z for simplicity ((Subscript[k, y])=0). We know this:
                                   ((Subscript[k, x])^2)+((Subscript[k, z])^2)  =  ((Subscript[k, 0])^2)μ ϵ =(((((ω )/c)))^2) μ ϵ.        (5.2.4)
Suppose ((Subscript[k, x])^2)  >  ((Subscript[k, 0])^2)μ ϵ, what is the solution?
Clearly,                          ((Subscript[k, z])^2)  =  ((Subscript[k, 0])^2)μ ϵ -  ((Subscript[k, x])^2) < 0.            (5.2.5)
                                      (Subscript[k, z])  = (Sqrt[((Subscript[k, 0])^2)μ ϵ -  ((Subscript[k, x])^2)])  =i κ.       (5.2.6a)
where                             κ = (Sqrt[((Subscript[k, x])^2)-((Subscript[k, 0])^2)μ ϵ ])                      (5.2.6b)
How about this solution:
                                   (e^(i ((Subscript[k, x]) x +(Subscript[k, z]) z ) ))=(e^(i (Subscript[k, x]) x ))(e^(-κ z ))                 (5.2.7)

Demo evanescent wave

Is it a physically meaningful wave? Where do you think they exist? What happens at z -> - ∞

The above wave is called evanescent wave in the z dimension. Discussion in class.

5.3 Wave in a medium: complex dielectric constant

In the above example, we see that plane wave can be made up of purely real and purely imaginary propagation constant (e^(i (Subscript[k, x]) x ))(e^(-κ z )). More generally, can we have wave (e^(i x. k))  where k is complex? Why not? we see that the only requirement is that
                      ((Subscript[k, x])^2) + ((Subscript[k, y])^2) + ((Subscript[k, z])^2)  = ((ω^(2 ))/(c^2))μ ϵ= ((Subscript[k, 0])^2)μ ϵ            (5.3.1)
- Case 1: we see that if ((Subscript[k, 0])^2)μ ϵ is real, but any component ((Subscript[k, x])^2), ((Subscript[k, y])^2), ((Subscript[k, z])^2) or the sum of 2 therof is greater than ((Subscript[k, 0])^2)μ ϵ, the remaining component(s) can be complex.
- Case 2: if a medium has a complex μ ϵ, then ((Subscript[k, 0])^2)μ ϵ  is complex and hence some or all of (Subscript[k, x]), (Subscript[k, y]), (Subscript[k, z]) have to be complex. 
The 2 cases are different however, and we must make distinction.

Fundamentally, what is complex μ, ϵ?
Remember: classical optics wave theory is ONLY a description of lightwave behavior based on macroscopic average or approximation of the interaction of between photons and matters (electrons in atoms, molecules,...). 
As it turns out, many of these interactions result in a simplified description of lightwave in media with quantity μ, ϵ that are complex.
A specifically and particularly important class of phenomena is absorption and/or stimulated emission. 
We are familiar with absorption all around us.
We will learn about stimulated emission and optical gain later in discussing about lasers.
In these phenomena, lightwave behavior can be well described with a complex dielectric constant.
But... remember again: the real fundamental theory here is the quantum interaction between photons and atoms, molecules. There is no fundamental meaning in complex ϵ other than that.
In other media, nanoscopic structures can also result in complex optical behaviors, e. g. metamaterials, negative index of refraction etc.

For Case 2, we write ϵ= (Subscript[ϵ, R]) +i (Subscript[ϵ, I ]) or ϵ= (Subscript[ϵ, 1]) +i (Subscript[ϵ, 2]).
If (Subscript[ϵ, I]) > 0, (k^2) = ((((Subscript[k, R])+i (Subscript[k, I])))^2)= ((Subscript[k, 0])^2)μ ((Subscript[ϵ, R]) +i (Subscript[ϵ, I])) entails that (Subscript[k, I])>0; then the wave:
(e^(i x ((Subscript[k, R])+i (Subscript[k, I])))) =(e^(i x (Subscript[k, R]))) (e^(- (Subscript[k, I]) x)) has both a spatially periodic term (wave) and a exponentially decaying term in the SAME direction. It means that the wave gets weaker as it propagates. 
The energy must be lost in the medium somehow (unless there is a reflection wave - but here, we don’t see any reflection wave, hence the energy is loss in the medium). On the reverse, the wave can gain energy from a medium.
In general, we don't have to make any distinction for a lossy or amplifying medium. All wave parameters are assumed to have complex component whenever physically possible.

Illustration

Make distinction between these 2 cases:

Demo abs/evanesc wave comparison

5.4 Wave in other geometries: spherical and cylindrical

5.4.1 Spherical wave in isotropic media

Let's start from the Helmholtz equation:
                                                     (∇^2)F+ ((Subscript[k, 0])^2)ϵ μ F=0. 
Earlier we write: (∇^2)=   (∂_(x,x)) + (∂_(y,y))+ (∂_(z,z)). What if we use a different coordinate?

Laplacian[F[r, θ, ϕ],{r,θ,ϕ},"Spherical"]

Let’s consider a wave that is purely spherical, it has no dependence on the angles:

Laplacian[F[r],{r,θ,ϕ},"Spherical"]

k0= .;
DSolve[((2 F'(r))/r)+F''(r)+(k0^2) ϵ μ F[r]==0,F[r],r]

The solution is
                                                ((e^(-i k r  ))/( r ))(Subscript[c, 1])+((e^(i k r  ))/( r ))(Subscript[c, 2])                (5.4.1)
where                                         k=(Subscript[k, 0])(Sqrt[ϵ μ])
So, we have a spherical wave solution also: ((e^(±i k r- i ω t  ))/( r )). What do we do with this solution? What is the physical meaning? what happens at r = 0?

Demo spherical wave 1

Let’s cut across the equator plane and plot the wave amplitude in density plot:

Demo spherical wave 2

Demo spherical wave 3

If the wave decays in amplitude, where does its energy (or power) go?

Demo spherical wavefront and source distance

Where (in which phenomena) do you expect to see this type of wave?

Remember: like plane wave, an ideal or perfect spherical wave doesn’t exist in nature. But they can be very useful approximation to numerous problems. For example, if we consider a very small source of light from a distance much larger than the extent of the source (far away) and many (a large number) wavelengths, it is as good as a spherical wave for any practical purpose calculation. The point source can be just a small aperture, or even a single atom or molecule emitting the photon. 
Of course each photon has a finite momentum k, hence a specific direction. But if k is purely random with  equal probability in any direction, then the atom or the molecule is practically a source of spherical wave.
We will see that in dealing with various problems in diffraction, spherical wave is a useful model.

Spherical wave in non-isotropic media can be quite more complicated (mathematically)

5.4.2 Cylindrical wave

We choose a different coord syst and we get an answer not expected from plane wave solutions. What if we choose another one like cylindrical?

Laplacian[F[ρ], {ρ,θ,z},"Cylindrical" ]

DSolve[F''(ρ)+((F'(ρ))/ρ)+(k0^2) ϵ  F[ρ]==0,F[ρ],ρ]

Again, it can be an approximated model for some some waves in some problem.

The solution now is Bessel J and Y function!. What are other coordinates? what are other solutions?

Demo cylindrical wave

Notice that it behaves like a standing wave! which means that it must be a sum of incoming and outgoing wave. That's the only way to have finiteness at the center. Can we construct a radiating wave?

Where (in which phenomena) do you expect to see this type of waves?

5.5 Summary

There are many many (infinite) solutions to the Helmholtz wave equation. Many of them, if not all, have physical meaning one way or another. Our objective is to understand, gain insight, and have the ability to select the right solution (including most appropriate approximated solutions) to the physical problem of our interest. Use various wave models intelligently to do what we want with the given problems.

6. Some key light wave parameters

6.1 Review basic concepts, power and energy

Light carries energy as we know. We learn from modern (quantum) physics that light consists of quantized package energy known as photons. (we will have a review of quantum physics). Each photon has an energy:
                                                     E=h ν                        (6.1.1)
where h is the Planck's constant.
If we have a spatial density of (Subscript[N, ph])[ν] photons in a medium with the speed (c/(Sqrt[ϵ μ])), the power crossing a surface area S in the direction of light propagation is:
                                       P=(lim_(Δt-> 0))(1/(Δt))(E S (c/(Sqrt[ϵ μ]))Δt)=(c/(Sqrt[ϵ μ])) (Subscript[N, ph])[ν] h ν S          (6.1.2)
In the limit S-> 0, we have power density or intensity:
                                         (lim_(ΔS-> 0))((ΔP)/(ΔS)) =I=(c/(Sqrt[ϵ μ])) (Subscript[N, ph])[ν] h ν              (6.1.3) 
Here, we study classical wave optics, where we deal with a large number of photons, not having to worry about the quantum behavior of light because we can treat them statistically as ensemble average. Their ensemble statistical behavior is described by the Maxwell's equations. Thus, in analogy with Eqs above:

The instantaneous energy density of an EM field is generally given by:
                                         u[x;t]=(1/(8 π))(ϵ (E^2)[x;t]+μ (H^2)[x;t])     (6.1.4)
For an EM wave, the energy is carried by the wave and we do not think of energy density as something stationary like for static E and H field in a capacitor or a Helmholtz coil. As EM field forms a wave, the energy flows with the wave. For harmonic field, the physically meaningful quantity is the time-averaged energy:
                                          ⟨u[x]⟩=(1/(16 π))(ϵ ⟨(E^2)[x]⟩+μ ⟨(H^2)[x]⟩)     (6.1.5)
over one cycle.

For a plane wave or approximately a plane wave at an infinitesimal location:
                                            H[x]=(Sqrt[(ϵ/μ)])(OverHat[k])×E[x]               (6.1.5)
hence:                          ⟨u[x]⟩=(1/(16 π))(ϵ ⟨(E^2)[x]⟩+ϵ ⟨(E^2)[x]⟩)=(1/(8 π))ϵ ⟨(E^2)[x]⟩     (6.1.6)     
Similar to (6.1.3) we have:     
                                            I=((c )/(Sqrt[ϵ μ]))⟨u[x]⟩=((c )/(8 π (Sqrt[ϵ μ])))ϵ ⟨(E^2)[x]⟩ 
Or                                           I= (c/(8π))(Sqrt[(ϵ/μ)]) ⟨(E^2)[x]⟩                   (6.1.7) 
which, we will see can be described with Poynting vector.

Excersise: What is the energy density of a plane wave

Excersise: What is the power density of a plane wave

6.2 Poynting vector

Poynting vector is defined as the energy flux per unit area:
                                               S=(c/(4 π))E×H                        (6.2.1)
(note: use real representatives of the fields, i. e. S=(c/(4 π))Re[E]×Re[H] ). 
But instantaneous Poynting vector has no signigicant physical meaning, we care only its time-average (over one period).
It can be shown that the Poynting vector represents the power density (((dP)/(dS))) or intensity I of a light beam by taking its time-averaged, which, for harmonic wave, is:
                                         ⟨(Overscript[S, →])⟩=(c/(8 π))Re[ (Overscript[E, ->]) ×(Overscript[(SuperStar[H]), ->])]                (6.2.2a)
(the brackets ⟨ ⟩ here is used to denote time-averaged quantity).
                                        I=|⟨(Overscript[S, →])⟩|=(c/(8 π))|Re[ (Overscript[E, ->]) ×(Overscript[(SuperStar[H]), ->])]|          (6.2.2b)
Plane wave in a medium:       I=|⟨(Overscript[S, →])⟩|=|(c/(8 π))Re[ (Overscript[E, ->]) ×(Overscript[(SuperStar[H]), ->])]|   
                                               I=(c/(8 π))(Sqrt[(ϵ/μ)])((|(Overscript[E, ->])|)^2)                  (6.2.2c)

A best illustrative example is to apply the Poyting vector concept to a plane wave (HW).

Calculate the Poynting vector for a plane wave

See supplementary materials on plane wave.

6.3 A note on polarization (more in a different chapter)

Remember that EM wave is vectorial in nature. It is not a scalar wave. Hence, although we generally care more about the absolute magnitude of the field rather than individual components such as (Subscript[E, x]),  (Subscript[E, y]), (Subscript[E, z]). The orientation of E and H vectors area also an important parameters of a light wave. Polarization is a concept used to describe this orientation. We will study this in another chapter.

7. Principle of linear superposition (spatial)

7.1 Discussion - concept

We know that light waves all around us, in various technologies are certainly more complicated than those plane wave, spherical or cylindrical waves above. How do we describe these waves?
Just like in differential equation, we know that there are infinite solutions, but all of them can be constructed from the basis set of solutions, the same is applicable to light waves. The solutions we discuss above are “basis” solution and using the linear superposition principle, we can construct any light wave.

Example:

DSolve[y''[x] +(k^2) y[x]==0, y[x],x]

Source code

Manipulate[ Plot[{Sin[2 π x],Cos[2 π x], q . {Sin[2 π x],Cos[2 π x]} },{x,0,1},
PlotStyle->{ {Thick,Red},{Thick,Blue}, {Thick,Black}}, PlotRange-> {-1.5,1.5}
,Filling->{{3-> {Axis, Hue[0.8,1,1,0.2]}}}
,Frame-> True, GridLines-> Automatic,FrameLabel-> {Style["x",20 ],Style["y",20]}  ,LabelStyle-> {20, FontFamily-> "Arial"}]
,{{q,{0.5,0.5}},{-1,-1},{1,1}}]

The black curve shows that any solution to the differential equation is a linear combination of the red (Sin) and blue (Cos) basis functions.

Coefficients (Subscript[c, 1]) and (Subscript[c, 2]) are arbitrary and we can construct infinite solutions. Sin and Cos are basis functions. In many systems, the basis functions are orthogonal, and they can be normalized to become orthonormal.

A basis set may have an infinite number of orthonormal functions, such as Sin[n ω x], Cos[ n ω x], n= 0,1, 2,..
An example of linear combination of an infinite number of basis function is the Fourier series theorem, which allows us to construct any arbitrary-shape periodic function.

For discrete basis set, we write a general function as a sum:
                      Ψ[x,t]=(∑_(n=1)^∞)(Subscript[c, n]) (Subscript[ψ, n])[x,t]                   (7.1.1)
For continuous basis set, we use integration: 
                     Ψ[x,t]=(∫_□^□)c[k] (Subscript[ψ, n])[k;x,t] dk          (7.1.2)

How many plan wave solutions do we have for the Helmholtz Eq.?

(∇^2)F+ ((ω^(2 ))/(c^2))μ ϵ F=0       ;    F[x] =(e^(i (x. k - ω t) ))

The answer is: continuously infinite.
We know that plane waves form a continuous set, so in principle, we can form a solution like this for a monochromatic wave
             (∫_□^□)c[k] (e^(i (x. k-ω t))) dk =(e^(-i ω t))(∫_□^□)c[k] (e^(i x. k)) dk         (7.1.3)
So, this is the principle, let’s implement it in an example:

7.2 Example of linear superposition: 2D Gaussian beam

Link to Assignment

Suppose we have a monochromatic 2D light beam in x and z. The beam is constant along y direction. For example, if we put a screen perpendicular to z (i. e. x, y plane), it looks like this:

Plot[(E^(- 2 (x^2)/(w^2))) /. w-> 1 ,{x,-3,3}, Filling-> Axis, FillingStyle-> Blue]

Source code

Manipulate[Show[DensityPlot[(E^(- 2 (x^2)/(w^2)))  ,{y,-10,10},{x,-3,3},ColorFunction->GrayLevel, FrameLabel->{y,x}, AspectRatio->0.3 ]
,ParametricPlot[{11+4(E^(- 2 (x^2)/(w^2))),x},{x,-3,3}]
,Graphics[{Red,Thickness[0.01],Dashed,Line[{{0,-5},{0.,5}}]}]
,PlotRange-> {{-10.,15},{-3,3}}
]

,{{w,1.,Button["Beamwaist w",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.344,1.,0.634,1]]},0.5,3}]

This is what we call “Gaussian” profile (along x-direction). How do we describe the E field wave E[x,z,t] of this light beam? (the above give us ONLY a profile in x at z=0, and time-averaged.

Since it is monochromatic, we know that the time factor is simply (e^(-i ω t)). It’s up to us to find out the z part.

We know that  (k^2) =((Subscript[k, x])^2) + ((Subscript[k, z])^2)  = ((ω^(2 ))/(c^2))μ ϵ= ((Subscript[k, 0])^2)μ ϵ . So only one of the 2 k's can be independent. Let's pick (Subscript[k, x]) to be independent and the integrating variable. Then, the formula (7.1.3) becomes:
             (e^(-i ω t))(∫_□^□)c[(Subscript[k, x])] (e^(i (x (Subscript[k, x]) + z (Subscript[k, z])))) d(Subscript[k, x])                     (7.2.1a)
                (e^(-i ω t))(∫_(-∞)^∞)c[(Subscript[k, x])] (e^(i (x (Subscript[k, x]) + z (Sqrt[((Subscript[k, 0])^2)μ ϵ- ((Subscript[k, x])^2)])))) d(Subscript[k, x])      (7.2.1b)
Now, what we really need is c[(Subscript[k, x])]. This is the “arbitrary” coefficient that allows us to construct any solution. How do we choose c[(Subscript[k, x])]  in this case?

7.2.1 Integral expression for Gaussian profile

For the Gaussian profile above:

Plot[(E^(- (x^2)/(w^2))) /. w->1 ,{x,-3,3}, Filling-> Axis, FillingStyle-> Blue]

This is what we have:  (the solution at z=0) = (e^(- (x^2)/(w^2))).
We want:  (∫_(-∞)^∞)c[(Subscript[k, x])] (e^(i (x (Subscript[k, x]) + z (Sqrt[((Subscript[k, 0])^2)μ ϵ - ((Subscript[k, x])^2)])))) d(Subscript[k, x]) at z = 0 = (e^(- (x^2)/(w^2)))
Or:                                     ~(∫_(-∞)^∞)c[(Subscript[k, x])] (e^(i x (Subscript[k, x]) )) d(Subscript[k, x]) = (e^(- (x^2)/(w^2)))     (7.2.2)        
How do we solve for c[(Subscript[k, x])] ? We need inverse Fourier transform:
                     c[(Subscript[k, x])]=(1/(2 π))(∫_(-∞)^∞) (e^(- (x^2)/(w^2))) (e^(-i x (Subscript[k, x]) ))dx               (7.2.3)

k= . ;
(1/(2 π))(∫_(-∞)^∞) (E^(-I x k )) (E^(- (x^2)/(w^2)))\[DifferentialD]x

So the solution is                     c[(Subscript[k, x])] = (( w)/(2 (Sqrt[π])))(e^(-(1/4) ((Subscript[k, x])^2) (w^2) ))   (7.2.4)    
What is the meaning of  c[(Subscript[k, x])]?

c[k] describes the beam profile in terms of a "package" of several beams with different spatial frequency. Key concept: spatial frequency k. c[k] is the envelop or coefficient of the "package". (we'll study a time wave package later on). It can also be called the "angular spectrum".

Now we have our solution, complete with the time factor:
            E[x,z,t]= A (( w)/(2 (Sqrt[π]))) (e^(-i ω t)) (∫_(-∞)^∞)(e^(-(1/4) (k^2) (w^2))) (e^(i (x k + z (Sqrt[((Subscript[k, 0])^2) ϵ- (k^2)])))) dk   (7.2.5)

This cannot be integrated in a closed form, but we can numerically evaluate it, or later, do an approximation.

7.2.2 Graphic illustration

Numerical integration code:

Source code

Beam[x_,z_, n_,λ_,w_]:=Module[{kn,bR,bE }, kn=((2  π)/λ)n ;
bR=w *NIntegrate[(E^(I  k x ))  (E^(I  (Sqrt[(kn^2)-(k^2)]) z ))   (E^(-(1/4) (k^2) (w^2)))   ,{k,-Min[kn,4./w], Min[kn,4./w]}]  ;
bE= If[ kn < 4./w, 
		 2* w * NIntegrate[Re[(E^(I  k x )) (E^(I  (Sqrt[(kn^2)-(k^2)]) z )) ]  (E^(-(1/4) (k^2) (w^2)))   ,{k,kn ,4./w}] ,
		                                                                   0] ; {bR,bE} ]

We can plot the E field:

λ=0.6328 ;n=1;w= 0.5 ;
beamplot=Plot3D[ Re[Beam[x,z,n ,λ ,w][[1]] ],{x,-5,5},{z,-5,5},PlotPoints->{25,50},Mesh->False,FaceGrids->All, BoxRatios-> {1,1,1}, PlotRange-> All]

This is one half of it:

7.2.3 Paraxial approximation

Let’s find an approximation to do the ingration for the beam. Let's look at the integral:
                                             (∫_(-∞)^∞)(e^(-(1/4) (k^2) (w^2)))(e^(i  k x ))(e^(i  (Sqrt[((Subscript[k, 0])^2)ϵ-(k^2)]) z ))dk
If z is very large (far field) and the beam is paraxial (mostly in the propagation axis), we can think that (Subscript[k, z]) is the dominant term and (Subscript[k, x]) is small. Then, we can approximate:
(Sqrt[((Subscript[k, 0])^2)ϵ-(k^2)]) =(Sqrt[((Subscript[k, n])^2)-(k^2)]) =(Subscript[k, n])(Sqrt[1-((k^2)/((Subscript[k, n])^2))]) ≈(Subscript[k, n]) (1-((k^2)/(2 ((Subscript[k, n])^2)))) = (Subscript[k, n]) - ((k^2)/(2 (Subscript[k, n]))) 
(∫_(-∞)^∞)(e^(-(1/4) (k^2) (w^2)))(e^(i  k x ))(e^(i  (Sqrt[((Subscript[k, 0])^2)ϵ-(k^2)]) z ))dk ≈ (∫_(-∞)^∞)(e^(-(1/4) (k^2) (w^2)))(e^(i  k x ))(e^(i  (Subscript[k, n]) z -i ((k^2)/(2 (Subscript[k, n]))) z ))dk 
                                                           = (e^(i  (Subscript[k, n]) z  ))(∫_(-∞)^∞)(e^(-(1/4) (k^2) (w^2)))(e^(i  k x ))(e^( -i ((k^2)/(2 (Subscript[k, n]))) z ))dk
                                                           = (e^(i  (Subscript[k, n]) z  ))(∫_(-∞)^∞)(e^(-(1/4) (k^2) ((w^2) + i ((2 z)/( (Subscript[k, n])))  ) ))(e^(i  k x ))dk

Now, we can perform integration. Define   P=((w^2) + i ((2 z)/( (Subscript[k, n])))  )

(∫_(-∞)^∞)(E^(-(1/4) (k^2) P))(E^(I  k x ))\[DifferentialD]k

Hence:

w=.  ;n = . ; z= . ;
(E^(I  (Subscript[k, n]) z  )) ((2 (Sqrt[π]) (E^(-((x^2)/P))))/(Sqrt[P])) /. P-> ((w^2) + I ((2 z)/( (Subscript[k, n])))  )

Thus:  E[x,z,t]= A (( w)/(2 (Sqrt[π])))  ((2 (Sqrt[π]) (e^(i z (Subscript[k, n])-((x^2)/((w^2)+((2 i z)/(Subscript[k, n])))))))/(Sqrt[(w^2)+((2 i z)/(Subscript[k, n]))])) (e^(-i ω t)) =A (( 1)/(Sqrt[1+((2 i z)/((w^2) (Subscript[k, n])))]))   (e^(i z (Subscript[k, n])-((x^2)/((w^2)+((2 i z)/(Subscript[k, n])))))) (e^(-i ω t))    (7.2.6)

Show[,BoxRatios-> {1,1,0.1}]

We can do some more mathematical manipulation to give us some insight. We can rewrite:

E[x,z,t]=A (( 1)/(Sqrt[1+((2 i z)/((w^2) (Subscript[k, n])))]))   (e^(i z (Subscript[k, n])-((x^2)/((w^2)+((2 i z)/(Subscript[k, n])))))) (e^(-i ω t))    (7.2.6)
                                          E[x,z,t]=A (( w)/(W[z]))   (e^(i (z (Subscript[k, n])-ω t)))  (e^(-((x^2)/((W[z])^2))))         (7.2.7a)
 where                                W[z]=(Sqrt[(w^2)+((2 i z)/( (Subscript[k, n])))]) ;  ((W[z])^2)=(w^2)+((2 i z)/( (Subscript[k, n])))   (7.2.7b)
 then, this is what we see: 
 -The part:  (e^(i (z (Subscript[k, n])-ω t))) is just the plane wave in the z direction.
 - The part (e^(-((x^2)/((W[z])^2)))) is a Gaussian profile in the x-dimension that modifies the plane wave part (e^(i (z (Subscript[k, n])-ω t))). It has a Gaussian waist ((W[z])^2)
 - The only thing “unsual” is that the Gaussian waist is z-dependence and complex: ((W[z])^2)=(w^2)+((2 i z)/( (Subscript[k, n])))
 We will learn more about the “meaning” of complex beam waist!

7.3 Illustration of some properties of 2D Gaussian beam

7.3.1 Phase front.

The most obvious thing we see in the Gaussian beam calculation above is the circular “wave front”

The key insight to the wavefront is in this term: (e^(i z (Subscript[k, n])-((x^2)/((w^2)+((2 i z)/(Subscript[k, n]))))))   (7.3.1) 
Let’s do some reduction to separate the amplitude and the phase:

i z (Subscript[k, n])-((x^2)/((w^2)+((2 i z)/(Subscript[k, n])))) =  i  (Subscript[k, n]) ( z -((x^2)/(i (Subscript[k, n]) (w^2)-2 z))) = i  (Subscript[k, n]) ( z +(((x^2) (2 z+i (Subscript[k, n]) (w^2)))/((2 z-i (Subscript[k, n]) (w^2))(2 z+i (Subscript[k, n]) (w^2)))))
                          = i  (Subscript[k, n]) ( z +(((x^2) 2 z)/(4 (z^2)+ ((Subscript[k, n])^2) (w^4)))) -(x^2) (( ((Subscript[k, n ])^2)  (w^2))/(4 (z^2)+ ((Subscript[k, n])^2) (w^4)))
                          = i  (Subscript[k, n]) ( z +((x^(2 ))/(2 z+(( ((Subscript[k, n])^2) (w^4))/(2 z))))) -((x^2)/(w^2)) (( 1)/(1+(((((2 z)/((Subscript[k, n]) (w^2)))))^2)))    ≈  i  (Subscript[k, n]) ( z +((x^2)/(2 z))) -((x^2)/(w^2)) (( 1)/(1+(((((2 z)/((Subscript[k, n]) (w^2)))))^2) ))
But:  i  (Subscript[k, n]) ( z +((x^2)/(2 z))) ≈i  (Subscript[k, n])  (Sqrt[(z^2)+ (x^2)])=i  (Subscript[k, n])  r  
This is a circular (cylindrical) wavefront indeed with the term (e^(i (Subscript[k, n]) r ))

λ=0.6328 ;n=1;w= 0.5 ; kn=2 π n/λ  ;
ParametricPlot3D[{r Cos[θ], r Sin[θ], Re[(E^(I   kn r ))]},{r,0,5},{θ,0,2 π}, Mesh-> False, PlotPoints-> {40,15}, BoxRatios->{1,1,0.1}]

Let’s compare:

A natural concept in describing wave is "phase front". We commonly think of "wave front" in layman term. It's the same. A phase front is a surface of all points that have the same continuous phase. Mathematically, if we express a wave as:
(e^(i ϕ[x,y,z])) (where ϕ[x,y,z] is real), then, any surface defined by ϕ[x,y,z]=constant is a phase front. (other name: cophasal surfaces, wave surfaces)

Notice the phasefront, but what else do we need?

7.3.2 Beam divergence

For the Gaussian beam, the amplitude is NOT circularly uniform, but it is rather distribute in the center lobe:

We need a term that make a “lobe” in the center along the propagation and the lobe should have a constant angular divergence.

We now can examine the second term (amplitude in the x-dimension) of the exponent:
                   (e^(i z (Subscript[k, n])-((x^2)/((w^2)+((2 i z)/(Subscript[k, n])))))) ≈(e^(i  (Subscript[k, n])  r))  (e^(-((x^2)/((1+(((((2 z)/((Subscript[k, n]) (w^2)))))^2)))))) ≈  (e^(i  (Subscript[k, n])  r))  (e^(-((x^2)/((w^2) (((((2 z)/((Subscript[k, n]) (w^2)))))^2)(((((( (Subscript[k, n]) (w^2))/(2 z))))^2)+1))))) 
                                      ≈ (e^(i  (Subscript[k, n])  r))  (e^(-((((x/z))^2)/(((((2/((Subscript[k, n]) w))))^2)(((((( (Subscript[k, n]) (w^2))/(2 z))))^2)+1)))))  ≈ (e^(i  (Subscript[k, n])  r))  (e^(-(((Tan[θ])^2)/(((((2/((Subscript[k, n]) w))))^2)(((((( (Subscript[k, n]) (w^2))/(2 z))))^2)+1)))))
                                      ((≈e)^(i  (Subscript[k, n])  r))  (e^(-(((Tan[θ])^2)/((((((Subscript[λ, n])/(π w))))^2)(((((( (Subscript[k, n]) (w^2))/(2 z))))^2)+1)))))

What does the term  (e^(-(((Tan[θ])^2)/((((((Subscript[λ, n])/(π w))))^2)(((((( (Subscript[k, n]) (w^2))/(2 z))))^2)+1))))) tell us?  It is a Gaussian profile with farfield divergence angle θ, and for z-> ∞:
                                                      (e^(-(((Tan[θ])^2)/(((((Subscript[λ, n])/(π w))))^2)))) ~  (e^(-(((Tan[θ])^2)/((Tan[(Subscript[θ, 0])])^2))))              (7.3.2a) 
where:                                                 Tan[(Subscript[θ, 0])]=((Subscript[λ, n])/(π w))                   (7.3.2b)

In deed, this is the term that gives is the Gaussian profile along the propagation direction and with a divergence angle (Subscript[θ, 0])

λ=0.6328 ;n=1;w= 0.5 ; kn=2 π n/λ  ;
ParametricPlot3D[{r Cos[θ], r Sin[θ], (E^(-(((Tan[θ])^2)/(((((λ/(n π w))))^2)(((((( kn (w^2))/(2 r Sin[θ]))))^2)+1)))))},{r,0,5},{θ,0,2 π}, Mesh-> False, PlotPoints-> {100,15}, BoxRatios->{1,1,0.1}]

Now, we see that:

x  =>

Where (in which phenomena) do you expect to see this type of wave?

(We’ll see that at r=0, the approximation for the second plot is not correct because we don’t have z-> with correct limit).

But, now we can plot and verify a key basic property: divergence vs. beam waist and wavelength.

7.3.3 Beam divergence vs. beam waist

Recall Gaussian profile with farfield divergence angle θ for z-> ∞:
                                                      (e^(-(((Tan[θ])^2)/(((((Subscript[λ, n])/(π w))))^2)))) ~  (e^(-(((Tan[θ])^2)/((Tan[(Subscript[θ, 0])])^2)))) 
where:                                                 Tan[(Subscript[θ, 0])]=((Subscript[λ, n])/(π w))

Source code

Manipulate[
n=1; kn=2 π n/λ  ; T=ArcTan[ λ/(n π w)] ;
ParametricPlot3D[{r Cos[θ], r Sin[θ], (E^(-(((Tan[θ])^2)/((((λ/(n π w))))^2))))},{r,5,20},{θ,-4 T, 4 T}, Mesh-> False, PlotPoints-> {50,20}, PlotRange->{{0,20},{-20,20},{0,1}},BoxRatios->{1,2,0.2}]
,{{λ,0.5,Button["wavelength λ",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.036,1,0.888,1]]},0.5,1}
,{{w,1,Button["Beamwaist w",BaseStyle->{16, Bold,White,FontFamily-> "Calibri"},Background->Hue[0.344,1.,0.634,1]]},0.5,3}]

What do we see of the divergence as function of λ and w?

Demo Gaussian 1

We see how it moves vs. time with the last two plots.

Summary

Summary:
1. Any light packet (spatially as well as temporally) can be decomposed as a sum of many light components. (here, we will study spatial superposition first. Temporal superposition will be studied later).
2. The components can be pure harmonic plane waves (this is Fourier optics) or harmonic waves in other systems
3. The behavior of the beam is the net sum of those of the components
4. Most important to develop an  intuitive understand of the behavior of beam: "natural" description of the beam, not just any brute force solution