Interference in Space and Time

Two-Wave Interference

Interference arises whenever two or more wave solutions to the same linear wave equation are superimposed. The key results are familiar from experimental physics; we recall them here with an eye toward what is needed for multi-wave interference and Fourier optics.

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}\]

Writing \(U_j=\sqrt{I_j}\,e^{i\phi_j}\) and collecting terms, the total intensity is:

\[ I=I_1+I_2+2\sqrt{I_1 I_2}\cos(\Delta \phi) \tag{1}\]

where \(\Delta \phi=\phi_2-\phi_1\) is the phase difference. For equal intensities \(I_1=I_2=I_0\):

\[ I=4I_0\cos^2\!\left(\frac{\Delta \phi}{2}\right) \]

giving maximum intensity \(4I_0\) at \(\Delta\phi = 2\pi m\) and zero at \(\Delta\phi = (2m+1)\pi\).

Code 
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,axes=plt.subplots(2,3,figsize=get_size(14,6))

# Constructive
axes[0,0].plot(x/wavelength,y1); axes[0,0].set_ylabel(r"$U_1$")
axes[0,1].plot(x/wavelength,y1); axes[0,1].set_ylabel(r"$U_2$")
axes[0,2].plot(x/wavelength,2*y1); axes[0,2].set_ylabel(r"$U_1+U_2$")
axes[0,0].set_title(r"Constructive ($\Delta\phi=0$)")

# Destructive
axes[1,0].plot(x/wavelength,y1); axes[1,0].set_ylabel(r"$U_1$"); axes[1,0].set_xlabel(r"distance [$\lambda$]")
axes[1,1].plot(x/wavelength,y2); axes[1,1].set_ylabel(r"$U_2$"); axes[1,1].set_xlabel(r"distance [$\lambda$]")
axes[1,2].plot(x/wavelength,y1+y2); axes[1,2].set_ylabel(r"$U_1+U_2$"); axes[1,2].set_ylim(-1,1); axes[1,2].set_xlabel(r"distance [$\lambda$]")
axes[1,0].set_title(r"Destructive ($\Delta\phi=\pi$)")

plt.tight_layout()
plt.show()

Constructive interference (top) and destructive interference (bottom).

Phase and Path Difference

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\):

\[\phi_j(\mathbf{r},t) = \mathbf{k}_j\cdot\mathbf{r} - \omega t + \phi_{0j}\]

The time-independent phase difference is then:

\[ \Delta\phi(\mathbf{r}) = (\mathbf{k}_2-\mathbf{k}_1)\cdot\mathbf{r} + (\phi_{02}-\phi_{01}) \]

For waves traveling along the same direction, this reduces to \(\Delta\phi = k\Delta s + \Delta\phi_0\).

Phase Difference and Path Difference

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 constructive interference.

Interference of Waves in Space

Code 
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()

Interference of two plane waves propagating under an angle of 45°. The two left graphs show the original waves. The two right show the total amplitude and the intensity pattern.

A particularly important case is the interference between a known reference wave \(U_R\) (e.g. a plane wave) and an unknown object wave \(U_O\) with complicated wavefronts. Expanding the interference formula:

\[ I = |U_R + U_O|^2 = \underbrace{|U_R|^2 + |U_O|^2}_{\text{backgrounds}} + \underbrace{U_R^*\, U_O + U_R\, U_O^*}_{\text{cross terms}} \]

For a plane-wave reference with uniform amplitude \(|U_R| = 1\), the cross term \(U_R^* U_O\) is simply \(U_O\) multiplied by a known phase \(e^{-i\mathbf{k}_R\cdot\mathbf{r}}\). The intensity pattern \(I(\mathbf{r})\) therefore encodes the full complex amplitude of \(U_O\) – both its magnitude and its phase – in the spacing and shape of the fringes.

Code 
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()

Interference of a spherical wave and a plane wave. The top graphs show the original waves. The two bottom show the total amplitude and the intensity pattern.

The fringe pattern in the intensity image traces the wavefronts of the spherical wave: where the curved wavefronts of \(U_O\) coincide with the flat wavefronts of \(U_R\), the path difference is zero and fringes are bright; where they differ by \(\lambda/2\), fringes are dark. The pattern is therefore a geometric record of \(U_O\)’s phase structure.

The Holography Principle

Recording. Illuminate the object with coherent light; let the scattered object wave \(U_O\) interfere with a reference wave \(U_R\) on a recording medium. The recorded intensity is

\[I = |U_R|^2 + |U_O|^2 + U_R^*\,U_O + U_R\,U_O^*\]

Reconstruction. Illuminate the developed recording with \(U_R\) again. The transmitted field contains the term

\[U_R \cdot U_R^*\,U_O = |U_R|^2\, U_O \propto U_O\]

which is the original object wave, propagating in exactly its original direction. An observer looking through the hologram sees the object as if it were physically present, including parallax. A fourth term \(U_R^2 U_O^*\) produces a conjugate (pseudoscopic) image traveling in the mirror direction.

The key requirement is coherence: the path difference between reference and object beams must be smaller than the coherence length \(L_c\), so that the fringes have good contrast and faithfully encode the phase.

