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.

The two-wave interference formula Eq. 1 is the working principle of the Mach–Zehnder interferometer (MZI), one of the most important building blocks of integrated photonics. An input wave is split equally into two waveguide arms by a 3 dB coupler (or Y-junction), accumulates a phase \(\phi_1\) and \(\phi_2\) along each arm, and is then recombined in a second coupler. The intensity at one of the two output ports follows directly from Eq. 1,

\[ I_\text{out}=\frac{I_\text{in}}{2}\left[1+\cos(\Delta\phi)\right],\qquad \Delta\phi=\phi_2-\phi_1, \]

so that a small change in the optical path length of one arm switches the device between bright and dark output. Because \(\Delta\phi\) can be tuned electrically — through the thermo-optic effect (a microheater changing \(n\) via temperature) or the electro-optic effect (carrier injection or the Pockels effect in LiNbO₃ and thin-film lithium niobate) — the MZI acts as a compact, high-speed intensity modulator.

On silicon photonic chips, MZIs are routinely used as:

  • electro-optic modulators for optical communication (data rates of 50–100 Gbit/s per channel),
  • switches and routers in reconfigurable photonic circuits,
  • biosensors, where a change of refractive index on the surface of one arm shifts the output fringes,
  • programmable mesh elements for photonic neural networks and quantum information processing, where a 2D lattice of MZIs implements an arbitrary unitary transformation on many optical modes.

The MZI is therefore the integrated-photonics counterpart of the bulk two-beam interferometer: the same \(\cos^2(\Delta\phi/2)\) transfer function, realized in a waveguide footprint of a few hundred micrometers.

Figure 1— Monolithically integrated Mach–Zehnder interferometer biochip. (a) Cross-section of the silicon chip with an avalanche-LED light source coupled into a Si₃N₄ rib waveguide; the inset shows the canonical MZI layout (Y-splitter → two arms of length \(L\) → Y-combiner) that converts a phase difference into an oscillating output intensity. (b) Layout of an array of 10 MZIs sharing a common output toward a spectrometer. (c) Photograph of the fabricated chip with the microfluidic cover for biosensing. Image from Psarouli et al., Sci. Rep. 5, 17600 (2015), doi:10.1038/srep17600, reproduced under CC BY 4.0.

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 2— 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\).

Phase-Modulated Grating and Wavefront Sensing (QLSI)

A standard grating treats all slits equally. A richer situation arises when neighboring slits carry an alternating phase shift \(\phi_0\): even-indexed slits transmit with phase 0, odd-indexed slits with phase \(\phi_0\). This is the operating principle of the quadriwave lateral shearing interferometer (QLSI), where \(\phi_0 = \pi\) and the grating is extended to 2D.

For \(M = 2N\) slits with slit spacing \(d\), splitting the sum over even and odd slits gives:

\[ U = \sqrt{I_0}\!\left[\sum_{q=0}^{N-1} e^{i2q\Delta\phi} + e^{i\phi_0}\sum_{q=0}^{N-1} e^{i(2q+1)\Delta\phi}\right] = \sqrt{I_0}\bigl(1 + e^{i(\Delta\phi+\phi_0)}\bigr)\sum_{q=0}^{N-1}e^{i2q\Delta\phi} \]

Taking \(|U|^2\) and using \(|1+e^{i\alpha}|^2 = 4\cos^2(\alpha/2)\):

\[ \boxed{I(\Delta\phi) = 4I_0\cos^2\!\!\left(\frac{\Delta\phi+\phi_0}{2}\right)\cdot\frac{\sin^2(N\Delta\phi)}{\sin^2(\Delta\phi)}} \tag{3}\]

The intensity is a product of two factors:

  • \(\sin^2(N\Delta\phi)/\sin^2(\Delta\phi)\): the usual multi-beam pattern with principal maxima at \(\Delta\phi = m\cdot 2\pi\) (i.e. \(d\sin\theta = m\lambda\))
  • \(4\cos^2((\Delta\phi+\phi_0)/2)\): a slowly varying envelope controlled entirely by \(\phi_0\)

For \(\phi_0=0\) the envelope equals 1 at all principal maxima and the standard result is recovered. For \(\phi_0=\pi\) (the QLSI case):

