Oil, Soap, and Rainbows: Modeling Wave Interference with Thin Films
MIT Student
3.016 Fall 2017

Thin Films in Everyday Life

While we often don’t notice them, there are rainbows all around us as we go about our daily lives. They are hidden in the reflection of soap bubbles, the shine on a CD, and even in oily puddles on the street. This notebook explores the phenomenon of light wave interference and explains how it creates the colors we see on the surfaces of thin films, starting with the laws of light reflection and refraction, then perusing through some visual simulations of these fundamental principles.

GraphicsRow[{,,},ImageSize->1000]

Constructive Interference

When waves travel through space and hit an interface or a surface, some of the power of the wave reflects off the surface and the rest refracts, continuing through the new medium at a different angle. When the refracted wave hits another surface, part of it reflects back out of the medium and “adds” to the first reflected wave, or “interferes” with the first wave.

image=Import["https://upload.wikimedia.org/wikipedia/commons/thumb/5/51/Thin_film_interference.svg/300px-Thin_film_interference.svg.png"];
Show[image,ImageSize->500,Frame->True, FrameTicks-> False,FrameLabel->Style["Source: Wikipedia",10]]

Since the first and second reflected waves travel different paths, the second wave acquires a phase shift in comparison to the first wave. When it adds to the first wave, this phase shift may cause the two waves to completely cancel each other, resulting in destructive interference, or if the phases are perfectly aligned, the waves may add together to produce a resultant wave with twice the amplitude, constructive interference. Playing around with the simulation below, notice that the resultant wave disappears (destructive interference) when the phase shift is π radians, or 180°, and that it reaches maximum amplitude (constructive interference) when the phase shift is a multiple of 2π.

(* shows the interference between two sine waves with varying phase shift *)
Manipulate[
Module[{
firstSine=Sin[x],
secondSine=Sin[x-ϕ]},
Plot[{firstSine,secondSine,firstSine+secondSine},{x,-2Pi,2Pi},
ImageSize->Large,
PlotStyle->{{Dashed,Blue},{Dashed,Red},Purple},
PlotLabel->Style["Interference of Two Waves",20],
PlotRange->{-2,2}]
],
{{ϕ,0,"Phase Shift"},-2Pi,2Pi}
]

For light travelling from one medium (material 1) to another medium (material 2), i.e. from air to water, Snell’s law of refraction states that the angle at which light travels through two different mediums is proportional to the velocities of light through the two mediums ((Subscript[v, 1]) and (Subscript[v, 2])), which is inversely proportional to their refractive indices ((Subscript[n, 1]) and (Subscript[n, 2])):

((sin (Subscript[θ, 1]))/(sin (Subscript[θ, 2])))= ((Subscript[v, 1])/(Subscript[v, 2]))= ((Subscript[n, 2])/(Subscript[n, 1]))= ((Subscript[λ, 1])/(Subscript[λ, 2]))

This can be simplified to:

(Subscript[n, 1])sin (Subscript[θ, 1])=(Subscript[n, 2])sin (Subscript[θ, 2])

snell = (n1 Sin[θ1] == n2 Sin[θ2])

Snell’s law can be rearranged to calculate (Subscript[θ, 2]), the incident angle in the second medium, in terms of the refractive indices of the two materials and the incident angle of the light wave at the top surface.

calctheta2[theta1_,n1_,n2_]:= ArcSin[n1/n2 Sin[theta1]]

Geometrically, the path difference between the wave reflected off the top of the surface and the wave reflected off the bottom of the surface is given by (Subscript[n, 2])((OverBar[AB])+(OverBar[BC]))-(Subscript[n, 1])((OverBar[AD])), which is rewritten below using trigonometric substitutions for the line segment lengths.

geometricPathDiff=(pathDiff == 2 n2 (d/Cos[θ2])-n1 (2 d Tan[θ2]) ( Sin[θ1]))

By applying Snell’s Law to substitute for (Subscript[n, 1]), the path difference between the two reflected waves can be expressed in terms of the thin film medium only, only depending on (Subscript[n, 2]) and (Subscript[θ, 2]).