A super nice website to try out interference interactively is here.

Multi-Wave Interference

So far we looked at the interference of two waves. In practice, many partial waves contribute simultaneously – this is the generic situation for gratings, etalons, and the plane-wave decomposition in Fourier optics. We treat two cases:

  • multi-wave interference with constant amplitude (gratings, slits)
  • multi-wave interference with decreasing amplitude (Fabry-Perot etalon)

Multiple Wave Interference with Constant Amplitude

In the case of constant amplitude (realized by a grating, for example), the total wave amplitude is:

\[ U=U_1+U_2+U_3+\ldots+U_M \]

where we sum \(M\) partial waves, each separated from its neighbor by a phase difference \(\Delta \phi\). The \(p\)-th wave is:

\[ U_p=\sqrt{I_0}e^{i(p-1)\Delta \phi} \]

The total amplitude is a geometric sum:

\[ U=\sqrt{I_0}\left (1+h+h^2+\ldots +h^{M-1}\right), \qquad h=e^{i\Delta \phi} \]

Applying the geometric-series formula and computing \(|U|^2\):

\[ I=I_{0}\frac{\sin^2(M\Delta \phi/2)}{\sin^2(\Delta \phi/2)} \tag{2}\]

Code 
# 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.axvline(x=first_min/np.pi, color='r', linestyle='--', label=f'φ = 2π/M = {first_min/np.pi:.2f}π')

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()

Multiple wave interference of \(M=6\) waves with a phase difference of \(\phi=\pi/8\). The black arrows represent the individual waves, the red arrow the sum of all waves.
Figure 1— Multiple beam interference pattern for M=6 beams. The intensity distribution is shown as a function of the phase shift \(\phi\). The first minimum is at \(\phi=2\pi/M\). The intensity distribution is symmetric around \(\phi=0\).

The numerator oscillates \(M\) times faster than the denominator, creating sharp principal maxima. The first minimum after each maximum is at:

\[ \Delta \phi=\frac{2\pi}{M} \]

so the peaks narrow as \(M\) increases. At the principal maxima (\(\Delta\phi = m\,2\pi\), \(m\) integer), l’Hopital’s rule gives \(I = M^2 I_0\). For a grating with slit spacing \(d\), the phase difference is \(\Delta\phi = 2\pi d\sin\theta/\lambda\), so the angular width of the \(m\)-th order scales as \(\lambda/(Md)\).

Wavevector Representation

The constructive interference condition for the first order (\(m=1\)) reads:

\[ \frac{2\pi}{\lambda}\sin{\theta}= \frac{2\pi}{d} \quad \Longleftrightarrow \quad k \sin{\theta}= K \]

where \(K = 2\pi/d\) is the grating wavevector. Because \(|\mathbf{k}|\) is conserved (elastic scattering), the incident and diffracted wavevectors differ by exactly \(\mathbf{K}\):

Code 
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()

Wavevector summation for the diffraction grating. The wavevector of the incident light \(k\) and the wavevector of the light traveling along the direction of the first interference peak \(K\) form an equilateral triangle.

The grating supplies a wavevector \(\mathbf{K}\) that shifts the incident wavevector to a new direction. This picture – a periodic structure in real space providing discrete momentum kicks in \(\mathbf{k}\)-space – recurs throughout Fourier optics: a spatial frequency component at \(K = 2\pi/d\) in the aperture function routes light into the direction \(\sin\theta = K/k\).

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.

Intensity at Boundaries

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) \]

Code 
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()

Phase construction of a multiwave intereference with M waves with decreasing amplitude due to a reflection coefficient \(r=0.95\).
Code 
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()