\[ \cos^2\!\!\left(\frac{\Delta\phi+\pi}{2}\right) = \sin^2\!\!\left(\frac{\Delta\phi}{2}\right) \]

which is zero at every even order (\(\Delta\phi = 2\pi m\)) — including the 0th order — and one at every odd order (\(\Delta\phi = (2m+1)\pi\), i.e. \(d\sin\theta = (m+\tfrac{1}{2})\lambda\)). All energy is redirected from the even orders into the odd orders.

Code 
import numpy as np
import matplotlib.pyplot as plt

M = 20          # total slits
N = M // 2      # pairs
dphi = np.linspace(-3*np.pi, 3*np.pi, 8000)
phi0_values = [0, np.pi/2, np.pi]
labels = [r'$\phi_0=0$ (standard)', r'$\phi_0=\pi/2$', r'$\phi_0=\pi$ (QLSI)']
colors = ['steelblue', 'darkorange', 'crimson']

def I_phase_grating(dphi, N, phi0):
    envelope = 4 * np.cos((dphi + phi0) / 2)**2
    num = np.sin(N * dphi)**2
    den = np.sin(dphi)**2
    fast = np.where(np.abs(den) > 1e-12, num / den, N**2)
    return envelope * fast

fig, axes = plt.subplots(len(phi0_values), 1, figsize=get_size(14, 12), sharex=True)

for ax, phi0, label, color in zip(axes, phi0_values, labels, colors):
    I = I_phase_grating(dphi, N, phi0)
    ax.plot(dphi / np.pi, I / M**2, color=color, lw=0.9)
    ax.set_ylabel(r'$I / I_{\max}$')
    ax.set_title(label, fontsize=9)
    ax.set_ylim(-0.05, 1.15)
    # Mark integer orders
    for m in range(-2, 3):
        ax.axvline(m * 2, color='gray', lw=0.5, ls='--', alpha=0.5)
    # Mark half-integer orders
    for m in [-3, -1, 1, 3]:
        ax.axvline(m, color='gray', lw=0.5, ls=':', alpha=0.4)

axes[-1].set_xlabel(r'phase $\Delta\phi\,/\,\pi$')
axes[-1].set_xticks(np.arange(-6, 7))
axes[-1].set_xticklabels([f'{x}' for x in np.arange(-6, 7)])

plt.tight_layout()
plt.show()
Figure 3— Intensity patterns for a 2N=20-slit grating with alternating phase shift φ₀. For φ₀=0 (standard grating) principal maxima fall at integer orders. For φ₀=π (QLSI) the 0th and all even orders are completely suppressed and energy concentrates in the ±1st orders.

Wavevector picture. The alternating phase makes every two slits a unit cell of period \(2d\). The effective grating wavevector is \(K_\text{eff} = 2\pi/(2d) = K/2\), so the first-order condition becomes \(k\sin\theta = K/2\), placing the primary beams at \(\sin\theta = \lambda/(2d)\) — half the angle of the standard first order. The physical slit spacing \(d\) provides the modulation envelope that kills the even orders.

Connection to QLSI Wavefront Sensing

In Quadriwave Lateral Shearing Interferometry, a 2D crossed grating with a checkerboard \(\pi\) phase pattern (period \(d\) in \(x\) and \(y\), sign alternating as \((-1)^{m+n}\) for slit \((m,n)\)) is placed in the beam. The \(\phi_0=\pi\) condition suppresses the \(0^\text{th}\) order and concentrates diffracted light into four spots at \((\pm 1, 0)\) and \((0, \pm 1)\) first orders. Each spot is a laterally sheared copy of the incident wavefront, shifted by \(d/2\) in the corresponding direction. Interfering these four copies produces a fringe pattern whose local phase encodes the wavefront gradient \(\nabla W(x,y)\). Integration of the gradient then recovers the wavefront \(W(x,y)\) – the same Fourier-domain operation we will study in the Fourier optics lectures.

The QLSI result is a special case of a much more general idea: an extra phase profile \(\phi(x)\) imposed within each period of a periodic structure controls how the energy is distributed among the diffraction orders, while the period \(d\) alone fixes their angular positions through \(d\sin\theta_m = m\lambda\).

