Interference is a fundamental physical phenomenon that demonstrates the superposition principle for linear systems. This principle, which states that the net response to multiple stimuli is the sum of the individual responses, is central to our understanding of wave physics. Interference appears across many domains of physics: in optics where it enables high-precision measurements and holography, in quantum mechanics where it reveals the wave nature of matter, and in acoustics where it forms the basis for noise cancellation technology. The ability of waves to interfere constructively (amplifying each other) or destructively (canceling each other) has profound practical applications, from the anti-reflective coatings on optical elements to the operational principles of interferometric gravitational wave detectors like LIGO. Understanding interference is therefore not just of theoretical interest but crucial for modern technology and experimental physics.
When two wave solutions \(U_1(\mathbf{r})\) and \(U_2(\mathbf{r})\) combine, their superposition gives:
\[ U(\mathbf{r})=U_1(\mathbf{r})+U_2(\mathbf{r}) \]
The resulting intensity is:
\[\begin{eqnarray} I &= &|U|^2\\ &= &|U_1+U_2|^2\\ &= &|U_1|^2+|U_2|^2+U^{*}_1 U_2 + U_1 U^{*}_2 \end{eqnarray}\]
The individual wave intensities are given by \(I_1=|U_1|^2\) and \(I_2=|U_2|^2\). Using this, we can express each complex wave amplitude in polar form, separating its magnitude (related to intensity) and phase:
\[ U_1=\sqrt{I_1}e^{i\phi_1} \] \[ U_2=\sqrt{I_2}e^{i\phi_2} \]
Substituting these expressions back into our interference equation and performing the algebra, the total intensity becomes:
\[ I=I_1+I_2+2\sqrt{I_1 I_2}\cos(\Delta \phi) \]
where \(\Delta \phi=\phi_2-\phi_1\) is the phase difference between the waves. This equation is known as the interference formula and contains three terms:
A particularly important special case occurs when the interfering waves have equal intensities (\(I_1=I_2=I_0\)). The equation then simplifies to:
\[ I=2I_0(1+\cos(\Delta \phi))=4I_0\cos^2\left(\frac{\Delta \phi}{2}\right) \]
This last form clearly shows that:
Occurs when \(\Delta \phi=2\pi m\) (where \(m\) is an integer), resulting in \(I=4I_0\)
x=np.linspace(0,2,1000)
wavelength=0.532
k=2*np.pi/0.532
y1=np.cos(k*x)
fig,[ax1,ax2,ax3]=plt.subplots(3,1,figsize=get_size(10,8))
ax1.plot(x/wavelength,y1,label='Wave 1')
ax2.plot(x/wavelength,y1,label='Wave 1')
ax3.plot(x/wavelength,2*y1,label='Wave 1')
ax3.set_xlabel(r"distance [$\lambda$]")
ax1.set_ylabel(r"$U_1$")
ax2.set_ylabel(r"$U_2$")
ax3.set_ylabel(r"$U_1+U_2$")
plt.tight_layout()
plt.show()
Occurs when \(\Delta \phi=(2m-1)\pi\) (where \(m\) is an integer), resulting in \(I=0\)
x=np.linspace(0,2,1000)
wavelength=0.532
k=2*np.pi/0.532
y1=np.cos(k*x)
y2=np.cos(k*x+np.pi)
fig,[ax1,ax2,ax3]=plt.subplots(3,1,figsize=get_size(10,8))
ax1.plot(x/wavelength,y1,label='Wave 1')
ax2.plot(x/wavelength,y2,label='Wave 1')
ax3.plot(x/wavelength,y1+y2,label='Wave 1')
ax3.set_xlabel(r"distance [$\lambda$]")
ax1.set_ylabel(r"$U_1$")
ax2.set_ylabel(r"$U_2$")
ax3.set_ylabel(r"$U_1+U_2$")
ax3.set_ylim(-1,1)
plt.tight_layout()
plt.show()
The phase difference \(\Delta \phi\) can be related to the path difference \(\Delta s\) between the two waves. For two waves with the same frequency \(\omega\), we can write their complete phase expressions as:
\[\phi_1(\mathbf{r},t) = \mathbf{k}_1\cdot\mathbf{r} - \omega t + \phi_{01}\] \[\phi_2(\mathbf{r},t) = \mathbf{k}_2\cdot\mathbf{r} - \omega t + \phi_{02}\]
where:
The instantaneous phase difference is then:
\[ \Delta\phi(\mathbf{r},t) = \phi_2(\mathbf{r},t) - \phi_1(\mathbf{r},t) = (\mathbf{k}_2-\mathbf{k}_1)\cdot\mathbf{r} + (\phi_{02}-\phi_{01}) \]
For stationary interference patterns, we typically observe the time-independent phase difference. When the waves travel along similar paths (same direction), this reduces to:
\[\Delta\phi = k\Delta s + \Delta\phi_0\]
where \(\Delta s\) is the path difference and \(\Delta\phi_0\) is any initial phase difference between the sources.
A path difference \(\Delta s\) corresponds to a phase difference \(k\Delta s=2\pi\Delta s/\lambda\). Path differences of integer multiples of \(\lambda\) result in phase differences of integer multiples of \(2\pi\).
def plane_wave(k,omega,r,t):
return(np.exp(1j*(np.dot(k,r)-omega*t)))
wavelength=532e-9
k0=2*np.pi/wavelength
c=299792458
omega0=k0*c
vec=np.array([0.0,0.,1.])
vec1=np.array([1.0,0.,1.])
vec=vec/np.sqrt(np.dot(vec,vec))
vec1=vec1/np.sqrt(np.dot(vec1,vec1))
k=k0*vec
k1=k0*vec1
x=np.linspace(-2.5e-6,2.5e-6,300)
z=np.linspace(0,5e-6,300)
X,Z=np.meshgrid(x,z)
r=np.array([X,0,Z],dtype=object)
fig,ax=plt.subplots(2,2,figsize=get_size(10,10))
field=plane_wave(k,omega0,r,0)
field1=plane_wave(k1,omega0,r,0)
extent = np.min(z)*1e6, np.max(z)*1e6,np.min(x)*1e6, np.max(x)*1e6
ax[0,0].imshow(np.real(field.transpose()),extent=extent,vmin=-1,vmax=1,cmap='seismic')
ax[0,0].set_title('wave 1')
ax[0,1].imshow(np.real(field1.transpose()),extent=extent,vmin=-1,vmax=1,cmap='seismic')
ax[0,1].set_title('wave 2')
ax[1,0].imshow(np.real(field.transpose()+field1.transpose()),extent=extent,vmin=-1,vmax=1,cmap='seismic')
ax[1,0].set_title('sum')
ax[1,1].imshow(np.abs(field.transpose()+field1.transpose())**2,extent=extent,cmap='gray')
ax[1,1].set_title('intensity')
ax[1,1].set_xlabel('z-position [µm]')
ax[1,0].set_xlabel('z-position [µm]')
ax[1,0].set_ylabel('x-position [µm]')
ax[0,0].set_ylabel('x-position [µm]')
plt.tight_layout()
plt.show()
def spherical_wave(k,omega,r,r0,t):
k=np.linalg.norm(k)
d=np.linalg.norm(r-r0)
return( np.exp(1j*(k*d-omega*t))/d)
x=np.linspace(-5e-6,5e-6,300)
z=np.linspace(-5e-6,5e-6,300)
X,Z=np.meshgrid(x,z)
r=np.array([X,0,Z],dtype=object)
wavelength=532e-9
k0=2*np.pi/wavelength
c=299792458
omega0=k0*c
k=k0*np.array([0,1.0,0])
r0=np.array([0,2e-6,0])
field=spherical_wave(k,omega0,r,r0,0)
field1=plane_wave(k,omega0,r,0)
extent = np.min(z)*1e6, np.max(z)*1e6,np.min(x)*1e6, np.max(x)*1e6
fig,ax=plt.subplots(2,2,figsize=get_size(12,12))
ax[0,0].imshow(np.real(field.transpose()+0*field1.transpose()),extent=extent,cmap='seismic')
ax[0,0].set_title('Spherical wave')
ax[0,1].imshow(np.real(0*field.transpose()+field1.transpose()),extent=extent,cmap='seismic')
ax[0,1].set_title('Plane wave')
ax[1,0].imshow(np.real(0.00001*field.transpose()+field1.transpose()),extent=extent,cmap='seismic')
ax[1,0].set_title('Sum')
ax[1,1].imshow(np.abs(0.00001*field.transpose()+field1.transpose())**2,extent=extent,cmap='gray')
ax[1,1].set_title('Intensity')
ax[1,0].set_xlabel('z [µm]')
ax[1,1].set_xlabel('z [µm]')
ax[1,0].set_ylabel('x [µm]')
ax[0,0].set_ylabel('x [µm]')
plt.tight_layout()
plt.show()
The interference of the spherical and the plane wave (also the one of the two plane waves) give also an interesting result. The intensity resembles to be a snapshot of the shape of the wavefronts of the spherical wave. We can therefore measure the wavefronts of the spherical wave by interfering it with a plane wave. This is also the basic principle behind holography. There we use a reference wave to interfere with the wave that we want to measure. The interference pattern is recorded and can be used to reconstruct the wavefronts of the wave.
A super nice website to try out interference interactively is here.
In the earlier consideration we obtained a general description for the phase difference between two waves. TIt is given by and contains the pathlength difference \(\Delta s\) and some intrinsic phase \(\Delta\phi_0\) that could be part of the wave generation process.
\[\Delta\phi = k\Delta s + \Delta\phi_0\]
To observe stationary interference, it is important that these two quantities are also stationary, i.e. the phase relation between the two waves is stationary. This relation between the phase of two waves is called coherence and was assumed in all the examples before.
The above image shows the timetrace of the amplitude of two wave with slightly different frequency. Due to the frequency, the waves run out of phase and have acquired a phase different of \(\pi\) after \(40\) fs.
The temporal coherence of two waves is now defined by the time it takes for the two waves to obtain a phase difference of \(2\pi\). The phase difference between two wave of frequency \(\nu_1\) and \(\nu_2\) is given by
\[ \Delta \phi = 2\pi (\nu_2-\nu_1)(t-t_0) \]
Here \(t_0\) refers to the time, when thw two waves were perfectly in sync. Lets assume that the two frequencies are seperarated from a central frequency \(\nu_0\) such that
\[ \nu_1=\nu_0-\Delta \nu/2 \] \[ \nu_2=\nu_0+\Delta \nu/2 \]
Inserting this into the first equation yields
\[ \Delta \phi = 2\pi \Delta \nu \Delta t \]
with \(\Delta t=t-t_0\). We can now define the coherence time as the time interval over which the phase shift \(\Delta \phi\) grows to \(2\pi\), i.e. \(\Delta \phi=2\pi\). The coherence time is thus
\[ \tau_{c}=\Delta t =\frac{1}{\Delta \nu} \]
Thus the temporal coherence and the frequency distribution of the light are intrisincly connected. Monochromatic light has \(\Delta nu=0\) and thus the coherence time is infinitely long. Light with a wide spectrum (white light for example) therefore has and extremly short coherence time.
The coherence time is also connected to a coherence length. The coherence length \(L_c\) is given by the distance light travels within the coherence time \(\tau_c\), i.e.
\[ L_c=c\tau_c \]
Two waves are called coherent, if they exihibit a fixed phase relation in space or time relation over time. It measures their ability to interfer. The main types of coherence are
Coherence is a property of the light source and is connected to the frequency distribution of the light. Sources can be:
While the above definition provides an intuitive picture based on frequency spread, we can describe coherence more rigorously using correlation functions. These functions measure how well a wave maintains its phase relationships:
In real physical systems, perfect coherence (constant phase relationship) between waves is rare. Partial coherence describes the degree to which waves maintain a consistent phase relationship over time and space. We can characterize this using correlation functions:
\[g^{(1)}(\tau) = \frac{\langle U(t)U^*(t+\tau)\rangle}{\sqrt{\langle|U(t)|^2\rangle\langle|U(t+\tau)|^2\rangle}}\]
where:
\[g^{(1)}(\mathbf{r}_1,\mathbf{r}_2) = \frac{\langle U(\mathbf{r}_1)U^*(\mathbf{r}_2)\rangle}{\sqrt{\langle|U(\mathbf{r}_1)|^2\rangle\langle|U(\mathbf{r}_2)|^2\rangle}}\]
The obtained correlation functions can be used to calculate the coherence time and length and have the following properties:
A finite coherence time and length is leads to partial coherence affects interference visibility through:
omega0 = 2.0
delta_omega = 0.05 # frequency difference
tau_c = np.pi/delta_omega # coherence time (corrected)
beat_period = 2*np.pi/delta_omega # time for full beat cycle
t = np.linspace(0, 1000, 10000)
tau = np.linspace(0, 500, 200)
def generate_waves(t):
wave1 = np.exp(1j * omega0 * t)
wave2 = np.exp(1j * (omega0 + delta_omega) * t)
return wave1, wave2
def calc_correlation(wave, tau):
g = np.zeros(len(tau), dtype=complex)
N = len(wave)
for i, dt in enumerate(tau):
shift = int(dt * 10)
if shift >= N:
g[i] = 0
else:
g[i] = np.mean(wave[:(N-shift)] * np.conj(wave[shift:]))
return g / np.abs(g[0])
wave1, wave2 = generate_waves(t)
wave_total = wave1 + wave2
# Calculate correlation
g = calc_correlation(wave_total, tau)
# Plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=get_size(10, 8))
# Plot waves
ax1.plot(t[:500], np.real(wave1[:500]), label='Wave 1', alpha=0.7)
ax1.plot(t[:500], np.real(wave2[:500]), label='Wave 2', alpha=0.7)
ax1.plot(t[:500], np.real(wave_total[:500]), 'k', label='Sum', alpha=0.7,lw=0.5)
ax1.set_title('wave superposition')
ax1.set_xlabel('time')
ax1.set_ylabel('amplitude')
ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
# Plot correlation
ax2.plot(tau, np.abs(g))
ax2.axvline(x=tau_c, color='r', linestyle='--', label=r'$\tau_c$ ')
ax2.axvline(x=beat_period, color='g', linestyle=':', label=f'Beat period')
ax2.set_title('|g⁽¹⁾(τ)|')
ax2.set_xlabel('τ')
ax2.set_ylabel('|g⁽¹⁾(τ)|')
ax2.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.subplots_adjust(right=0.85)
plt.show()
Besides different frequencies the coherence time can also be affected by phase jumps. The following example shows two waves with the same frequency but multiple phase jumps. The temporal correlation function shows the decoherence due to the phase jumps.
import numpy as np
import matplotlib.pyplot as plt
omega0 = 1.0 # same frequency for both waves
tau = np.linspace(0, 500, 200)
t = np.linspace(0, 1000, 10000)
def generate_waves_with_jumps(t, n_jumps=10):
# Create two identical waves
wave1 = np.exp(1j * omega0 * t)
wave2 = np.exp(1j * omega0 * t) # same frequency
# Create regularly spaced jumps within first 500 time units
jump_positions = np.linspace(0, 500, n_jumps+1)[:-1] # exclude last point
jump_indices = [int(pos * len(t)/t[-1]) for pos in jump_positions]
phase_shifts = np.random.uniform(0, 2*np.pi, n_jumps)
# Apply phase shifts to wave2
wave2_with_jumps = wave2.copy()
current_phase = 0
for i in range(n_jumps):
start_idx = jump_indices[i]
if i < n_jumps-1:
end_idx = jump_indices[i+1]
else:
end_idx = len(t)
current_phase += phase_shifts[i]
wave2_with_jumps[start_idx:end_idx] *= np.exp(1j * current_phase)
return wave1, wave2_with_jumps, jump_positions
def calc_correlation(wave, tau):
g = np.zeros(len(tau), dtype=complex)
N = len(wave)
for i, dt in enumerate(tau):
shift = int(dt * 10)
if shift >= N:
g[i] = 0
else:
g[i] = np.mean(wave[:(N-shift)] * np.conj(wave[shift:]))
return g / np.abs(g[0])
# Generate waves with 30 jumps
wave1, wave2, jump_positions = generate_waves_with_jumps(t, n_jumps=30)
wave_total = wave1 + wave2
g = calc_correlation(wave_total, tau)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=get_size(10, 8))
ax1.plot(t[:2000], np.real(wave1[:2000]), label='Wave 1', alpha=0.9)
ax1.plot(t[:2000], np.real(wave2[:2000]), label='Wave 2', alpha=0.9)
ax1.plot(t[:2000], np.real(wave_total[:2000]), 'k-', label='Sum', lw=0.5)
ax1.set_xlim(0, 200)
# Add vertical lines for phase jumps in wave plot
for pos in jump_positions:
ax1.axvline(x=pos, color='r', linestyle='--', alpha=0.3)
ax1.set_title('Superposition with Multiple Phase Jumps')
ax1.set_xlabel('time')
ax1.set_ylabel('amplitude')
ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
# Plot correlation
ax2.plot(tau, np.abs(g))
ax2.set_title('|g⁽¹⁾(τ)|')
ax2.set_xlabel('τ')
ax2.set_ylabel('|g⁽¹⁾(τ)|')
ax2.set_xlim(0, 200)
ax2.set_ylim(0, 1)
# Adjust layout
plt.tight_layout()
plt.subplots_adjust(right=0.85)
plt.show()
Multiple Wave Interference
So far we looked at the interference of two waves, which was a simplification as I mentioned already earlier. Commonly there will be a multitude of partial waves contribute to the oberved intereference. This is what we would like to have a look at now. We will do that in a quite general fashion, as the resulting formulas will appear several times again for different problems.
Nevertheless we will make a difference between
Especially the latter is often occuring, if we have multiple reflections and each reflection is only a fraction of the incident amplitude.
In the case of constant amplitude (for example realized by a grating, which we talk about later), the total wave amplitude is given according to the picture below by
\[ U=U_1+U_2+U_1+U_3+\ldots+U_M \]
where we sum the amplitude over \(M\) partial waves. Between the neighboring waves (e.g. \(U_1\) and \(U_2\)), we will assume a phase difference (because of a path length difference for example), which we denote as \(\Delta \phi\).
The amplitude of the p-th wave is then given by
\[ U_p=\sqrt{I_0}e^{i(p-1)\Delta \phi} \]
with the index \(p\) being an interger \(p=1,2,\ldots,M\), \(h=e^{i\Delta \phi}\) and \(\sqrt{I_0}\) as the amplitude of each individual wave. The total amplitude \(U\) can be then expressed as
\[ U=\sqrt{I_0}\left (1+h+h^2+\ldots +h^{M-1}\right) \]
which is a geometric sum. We can apply the sum formula for geometric sums to obtain
\[ U=\sqrt{I_0}\frac{1-h^M}{1-h}=\sqrt{I_0}\frac{1-e^{iM\Delta \phi}}{1-e^{i\Delta \phi}} \]
We now have to calculate the intensity of the total amplitude
\[ I=|U|^2=I_{0}\left | \frac{e^{-iM\Delta \phi/2}-e^{iM\Delta \phi/2}}{e^{-i\Delta \phi/2}-e^{i\Delta \phi/2}}\right |^2 \]
which we can further simplify to give
\[ I=I_{0}\frac{\sin^2(M\Delta \phi/2)}{\sin^2(\Delta \phi/2)} \]
# Parameters
M = 6 # number of phasors
phi = np.pi/8 # example phase difference between successive phasors
def plot_angle(ax, pos, angle, length=0.95, acol="C0", **kwargs):
vec2 = np.array([np.cos(np.deg2rad(angle)), np.sin(np.deg2rad(angle))])
xy = np.c_[[length, 0], [0, 0], vec2*length].T + np.array(pos)
ax.plot(*xy.T, color=acol)
return AngleAnnotation(pos, xy[0], xy[2], ax=ax, **kwargs)
# Calculate phasor positions
def calculate_phasors(phi, M):
# Initialize arrays for arrow start and end points
x_start = np.zeros(M)
y_start = np.zeros(M)
x_end = np.zeros(M)
y_end = np.zeros(M)
# Running sum of phasors
x_sum = 0
y_sum = 0
for i in range(M):
# Current phasor
x = np.cos(i * phi)
y = np.sin(i * phi)
# Store start point (end of previous phasor)
x_start[i] = x_sum
y_start[i] = y_sum
# Add current phasor
x_sum += x
y_sum += y
# Store end point
x_end[i] = x_sum
y_end[i] = y_sum
return x_start, y_start, x_end, y_end
x_start, y_start, x_end, y_end = calculate_phasors(phi, M)
plt.figure(figsize=get_size(6, 6))
ax = plt.gca()
for i in range(M):
plt.arrow(x_start[i], y_start[i],
x_end[i]-x_start[i], y_end[i]-y_start[i],
head_width=0.15, head_length=0.2, fc='k', ec='k',
length_includes_head=True,
label=f'E{i+1}' if i == 0 else "")
plt.arrow(0, 0, x_end[-1], y_end[-1],
head_width=0.15, head_length=0.2, fc='r', ec='r',
length_includes_head=True, label='Resultant')
ax.set_aspect('equal')
xx = np.linspace(-1, 3, 100)
ax.plot(xx,(xx-1)*np.tan(phi),'k--',lw=0.5)
ax.plot([1,3],[0,0],'k--',lw=0.5)
kw = dict(size=195, unit="points", text=r"$\Delta \phi$")
plot_angle(ax, (1.0, 0), phi*180/np.pi, textposition="inside", **kw)
plt.axis('off')
max_range = max(abs(x_end[-1]), abs(y_end[-1])) * 1.2
plt.xlim(-0, max_range/1.5)
plt.ylim(-0.1, max_range/1.)
plt.show()
# Parameters
M = 6
phi = np.linspace(-4*np.pi, 4*np.pi, 10000) # increased resolution
I0 = 1
def multiple_beam_pattern(phi, M):
numerator = np.sin(M * phi/2)**2
denominator = np.sin(phi/2)**2
I = np.where(denominator != 0, numerator/denominator, M**2)
return I
I = I0 * multiple_beam_pattern(phi, M)
first_min = 2*np.pi/M # theoretical value
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return array[idx], idx
half_max = M**2/2
phi_positive = phi[phi >= 0] # only positive values
I_positive = I[phi >= 0]
_, idx_half = find_nearest(I_positive, half_max)
half_width = phi_positive[idx_half]
# Create plot
plt.figure(figsize=get_size(10, 6))
plt.plot(phi/np.pi, I, 'b-', label=f'M={M}')
#plt.plot(first_min/np.pi, multiple_beam_pattern(first_min, M), 'ro')
#plt.annotate(f'First minimum\nφ = 2π/M = {first_min/np.pi:.2f}π',
plt.axvline(x=first_min/np.pi, color='r', linestyle='--', label=f'φ = 2π/M = {first_min/np.pi:.2f}π')
#plt.plot(half_width/np.pi, half_max, 'go')
plt.xlabel(r'phase $\Delta \phi/\pi$')
plt.ylabel('intensity I/I₀')
plt.title(f'Multiple Beam Interference Pattern (M={M})')
plt.ylim(0, M**2 + 15)
plt.legend()
plt.show()
The result is therefore an oscillating function. The numerator \(\sin^2(M\Delta \phi/2)\) shows and oscillation frequency, which is by a factor of \(M\) higher than the one in the denominator \(\sin^2 (\Delta \phi/2)\). Therefore the intensity pattern is oscillating rapidly and creating a first minimum at
\[ \Delta \phi=\frac{2\pi}{M} \]
This is an important result, since it shows that the number of sources \(M\) determines the position of the first minimum and the interference peak gets narrower with increasing \(M\). Since the phase difference \(\Delta \phi\) between neighboring sources is the same as for the double slit experiment, i.e. \(\Delta \phi=2\pi d/\lambda \sin(\theta)\), we can also determine the angular position of the first minimum. This is given by
\[ \sin(\theta_\textrm{min})=\frac{1}{M}\frac{\lambda}{d} \]
This again has the common feature that it scales as \(\lambda/d\). A special situation occurs, whenever the numerator and the denominator become zero. This will happen whenever
\[ \Delta \phi=m 2\pi \]
where \(m\) is an integer and denotes the interference order, i.e. the number of wavelength that neighboring partial waves have as path length difference. In this case, the intensity distributiion will give us
\[ I=I_0 \frac{0}{0} \]
and we have to determine the limit with the help of l’Hospitals rule. The outcome of this calculation is, that
\[ I(\Delta \phi=m2\Delta \pi)=M^2 I_0 \]
which can be also realized when using the small angle approximation for the sine functions.
We would like to introduce a different representation of the multiple wave interference of the grating, which is quite insightful. The first order (\(m=1\)) constructive interference condition is given by
\[ \frac{1}{\lambda}\sin{\theta}= \frac{1}{d} \]
which also means that
\[ \frac{2\pi}{\lambda}\sin{\theta}= \frac{2\pi}{d} \]
This can be written as
\[ k \sin{\theta}= K \]
where \(k\) is the magnitude of the wavevector of the light and \(K\) is the wavevector magnitude that corresponds to the grating period \(d\). As the magnitude of the wavevector of the light is conserved, the wavevectors of the incident light and the light traveling along the direction of the first interence peak form the sides of an equilateral triangle. This is shown in the following figure.
import numpy as np
import matplotlib.pyplot as plt
k = 1 # Magnitude of k₁ and k₂
origin = np.array([0, 0])
k1 = np.array([k, 0])
theta_deg = 30 # θ = 30 degrees
theta_rad = np.deg2rad(theta_deg)
k2 = k * np.array([np.cos(theta_rad), np.sin(theta_rad)])
K = k2 - k1
point_O = origin
point_A = point_O + k1
point_B = point_O + k2
plt.figure(figsize=get_size(10, 10))
ax = plt.gca()
# Plot vector k₁
ax.arrow(point_O[0], point_O[1], k1[0], k1[1],
head_width=0.02, head_length=0.03, fc='k', ec='k', length_includes_head=True)
ax.arrow(point_A[0], point_A[1], K[0], K[1],
head_width=0.02, head_length=0.03, fc='b', ec='b', length_includes_head=True)
ax.arrow(point_O[0], point_O[1], k2[0], k2[1],
head_width=0.02, head_length=0.03, fc='k', ec='k', length_includes_head=True)
# Label vectors
ax.text(k1[0]/2 - 0.05, k1[1]/2 - 0.05, r'$\mathbf{k}$', fontsize=14, color='k')
ax.text(point_A[0] + K[0]/2 , point_A[1] + K[1]/2 + 0.05, r'$\mathbf{K}$', fontsize=14, color='b')
ax.text(k2[0]/2 + 0.0, k2[1]/2+0.1, r'$\mathbf{k}$', fontsize=14, color='k')
# Indicate angle θ between k₁ and k₂ at the origin
arc_radius = 0.3 # Radius of the arc representing θ
angle_range = np.linspace(0, theta_rad, 100)
arc_x = arc_radius * np.cos(angle_range)
arc_y = arc_radius * np.sin(angle_range)
ax.plot(arc_x, arc_y, color='k')
ax.text(arc_radius * np.cos(theta_rad / 2) + 0.02,
arc_radius * np.sin(theta_rad / 2) + 0.02,
r'$\theta$', fontsize=14)
# Set equal aspect ratio
ax.set_aspect('equal', adjustable='box')
all_x = [point_O[0], point_A[0], point_B[0]]
all_y = [point_O[1], point_A[1], point_B[1]]
margin = 0.2
ax.set_xlim(min(all_x) - margin, max(all_x) + margin)
ax.set_ylim(min(all_y) - margin, max(all_y) + margin)
plt.axis('off')
# Display the plot
plt.show()
This means that the diffraction grating is providing a wavevector \(K\) to alter the direction of the incident light. This is again a common feature reappearing in many situations as for example in the X-ray diffraction of crystals.
We will turn our attention now to a slight modification of the previous multiwave interference. We will introduce a decreasing amplitude of the individual waves. The first wave shall have an amplitude \(U_1=\sqrt{I_0}\). The next wave, however, will not only be phase shifted but also have a smaller amplitude.
\[ U_2=h U_1 \]
where \(h=re^{i\phi}\) with \(|h|=r<1\). \(r\) can be regarded as a reflection coefficient, which deminishes the amplitude of the incident wave. According to that the intensity is reduced by
\[ I_2=|U_2|^2=|h U_1|^2=r^2 I_1 \]
The intensity of the incident wave is multiplied by a factor \(r^2\), while the amplitude is multiplied by \(r\). Note that the phase factor \(e^{i\Delta\phi}\) is removed when taking the square of this complex number.
The amplitude of the reflected wave is diminished by a factor \(r\le 1\), which is called the reflection coefficient. The intensity is diminished by a factor \(R=|r|^2\le1\), which is the reflectance.
In the absence of absorption, reflectance \(R\) and transmittance \(T\) add to one due to energy conservation.
\[ R+T=1 \]
Consequently, the third wave would be now \(U_3=hU_2=h^2U_1\). The total amplitude is thus
\[ U=U_1+U_2+U_3+\ldots+U_M = \sqrt{I_0}(1+h+h^2+\ldots) \]
M = 18 # number of phasors
phi = np.pi/6 # example phase difference between successive phasors
r = 0.95 # reduction factor for each subsequent phasor
def plot_angle(ax, pos, angle, length=0.95, acol="C0", **kwargs):
vec2 = np.array([np.cos(np.deg2rad(angle)), np.sin(np.deg2rad(angle))])
xy = np.c_[[length, 0], [0, 0], vec2*length].T + np.array(pos)
ax.plot(*xy.T, color=acol)
return AngleAnnotation(pos, xy[0], xy[2], ax=ax, **kwargs)
def calculate_phasors(phi, M, r):
x_start = np.zeros(M)
y_start = np.zeros(M)
x_end = np.zeros(M)
y_end = np.zeros(M)
x_sum = 0
y_sum = 0
for i in range(M):
amplitude = r**i # exponential decrease
x = amplitude * np.cos(i * phi)
y = amplitude * np.sin(i * phi)
x_start[i] = x_sum
y_start[i] = y_sum
x_sum += x
y_sum += y
x_end[i] = x_sum
y_end[i] = y_sum
return x_start, y_start, x_end, y_end
x_start, y_start, x_end, y_end = calculate_phasors(phi, M, r)
plt.figure(figsize=get_size(6, 6),dpi=150)
ax = plt.gca()
for i in range(M):
plt.arrow(x_start[i], y_start[i],
x_end[i]-x_start[i], y_end[i]-y_start[i],
head_width=0.15, head_length=0.2,
fc='k', ec='k',
length_includes_head=True,
label=f'E{i+1}' if i == 0 else "")
plt.arrow(0, 0, x_end[-1], y_end[-1],
head_width=0.15, head_length=0.2, fc='r', ec='r',
length_includes_head=True, label='Resultant')
ax.set_aspect('equal')
xx = np.linspace(-1, 3, 100)
ax.plot(xx,(xx-1)*np.tan(phi),'k--',lw=0.5)
ax.plot([1,3],[0,0],'k--',lw=0.5)
kw = dict(size=195, unit="points", text=r"$\phi$")
plot_angle(ax, (1.0, 0), phi*180/np.pi, textposition="inside", **kw)
plt.axis('off')
max_range = max(abs(x_end[-1]), abs(y_end[-1])) * 1.2
plt.xlim(-max_range/1.8, max_range/0.8)
plt.ylim(-0.1, max_range/0.9)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Create phase array from -2π to 2π
phi = np.linspace(-2*np.pi, 2*np.pi, 1000)
def calculate_intensity(phi, F):
return 1/(1 + 4*(F/np.pi)**2 * np.sin(phi/2)**2)
plt.figure(figsize=get_size(10, 6))
finesse_values = [1, 4, 20]
styles = ['-', '--', ':']
for F, style in zip(finesse_values, styles):
I = calculate_intensity(phi, F)
plt.plot(phi/np.pi, I, style, label=f'$\\mathcal{{F}}={F}$')
plt.xlabel('Phase $\\phi/\\pi$')
plt.ylabel('$I/I_{\\mathrm{max}}$')
plt.grid(True, alpha=0.3)
plt.legend()
plt.ylim(0, 1.1)
plt.show()
This yields again
\[ U=\sqrt{I_0}\frac{(1-h^M)}{1-h}=\frac{\sqrt{I_0}}{1-r e^{i\Delta\phi}} \]
Calculating the intensity of the waves is giving
\[ I=|U|^2=\frac{I_{0}}{|1-re^{i\Delta\phi}|^2}=\frac{I_0}{(1-r)^2+4r\sin^2(\Delta\phi/2)} \]
which is also known as the Airy function. This function can be further simplified by the following abbrevations
\[ I_{\rm max}=\frac{I_0}{(1-r)^2} \]
and
\[ \mathcal{F}=\frac{\pi \sqrt{r}}{1-r} \]
where the latter is called the Finesse. With those abbrevations, we obtain
\[ I=\frac{I_{\rm max}}{1+4\left(\frac{\mathcal{F}}{\pi}\right)^2\sin^{2}(\Delta\phi/2)} \]
for the interference of multiple waves with decreasing amplitude.
This intensity distribution has a different shape than the one we obtained for multiple waves with the same amplitude.
We clearly observe that with increasing Finesse the intensity maxima, which occur at multiples fo \(\pi\) get much narrower. In addition the regions between the maxima show better contrast and fopr higher Finesse we get complete destructive interference.
Let us consider now interference in the time domain. We introduce two monochromatic waves of frequencies \(\nu_1\) and \(\nu_2\). We will denote their amplitudes by \(\sqrt{I_1}\) and \(\sqrt{I_2}\).
The total amplitude is thus
\[ U=U_1+U_2 = \sqrt{I_1} \exp(i2\pi\nu_1 t) + \sqrt{I_2} \exp(i2\pi\nu_2 t) \]
such that we obtain an Intensity
\[ I=|U|^2 = I_1 + I_2 + 2\sqrt{I_1I_2}\cos(2\pi(\nu_1-\nu_2)t) \]
The intensity is thus time dependent and oscillates at a frequency \(\nu_1-\nu_2\), which is the so-called beating frequency. Similar schemes are used in optical heterodyne detection but also in acoustics when tuning your guitar.
Consider now a whole set of \(M=2L+1\) each with an amplitude \(\sqrt{I_0}\). The frequencies of the waves are given by \(\nu_q=\nu_0+q\Delta\nu\) with \(q=-L,\dots,L\) with \(\nu_0\) beeing the center frequency of the spectrum and \(\Delta \nu\) the frequency spacing. We will assume that \(\Delta nu<<\nu_0\) such that the total amplitude of the waves is given by
\[ U=\sum_{q=-L}^L \sqrt{I_0} \exp(i2\pi(\nu_0+q\Delta\nu)t) \]
The total intensity can then be calculated in the same way as for the multiple source in space before. Using \(\phi=2\pi \Delta \nu t\) we obtained \[ I(t)=I_0 \frac{\sin^2(M\pi t/T)}{\sin^2(\pi t/T)} \]
with \(T=1/\Delta\nu\) and a maximum intensity of \(I_{\rm max}=M^2 I_0\).
import numpy as np
import matplotlib.pyplot as plt
# Parameters
Delta_nu = 1e9 # 1 GHz
M = 1000
I0 = 1 # Normalized initial intensity
T = 1/Delta_nu # 1 ns
pulse_width = T/M # 1 ps
# Create figure with two subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=get_size(15, 10),
gridspec_kw={'width_ratios': [2, 1]})
# Time array for main plot showing multiple pulses
t_main = np.linspace(-2*T, 2*T, 20000)
# Time array for inset showing single pulse detail
t_detail = np.linspace(-5*pulse_width, 5*pulse_width, 10000)
# Calculate intensity function (avoiding division by zero)
def intensity(t, M, T, I0):
# Small value to prevent division by zero
eps = 1e-10
# Calculate using the formula
num = np.sin(M * np.pi * t / T)**2
den = np.sin(np.pi * t / T)**2
# Handle points where denominator is near zero
near_zero = np.abs(np.sin(np.pi * t / T)) < eps
result = np.zeros_like(t, dtype=float)
# For normal points, use formula
mask = ~near_zero
result[mask] = I0 * num[mask] / den[mask]
# For points where denominator is near zero, use limit value
result[near_zero] = M**2 * I0
return result
# Calculate intensities
I_main = intensity(t_main, M, T, I0)
I_detail = intensity(t_detail, M, T, I0)
# Maximum intensity
max_intensity = M**2 * I0
# Main plot showing multiple pulses
ax1.plot(t_main*1e9, I_main, 'b-')
ax1.set_xlabel('time [ns]')
ax1.set_ylabel(r' $I/I_{0}$')
ax1.grid(True, alpha=0.3)
# Mark maximum intensity
ax1.axhline(y=max_intensity, color='r', linestyle='--')
# Mark period T
ax1.annotate('', xy=(T*1e9, 0.5e6), xytext=(0, 0.5e6),
arrowprops=dict(arrowstyle='<->', color='g'))
ax1.text(0.5*T*1e9, max_intensity*1.1, f'T = 1/Δν = {T*1e9:.1f} ns',
color='g', ha='center')
ax1.set_ylim(0, max_intensity*1.2)
# Detail plot showing single pulse
ax2.plot(t_detail*1e12, I_detail, 'b-')
ax2.set_xlabel('time [ps]')
ax2.set_ylabel(r'$I/I_{0}$')
# Mark pulse width
ax2.annotate('', xy=(-pulse_width/1*1e12, max_intensity*1.05),
xytext=(pulse_width/1*1e12, max_intensity*1.05),
arrowprops=dict(arrowstyle='<->', color='m'))
ax2.text(0, max_intensity*1.1, f'T/M = {pulse_width*1e12:.1f} ps',
color='m', ha='center')
ax2.set_ylim(0, max_intensity*1.2)
plt.tight_layout()
plt.show()
The pulse generation we just examined leads us to an important concept in modern optics: frequency combs. A frequency comb is a spectrum consisting of a series of discrete, equally spaced frequency lines that results from a train of phase-coherent pulses in the time domain.
Let’s extend our understanding of multiple wave beating to include phase coherence. When we have a set of equidistant frequency components that maintain a fixed phase relationship, the resulting time-domain signal is a periodic train of pulses. Conversely, a periodic train of pulses in the time domain corresponds to a frequency comb in the spectral domain.
The relationship between these domains is described by the Fourier transform. For a pulse train with repetition rate \(f_{rep} = \Delta\nu\) (the spacing between frequency components), the frequency spectrum consists of lines at:
\[f_n = f_0 + n \cdot f_{rep}\]
where \(f_0\) is the carrier-envelope offset frequency and \(n\) is an integer.
Consider a train of pulses described by the electric field:
\[E(t) = \sum_{m=-\infty}^{\infty} A(t - m T_{rep}) e^{i[\omega_c(t - m T_{rep}) + \phi_{CE} \cdot m + \phi_0]}\]
where: - \(A(t)\) is the pulse envelope - \(T_{rep} = 1/f_{rep}\) is the pulse repetition period - \(\omega_c\) is the carrier frequency - \(\phi_{CE}\) is the carrier-envelope phase slip (the phase shift from pulse to pulse) - \(\phi_0\) is a constant phase
The Fourier transform of this pulse train gives us a frequency comb with:
\[f_0 = \frac{\phi_{CE}}{2\pi} \cdot f_{rep}\]
This carrier-envelope offset frequency (\(f_0\)) is crucial for determining the absolute positions of the comb lines.
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq, fftshift
# Parameters
f_rep = 0.1e12 # 100 GHz repetition rate
f_0 = 0.2 * f_rep # Carrier-envelope offset frequency (20% of f_rep)
f_c = 2e12 # 2 THz carrier frequency (simplified for visualization)
pulse_fwhm = 50e-15 # 50 fs pulse width
# Time domain
N = 2**14 # Number of points for good FFT resolution
T_total = 100e-12 # 100 ps total window (multiple pulses)
dt = T_total / N
t = np.linspace(0, T_total, N, endpoint=False)
# Pulse envelope function (Gaussian)
def gaussian_envelope(t, t0, fwhm):
# Convert FWHM to standard deviation
sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))
return np.exp(-((t - t0) ** 2) / (2 * sigma ** 2))
# Generate a train of pulses with carrier-envelope phase evolution
def pulse_train(t, f_rep, f_c, f_0, fwhm):
T_rep = 1/f_rep
omega_c = 2 * np.pi * f_c
phi_CE = 2 * np.pi * f_0 / f_rep # Phase slip per pulse
# Initialize field
E = np.zeros_like(t, dtype=complex)
# Add multiple pulses
for m in range(int(T_total * f_rep) + 5): # +5 to ensure we capture all pulses
t0 = m * T_rep
if t0 > T_total:
break
# Pulse envelope
env = gaussian_envelope(t, t0, fwhm)
# Carrier with phase evolution
carrier = np.exp(1j * (omega_c * (t - t0) + phi_CE * m))
# Add pulse to field
E += env * carrier
return E
# Generate the electric field
E_t = pulse_train(t, f_rep, f_c, f_0, pulse_fwhm)
# Calculate intensity in time domain
I_t = np.abs(E_t)**2
# Calculate the frequency spectrum
E_f = fftshift(fft(E_t))
f = fftshift(fftfreq(N, dt)) # Frequency axis
I_f = np.abs(E_f)**2
# Normalize for better visualization
I_t_norm = I_t / np.max(I_t)
I_f_norm = I_f / np.max(I_f)
# Create figure with two subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=get_size(16, 8))
# Time domain plot
ax1.plot(t*1e12, I_t_norm, 'b-')
ax1.set_xlabel('Time [ps]')
ax1.set_ylabel('Normalized Intensity')
ax1.set_title('Pulse Train (Time Domain)')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 50) # Show a subset of the time window
# Annotate T_rep
t_rep = 1/f_rep
ax1.annotate('', xy=(t_rep*1e12, 0.1), xytext=(0, 0.1),
arrowprops=dict(arrowstyle='<->', color='red'))
ax1.text(t_rep*1e12/2, 0.15, f'$T_{{rep}}$ = {t_rep*1e12:.1f} ps',
ha='center', color='red')
# Frequency domain plot - zoomed around carrier frequency
f_min = f_c - 5*f_rep
f_max = f_c + 5*f_rep
f_mask = (f > f_min) & (f < f_max)
ax2.plot(f[f_mask]*1e-12, I_f_norm[f_mask], 'g-')
ax2.set_xlabel('frequency [THz]')
ax2.set_ylabel('normalized spectral power')
ax2.set_title('Frequency Comb')
ax2.grid(True, alpha=0.3)
# Annotate frequency spacing
max_peaks = np.where(I_f_norm[f_mask] > 0.5)[0]
if len(max_peaks) >= 2:
peak_indices = np.sort(max_peaks)[:2]
peak_freqs = f[f_mask][peak_indices]
ax2.annotate('', xy=(peak_freqs[1]*1e-12, 0.5), xytext=(peak_freqs[0]*1e-12, 0.5),
arrowprops=dict(arrowstyle='<->', color='red'))
ax2.text((peak_freqs[0] + peak_freqs[1])*1e-12/2, 0.6,
f'$f_{{rep}}$ = {f_rep*1e-9:.0f} GHz', ha='center', color='red')
# Mark f_0
ax2.annotate('$f_0$', xy=(f_0*1e-12, 0.3), xytext=(f_0*1e-12, 0.4),
arrowprops=dict(arrowstyle='->', color='purple'))
plt.tight_layout()
plt.show()
Frequency combs have revolutionized precision measurements in physics and enabled numerous applications:
Optical Clocks: Frequency combs provide a “gear mechanism” to count optical frequencies, enabling optical atomic clocks that are orders of magnitude more precise than conventional atomic clocks.
Precision Spectroscopy: The precise and stable frequency references allow for high-resolution molecular spectroscopy.
Calibration of Astronomical Spectrographs: Frequency combs enable the detection of Earth-like exoplanets by providing precise wavelength calibration for astronomical instruments.
Distance Measurements: They enable precise absolute distance measurements used in applications from gravitational wave detectors to satellite formation flying.
Telecommunications: Frequency combs can be used for wavelength-division multiplexing in optical communications.
In practice, frequency combs are often generated using mode-locked lasers. In such lasers, multiple longitudinal modes of the laser cavity oscillate with a fixed phase relationship, resulting in short pulses. The Fourier transform of these regularly spaced pulses is precisely the frequency comb.
The mode-locking can be achieved through various mechanisms:
The 2005 Nobel Prize in Physics was awarded to Theodor W. Hänsch and John L. Hall for their contributions to the development of laser-based precision spectroscopy, including the optical frequency comb technique.
Frequency combs beautifully illustrate the duality between time and frequency domains in physics. A perfectly periodic sequence of events in time (the pulse train) corresponds to discrete, equally spaced frequencies. The more precise and stable the temporal pattern, the more precise and stable the frequency components.
This duality is fundamental to many areas of physics and engineering, from quantum mechanics to signal processing, and frequency combs represent one of its most elegant and useful manifestations in optics.