snellRules=Solve[snell,n1]

Simplify[Solve[geometricPathDiff /.{n1 -> n2 Sin[θ2]/Sin[θ1]},pathDiff]]

Thus, the interference between the two reflected waves will be constructive if this path difference is an integer multiple of the incident light’s wavelength, λ.

constructiveCondition = (m λ == 2 d n2 Cos[θ2])

When medium 2 has a refractive index that is greater than the index for medium 1 (that is, (Subscript[n, 2])>(Subscript[n, 1])), there is one caveat: the reflected wave has a 180° (or π radians) phase shift. In that case, the condition for constructive interference is that the path difference must be a half-integer multiple of the incident light’s wavelength.

constructiveConditionCaveat = ((m-1/2)λ==2d n2 Cos[θ2])

For a typical thin film, which has a higher refractive index than its surroundings (usually air), this is the condition for constructive interference between waves reflected off the top and bottom surfaces of the film.

λ (m - (1/2)) = 2 d (Subscript[n, 2]) cos (Subscript[θ, 2])

Great! Now, given the material and refractive index of the thin film and the outside medium, along with the incident angle of the light, we can calculate the wavelengths of visible light that will show the strongest signal, or constructive interference between the waves reflected off the top and bottom interfaces of the thin film. The following simulation shows a wave travelling through air, then reflecting off the top and bottom surfaces of a thin film, assuming that the refractive index of the film is greater than the refractive indices of the surrounding mediums (so that a 180° phase shift occurs at the top surface but not the bottom surface). The bold reflected wave represents the summation of the waves reflected off the top and bottom of the film, and grows largest when there is constructive interference and the two waves’ phases align.

(* Manipulate to change λ, θ1, d, and nFilm, showing the wave path and its interference *)
Manipulate[
Module[
{background=Graphics[
{Lighter[Blue],HalfPlane[{{0,0},{1,0}},{0,-1}], (* blue background *)
Lighter[LightBlue],Rectangle[{-5000,-d},{5000,0}],(* light blue thin film *)
Dashed,Thick,Black,Line[{{0,0},{0,5000}}], (* vertical line in the middle to compare incident wave angle to *)
Arrowheads[{-0.02,0.02}],Arrow[{{-4500,0},{-4500,-d}}]}, (* arrow showing thickness d of the film*)
Frame->True,
FrameTicks->False,
ImageSize->Large,
PlotRange->{{-5000,5000},{-5000,5000}},
Epilog->{Text[Style["(Subscript[n, air])=1",14],{-4500,4500}]}],
incidentWave=ParametricPlot[RotationMatrix[θ1 Degree].{100Sin[2Pi t/λ],t},{t,0,5000}, (* incident wave at upper surface *)
Axes->False,
PlotStyle->ColorData["VisibleSpectrum"][λ]],
upperReflectedWave=ParametricPlot[RotationMatrix[-θ1 Degree].{100Cos[2Pi t/λ],t},{t,0,5000}, (* reflected wave at upper surface, with a phase shift *)
Axes->False,
PlotStyle->{Opacity[0.3],ColorData["VisibleSpectrum"][λ]}],
refractedAngle= ArcSin[Sin[θ1 Degree]/nFilm], (* calculate the refracted wave's angle using Snell's law *)
refractedWavelength,
refractedWavePath,
waveColor=ColorData["VisibleSpectrum"][λ]},

refractedWavelength = λ*Sin[refractedAngle]/Sin[θ1 Degree];
refractedWavePath = d ArcCos[refractedAngle];
Show[
background,
incidentWave,
upperReflectedWave,
ParametricPlot[RotationMatrix[refractedAngle].{100Cos[2Pi t/refractedWavelength],t},{t,-refractedWavePath,0}, (* refracted wave before hitting the bottom surface *)
Axes->False,
PlotStyle->{Opacity[0.5],waveColor}],
ParametricPlot[RotationMatrix[-refractedAngle].{100Cos[2Pi t/refractedWavelength],t}+{3 refractedWavePath Sin[refractedAngle],refractedWavePath Cos[refractedAngle]},{t,-2 refractedWavePath,-refractedWavePath}, (* refracted wave reflecting off the bottom surface, without a phase shift *)
Axes->False,
PlotStyle->{Opacity[0.5],waveColor}],
ParametricPlot[RotationMatrix[-θ1 Degree].{100Cos[2Pi t/λ],t}+{2 refractedWavePath Sin[refractedAngle],0},{t,0,5000}, (* wave re-refracted back into the air, which interferes with the wave reflected off the top surface *)
Axes->False,
PlotStyle->{Opacity[0.3],waveColor}],
ParametricPlot[RotationMatrix[-θ1 Degree].{100Cos[2Pi(t-2 refractedWavePath Sin[refractedAngle]Sin[θ1])/λ]+100Sin[2Pi t/λ],t}+{refractedWavePath Sin[refractedAngle],0},{t,0,5000}, (* sum of the two reflected waves interfering, showing constructive or destructive interference *)
Axes->False,
PlotStyle->waveColor],
PerformanceGoal->{"Speed","Quality"}
]
],
{{λ,380},380,750},
{{θ1,45},0,90},
{{nFilm,1.5},1,3},
{{d,500},200,2000},
ControlPlacement->Left
]