In Eq. 3 the multi-slit factor \(\sin^2(N\Delta\phi)/\sin^2(\Delta\phi)\) sets the positions of the principal maxima, and the prefactor \(4\cos^2((\Delta\phi+\phi_0)/2)\) — the Fourier transform of the two-element phase pattern \(\{1, e^{i\phi_0}\}\) — sets their relative weights. Choosing \(\phi_0=\pi\) zeroes the even orders (QLSI). Choosing a different \(\phi(x)\) within a period redistributes the weights differently.

Blazed gratings. A blazed grating takes this idea to its extreme: instead of a binary \(0/\pi\) step, it imposes a linear phase ramp

\[ \phi(x) = \frac{2\pi}{d}\,x \sin\theta_B,\qquad 0\le x<d, \]

within each period, repeating sawtooth-like from period to period. This is exactly the phase profile of a plane wave tilted by the blaze angle \(\theta_B\). Its per-period Fourier transform is concentrated on a single diffraction order — the one whose angle matches \(\theta_B\) — so essentially all the diffracted energy is steered into that order, while all others (including the 0th) are suppressed. Physically the linear phase ramp is realised by tilting the reflecting (or transmitting) facets of each groove by \(\theta_B\), so that specular reflection from the facet and diffraction from the grating point in the same direction.

Acousto-optic modulators (AOMs). The phase profile does not have to be etched into the substrate — it can also be written dynamically by sound. In an AOM, an RF transducer launches an ultrasonic wave at frequency \(\Omega\) into a transparent crystal (TeO₂, fused silica, PbMoO₄…). Through the photoelastic effect, the strain modulates the refractive index as

\[ \Delta n(x,t) = \Delta n_0 \cos(Kx-\Omega t),\qquad K = \frac{2\pi}{\Lambda} = \frac{\Omega}{v_s}, \]

where \(\Lambda\) is the acoustic wavelength and \(v_s\) the speed of sound. A light beam traversing an interaction length \(L\) picks up a sinusoidal phase grating

\[ \phi(x,t) = k_0\,\Delta n_0\, L\,\cos(Kx-\Omega t), \]

which is exactly the same kind of per-period phase modulation as in QLSI or a blazed grating — only now the profile is sinusoidal, dynamic, and moving. The Fourier coefficients of \(e^{i\phi_0\cos\theta}\) are Bessel functions \(J_m(\phi_0)\), so the diffracted orders carry intensities \(\propto J_m^2(\phi_0)\) and their weights are tunable in real time via the RF drive amplitude.

Two features are unique to the acousto-optic case:

  • Tunable deflection angle. Because \(\Lambda = v_s/\Omega\), changing the RF frequency \(\Omega\) rescales the grating period and hence the diffraction angle \(\sin\theta_m = m\lambda/\Lambda\). AOMs are therefore widely used as fast, electronically steered beam deflectors, optical switches, and pulse pickers.
  • Doppler shift. Because the grating moves with the sound wave, photon–phonon energy conservation shifts the \(m\)-th order in frequency by \(m\Omega\): the diffracted light is frequency-shifted by the RF drive, which is the basis of laser frequency stabilisation, heterodyne detection, and laser cooling of atoms.

In the thick-grating (Bragg) regime only a single order survives — just as for a blazed grating — typically with \(>80\%\) diffraction efficiency.

The full hierarchy now reads:

Per-period phase \(\phi(x)\) Realisation Diffraction effect
Constant (\(\phi_0=0\)) Plain amplitude grating 0th order strongest, symmetric pattern
Two-step \(0,\pi\) Etched \(\pi\)-checkerboard Even orders suppressed (QLSI)
Linear ramp Tilted groove facets All energy into one order (blazed grating)
Sinusoidal, travelling Acoustic wave in a crystal Tunable, frequency-shifting AOM

In each case, an engineered sub-period phase profile decides where — and at what frequency — the diffracted light goes.

Multiple Wave Interference with Decreasing Amplitude

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.

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)

Light beating and frequency combs

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 4— 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 5— 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.