Multiple wave interference with decreasing amplitude. The graph shows the intensity distribution over the phase angle \(\phi\) for different values of the Finesse \(\mathcal{F}\).

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_{\mathrm{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_{\mathrm{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 of \(\pi\) get much narrower. In addition the regions between the maxima show better contrast and for higher Finesse we get complete destructive interference.

Coherence

In the earlier consideration we obtained a general description for the phase difference between two waves. It 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.

Two waves of different frequency over time.

The above image shows the timetrace of the amplitude of two waves with slightly different frequency. Due to the frequency difference, the waves run out of phase and have acquired a phase difference of \(\pi\) after \(40\) fs.

The temporal coherence of two waves is defined by the time it takes for them to accumulate a phase difference of \(2\pi\). For two waves at frequencies \(\nu_1 = \nu_0 - \Delta\nu/2\) and \(\nu_2 = \nu_0 + \Delta\nu/2\):

\[ \Delta \phi = 2\pi \Delta \nu \,\Delta t \]

Setting \(\Delta\phi = 2\pi\) defines the coherence time:

\[ \tau_{c}=\frac{1}{\Delta \nu} \]

and the coherence length \(L_c = c\tau_c\). Monochromatic light (\(\Delta\nu = 0\)) is perfectly coherent; broadband (white) light has an extremely short coherence length.

Coherence

Two waves are called coherent if they exhibit a fixed phase relation in space or time. It measures their ability to interfere. The main types of coherence are:

Temporal Coherence

  • Measures phase correlation of a wave with itself at different times
  • Characterized by coherence time \(\tau_c\) and coherence length \(L_c = c\tau_c\)
  • Related to spectral width: \(\tau_c = 1/\Delta\nu\)
  • Perfect for monochromatic waves (single frequency)
  • Limited for broad spectrum sources (like thermal light)

Spatial Coherence

  • Measures phase correlation between different points in space
  • Important for interference from extended sources
  • Determines ability to form interference patterns
  • Related to source size and geometry

Coherence is a property of the light source and is connected to the frequency distribution of the light. Sources can be:

  • Fully coherent: ideal laser
  • Partially coherent: real laser
  • Incoherent: thermal light

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:

  1. Temporal Coherence The complex degree of temporal coherence is given by:

\[g^{(1)}(\tau) = \frac{\langle U(t)U^*(t+\tau)\rangle}{\sqrt{\langle|U(t)|^2\rangle\langle|U(t+\tau)|^2\rangle}}\]

where:

  • \(\tau\) is the time delay
  • \(U(t)\) is the electric field
  • \(\langle...\rangle\) denotes time averaging
  1. Spatial Coherence Similarly, spatial coherence between two points is characterized by:

\[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:

  • \(|g^{(1)}| = 1\) indicates perfect coherence
  • \(|g^{(1)}| = 0\) indicates complete incoherence
  • \(0 < |g^{(1)}| < 1\) indicates partial coherence

A finite coherence time and length leads to partial coherence, affecting interference visibility through:

  • Reduced contrast in interference patterns
  • Limited coherence length/area
  • Spectral broadening
Code 
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', fontsize=7)


# 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', fontsize=7)
plt.tight_layout()
plt.subplots_adjust(right=0.85)

plt.show()

Temporal correlation for two waves with slightly different frequencies. The vertical line indicates the coherence time τc = π/Δω.

Besides different frequencies, coherence can also be destroyed by phase jumps. The following example shows two waves with the same frequency but multiple random phase jumps.

Code 
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', fontsize=7)

# 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()

Temporal correlation for two waves of same frequency showing decoherence due to multiple phase jumps. Vertical lines indicate positions of phase jumps.

Temporal Interference (supplementary)

Beating of two waves

Consider two monochromatic waves of frequencies \(\nu_1\) and \(\nu_2\). The total amplitude is:

\[ U=\sqrt{I_1} \exp(i2\pi\nu_1 t) + \sqrt{I_2} \exp(i2\pi\nu_2 t) \]

giving a time-dependent intensity:

\[ I=I_1 + I_2 + 2\sqrt{I_1I_2}\cos(2\pi(\nu_1-\nu_2)t) \]

The intensity oscillates at the beat frequency \(\nu_1-\nu_2\). This scheme underlies optical heterodyne detection and acoustic tuning.

Multiple wave beating and pulse generation

Consider \(M=2L+1\) waves with frequencies \(\nu_q=\nu_0+q\Delta\nu\), \(q=-L,\dots,L\). Using the same geometric-sum result as Eq. 2 with \(\phi=2\pi \Delta \nu t\):

\[ I(t)=I_0 \frac{\sin^2(M\pi t/T)}{\sin^2(\pi t/T)}, \qquad T=\frac{1}{\Delta\nu} \]

The intensity forms a periodic pulse train with period \(T\) and peak intensity \(M^2 I_0\). Each pulse has width \(\approx T/M\).

Code 
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()
Figure 2— Multiple wave beating with M=1000 monochromatic waves separated by Δν=1 GHz. The intensity oscillates with period T=1/Δν=1 ns. Each pulse has a width of approximately T/M=1 ps with maximum intensity I_max=M²I₀.
Frequency Combs: Phase-Coherent Temporal Interference

The pulse generation we just examined leads to an important concept in modern optics: frequency combs. A frequency comb is a spectrum consisting of discrete, equally spaced frequency lines arising from a phase-coherent pulse train in the time domain.

When the frequency components maintain a fixed phase relationship, the Fourier transform of the periodic pulse train gives lines at:

\[f_n = f_0 + n \cdot f_{rep}\]

where \(f_0\) is the carrier-envelope offset frequency and \(f_{rep} = \Delta\nu\) is the line spacing. The offset frequency is:

\[f_0 = \frac{\phi_{CE}}{2\pi} \cdot f_{rep}\]

where \(\phi_{CE}\) is the pulse-to-pulse carrier-envelope phase slip.

Code 
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()
Figure 3— Demonstration of a frequency comb. (Left) Time domain representation showing a train of phase-coherent pulses. (Right) Frequency domain representation showing equally spaced frequency lines forming a comb structure.

Frequency combs are the time-frequency Fourier transform pair made physical: a perfectly periodic pulse train in time produces discrete, equally spaced lines in frequency. Applications include optical atomic clocks, precision molecular spectroscopy, exoplanet Doppler calibration, and absolute distance measurement. The 2005 Nobel Prize in Physics was awarded to Theodor W. Hansch and John L. Hall for their role in developing the optical frequency comb technique.