Wavelength and Color

The visible light spectrum includes wavelengths ranging from 380 nm to 750 nm. This means that the light that we see making rainbow swirly patterns on the surface of soap bubbles, or in the oil sheen on water, has a minimum wavelength of 380 nm and a maximum wavelength of 750 nm.

Given a thin film with thickness d, refractive index (Subscript[n, film]), and an incident light angle of (Subscript[θ, i]), the particular colors of light that we can see reflecting off the surface of the film must have wavelengths that satisfy the constructive interference condition. Air has a refractive index of about 1, and a thin film of soap has a refractive index larger than 1. In this case, the constructive interference condition includes the phase-shift caveat. For example, a film with thickness 200 nm and refractive index 1.5, hit with white light at an angle of 60°, will lead to constructive interference for light with a wavelength of 600 nm.

(* takes in the film thickness d, refractive index n, and incident light angle θ, and return the visible wavelengths of light that satisfy the constructive interference condition for a soap bubble *)
reflectedWavelengths[d_,n_,theta_]:=Select[2 d n Cos[theta]/#&/@Range[Floor[2 d n Cos[theta]]],380<#<750&]

(* same as reflectedWavelengths, but for cases where (Subscript[n, 2])>(Subscript[n, 1]) *)
reflectedWavelengthsCaveat[d_,n_,theta_]:=Select[2 d n Cos[theta]/#&/@Range[0.5,Floor[2 d n Cos[theta]]+0.5,1],380<#<750&]

reflectedWavelengthsCaveat[200,1.5,60 Degree]

The human eye has light receptors, or cones, that actually only perceive three colors from the visible spectrum: red, green, and blue. The entire visible spectrum people see is really a superposition of the waves of these three colors, which mix in different ratios to provide all the other colors of the rainbow. Red has a wavelength of about 700 nm, green 546 nm, and blue 436 nm. 

The strength of each of these colors, or the intensity, that people perceive is proportional to the square of the wave’s amplitude. Intensity is also a representation of the wave’s energy density, and since the amplitude represents the energy (E) of the wave, the intensity (I) of a wave can be plotted as the wave’s function squared. The more intense (larger amplitude) a wave of a certain wavelength, the more that color dominates our perception.

I ∝|E((|)^2)

The interference and intensities of the red, green, and blue waves when reflected off a thin film with different thicknesses and incident light angles is visualized below.

interferingWave[d_,n_,theta_,wavelength_]:=Sin[2Pi x/wavelength]+Sin[2Pi( x-2 d n Cos[theta])/wavelength] (* sum of the waves reflected off the top and bottom surfaces of a film, given a wavelength of light, film thickness, refractive index, and incident light angle *)

Manipulate[
Module[{redWave=interferingWave[d,n,theta Degree,700],
greenWave=interferingWave[d,n,theta Degree,546],
blueWave=interferingWave[d,n,theta Degree,436]},
GraphicsRow[
{Plot[{redWave,
greenWave,
blueWave},
{x,0,200Pi},
PlotStyle->{Red,Green,Blue},
ImageSize->Large,
PlotLabel->Style["Interference of Red, Green and Blue Waves Off a Thin Film",18],
PlotRange->{-2,2},
AxesLabel->{"nm",Null}],
Plot[{(redWave)^2,(* red wave intensity *)
(greenWave)^2,(* green wave intensity *)
(blueWave)^2},(* blue wave intensity*)
{x,0,200Pi},
PlotStyle->{Red,Green,Blue},
ImageSize->Large,
PlotLabel->Style["Intensity of Red, Green and Blue Waves Off a Thin Film",18],
PlotRange->{0,4},
AxesLabel->{"nm","Intensity"}]}]
],
{{d,200,"Film Thickness (nm)"},100,1000},
{{n,1.5, "Refractive Index"},1,3},
{{theta,0,"Incident Light Angle"},0, 180}
]

The total color that we see reflected off a point on the film surface is a combination of these three waves proportional to their intensities.

(* finds the intensities of the red, green, and blue wavelengths, then divides by 4 (the max intensity possible for two added sine waves) to make all values between 0 and 1, to feed into the RGBColor function. Note: also takes into account the 180° phase shift from hitting the surface *)
reflectedColorIntensities[d_,n_,theta_]:=FindMaxValue[(Sin[2Pi x/#]+Sin[2Pi (x-2 d n Cos[theta]-Pi)/#])^2/4,x,PrecisionGoal->5,AccuracyGoal->5]&/@{700,546,436}

reflectedColor[d_,n_,theta_]:=RGBColor[reflectedColorIntensities[d,n,theta]] (* converts R, G, and B intensities into a color *)

Assuming that the thin film has a greater refractive index than the air, we can plot the colors of the wavelengths that constructively interfere when reflected off a thin film as a function of the film thickness, with a constant film refractive index of 1.5 and incident angle of 45°. Likewise, we can see how the reflected wavelength colors vary with the incident light wave angle when thickness and refractive index are constant at 200 nm and 1.5, respectively.

(* list of colors reflected off a thin film for film thicknesses 270 to 1000 *)
colorFunctionThickness=reflectedColor[#,1.5,Pi/4]&/@Range[270,1000];

colorFunctionAngle=reflectedColor[200,1.5,# Degree]&/@Range[90]; (* list of colors reflected off a thin film of thickness 200 through angles 0 to 90 *)

GraphicsRow[
{Plot[1,{x,100,800},
Filling->Axis,
ColorFunction-> (Blend[colorFunctionThickness,#1]&),
PlotRange->{0,1},
ImageSize->Large,
PlotLabel->Style["Reflected Colors Off a Thin Film With Varying Thickness", Black, 20],
Axes->{True,False},
Frame->{True,False},
FrameLabel->{Style["Film Thickness (nm)",14]}
],
Plot[1,{x,0,90},
Filling->Axis,
ColorFunction-> (Blend[colorFunctionAngle,#1]&),
PlotRange->{0,1},
ImageSize->Large,
PlotLabel->Style["Reflected Colors Off a Thin Film With Varying Incidence Angle", Black, 20],
Axes->{True,False},
Frame->{True,False},
FrameLabel->{Style["Angle of Incident Light Wave (°)",14]}
]}
]

Oil Film

Let’s envision a flat thin film of oil floating on a flat surface of water.

The refractive indices of oil, water, and air are respectively about 1.5, 1.33, and 1. As light reflects off the top of the oil, there is a 180° phase shift since (Subscript[n, oil])>(Subscript[n, air]), but there isn’t a phase shift when the wave reflects off the interface between the oil and the water because (Subscript[n, water])<(Subscript[n, oil]).

oilIndex = 1.5;
waterIndex = 1.33;
airIndex = 1;

(* calculates the color of wavelength of light reflected at a given point on the flat film surface in relation to the point source of light, returning the color of the reflected wave *)
waveColorPointLight[filmpoint_,lightpoint_,d_]:=Module[{theta=ArcTan[EuclideanDistance[{0,0},filmpoint[[1;;2]]]/(lightpoint[[3]]+5)]},
reflectedColor[d,1.5,theta]
]

The interactive graphic below shows how an oil film with varying thickness reflects light of a point source at varying heights above the oil’s surface. Notice that increased thicknesses of the oil layer lead to more muted colors, but a wider range of angles for which specific colors show.

Manipulate[
Module[
{background=Graphics3D[
{Opacity[0.1],Yellow,Cuboid[{-100,-100,-5-filmThickness/10},{100,100,-5}] (* oil film *),
Opacity[0.4],Lighter[Blue],Cuboid[{-100,-100,-5-filmThickness/10},{100,100,-100}] (* water *)},
PlotRange->{{-101,101},{-101,101},{-100,100}},
Lighting->{"Ambient",White},
ImageSize->Large],
coloring=Function[{x,y,z},{waveColorPointLight[{x,y,z},{0,0,lightHeight+5},filmThickness]}],
pointLight =Graphics3D[{PointSize[Large],Point[{0,0,lightHeight}]}]
},
Show[
background,
pointLight,
Plot3D[-5,{x,-100,100},{y,-100,100},
PlotStyle->Opacity[0.2],
ColorFunction-> coloring,
ColorFunctionScaling->False,
Mesh->False,
PerformanceGoal->{"Speed","Quality"}],
PerformanceGoal-> {"Speed","Quality"}
]
],
{{lightHeight,0},0,100},
{{filmThickness,100},100,1000}
]

In reality, the surface of a puddle isn’t really completely flat. What if the surface had waves; how would the light reflecting off the oil look then? The interactive graphic below shows thin film interference off oil on rippling water; try moving around the position of the point source light and changing the thickness of the oil layer to see how the reflected colors shift around.

(* calculates the wavelength of light reflected at a given point on the ripple film surface in relation to the point source of light, returning the color of the reflected wave *)
waveColorRipple[filmpoint_,lightpoint_,thickness_]:=Module[{theta=VectorAngle[Grad[c-10Sin[a/20]Cos[b/20],{a,b,c}],lightpoint-filmpoint]/.{a->filmpoint[[1]],b->filmpoint[[2]],c->filmpoint[[3]]}},
reflectedColor[thickness,1.5,theta]
]

Manipulate[
Module[
{background=
Show[
Plot3D[10Sin[x/20]Cos[y/20]-thickness/10,{x,-100,100},{y,-100,100}, (* wavy water background *)
PlotRange->{-100,200},
PlotStyle->{Lighter[Blue],Opacity[0.4]},
Lighting->{ "Ambient",White},
Filling->Bottom,
FillingStyle->{Lighter[Blue],Opacity[0.4]},
ImageSize->Large,
Mesh->False,
Axes->False],
Plot3D[10Sin[x/20]Cos[y/20],{x,-100,100},{y,-100,100}, (* oil film *)
PlotRange->{-100,100},
Mesh->False,
PlotStyle->{Yellow,Opacity[0.2]},
Lighting->{ "Ambient",White},
Filling->-thickness/10-5,
FillingStyle->{LightYellow,Opacity[0.2]}],
PerformanceGoal->{"Speed","Quality"}],
coloring=Function[{x,y,z},{waveColorRipple[{x,y,z},{a,b,c},thickness]}],
function= 10 Sin[x/20]Cos[y/20],
pointLight =Graphics3D[{PointSize[Large],Point[{a,b,c}]}] },

Show[
background,
pointLight,
Plot3D[function,{x,-100,100},{y,-100,100},
PlotStyle->Opacity[0.15],
ColorFunction-> coloring,
ColorFunctionScaling->False,
Mesh->False,
PerformanceGoal->{"Speed","Quality"}]
]
],
"Light Position",
{{a,0,"x"},-100,100},
{{b,0,"y"},-100,100},
{{c,100,"z"},0,200},
Delimiter,
{{thickness,100},100,1000},
ControlPlacement->Left
]

Soap Bubble

For a spherical bubble, the refractive index of soap varies with the recipe, but it’s always greater than the refractive index of air. This means that there is also a phase shift of 180° when light hits the outer surface surface of the bubble, but not when it reflects off the inner surface since the (Subscript[n, air])<(Subscript[n, soap]). Assuming the bubble has an even thickness of 100 nm all around and (Subscript[n, soap])≈1.4, the iridescence patterns on its surface should look pretty similar to the pattern we saw on a flat oil film, just wrapped around a ball shape.

The interactive graphic below shows the effect of multiple point source lights on a soap bubble. Note that in real life, the light sources around a soap bubble usually aren’t just point sources, and the bubble probably isn’t an even thickness throughout. The color patterns on bubble surfaces are usually far more complex than what we can see in this simulation!

(* gives the color of the reflected wavelength at a given point on a sphere in relation to a point source of light *)
waveColorSphere[spherePoint_,lightPoint_]:=Module[{theta=VectorAngle[spherePoint,lightPoint-spherePoint]},
reflectedColor[100,1.4,theta]]

(* gives the color of the reflected wavelength at a given point on a sphere in relation to 2 point sources of light *)
waveColorSphere2Lights[spherePoint_,light1_,light2_]:=Module[{theta1=VectorAngle[spherePoint,light1-spherePoint],theta2=VectorAngle[spherePoint,light2-spherePoint]},
Blend[{reflectedColor[100,1.4,theta1],reflectedColor[100,1.4,theta2]}]]

(* gives the color of the reflected wavelength at a given point on a sphere in relation to 3 point sources of light *)
waveColorSphere3Lights[spherePoint_,light1_,light2_,light3_]:=Module[{theta1=VectorAngle[spherePoint,light1-spherePoint],theta2=VectorAngle[spherePoint,light2-spherePoint],theta3=VectorAngle[spherePoint,light3-spherePoint]},
Blend[{reflectedColor[100,1.4,theta1],reflectedColor[100,1.4,theta2],reflectedColor[100,1.4,theta3]}]]

Manipulate[
Module[{sphere=10SpherePoints[10000],
coloring1=Function[{x,y,z},waveColorSphere[{x,y,z},{x1,y1,z1}]], (* color function if 1 light point is on *)
coloring2=Function[{x,y,z},waveColorSphere2Lights[{x,y,z},{x1,y1,z1},{x2,y2,z2}]],(* color function if 2 light points are on *)
coloring3=Function[{x,y,z},waveColorSphere3Lights[{x,y,z},{x1,y1,z1},{x2,y2,z2},{x3,y3,z3}]]},(* color function if 3 light points are on *)

ListSurfacePlot3D[sphere,
ColorFunctionScaling->False,
ImageSize->Large,
PlotStyle->{Opacity[0.3],White},
Mesh->None,
Lighting->{"Ambient",White},
ColorFunction->Which[numberOfLights=="1 light",coloring1,numberOfLights=="2 lights",coloring2,numberOfLights=="3 lights",coloring3],
PerformanceGoal->"Speed"]],

{numberOfLights,{"1 light","2 lights","3 lights"}},
Delimiter,
"Light 1 Position",
{{x1,0},-100,100},
{{y1,0},-100,100},
{{z1,50},-100,100},
Delimiter,
"Light 2 Position",
{{x2,0},-100,100},
{{y2,0},-100,100},
{{z2,50},-100,100},
Delimiter,
"Light 3 Position",
{{x3,0},-100,100},
{{y3,0},-100,100},
{{z3,50},-100,100},
ControlPlacement->Left]

References and Thanks

Special thanks to the Professor Carter, Dr. Keane, and the amazing 3.016 TAs Emma, Felipe, and Burhan for their guidance and help on this project.

References:

https://light-measurement.com/perception-of-color/
https://en.wikipedia.org/wiki/Thin-film_interference
https://en.wikipedia.org/wiki/Snell%27s_law
https://en.wikipedia.org/wiki/Intensity_(physics)