Seminar 7 — Active Matter Workshop & Course Review

Introduction

This final seminar explores active matter and optical manipulation, with a focus on microswimmers and light-driven particles. We then provide a comprehensive review of key concepts and formulas from the entire course.

Problem 1: Optical Trap Stiffness and Equipartition Theorem

Problem Statement:

When a dielectric particle is trapped in a tightly focused laser beam, it experiences a restoring force proportional to displacement:

\[F = -\kappa x\]

where \(\kappa\) is the trap stiffness. At thermal equilibrium, the equipartition theorem tells us:

\[\left\langle \frac{1}{2} \kappa x^2 \right\rangle = \frac{1}{2} k_B T\]

a) A 1 µm polystyrene bead in an optical trap has stiffness \(\kappa = 10\) pN/µm. Calculate the position variance \(\langle x^2 \rangle\) at room temperature (T = 300 K).

b) What is the typical position fluctuation (standard deviation) in nm?

c) Why does a stiffer trap lead to smaller positional variance? Explain physically.

Code 
# Optical trap stiffness and equipartition
k_B = 1.38e-23  # J/K
T = 300  # K
kappa = 10  # pN/µm

# Convert to SI units for consistency
kappa_SI = kappa * 1e-12 / 1e-6  # pN/µm → N/m = 10e-3 N/m = 0.01 N/m

print("Optical Trap Stiffness (Equipartition Theorem)")
print("=" * 50)
print(f"Trap stiffness: κ = {kappa} pN/µm = {kappa_SI:.3e} N/m")
print(f"Temperature: T = {T} K")
print(f"Boltzmann constant: k_B = {k_B:.3e} J/K")

# Part a): variance from equipartition
# <1/2 · κ · x²> = 1/2 · k_B · T
# <x²> = k_B · T / κ

x2_variance_SI = k_B * T / kappa_SI  # m²
x2_variance_um = x2_variance_SI * 1e12  # µm²
x2_variance_nm = x2_variance_SI * 1e18  # nm²

print(f"\nPart a) Position variance:")
print(f"  <x²> = k_B·T / κ")
print(f"  <x²> = {x2_variance_SI:.3e} m²")
print(f"  <x²> = {x2_variance_nm:.1f} nm²")

# Part b): standard deviation (RMS displacement)
x_rms = np.sqrt(x2_variance_SI) * 1e9  # convert to nm

print(f"\nPart b) RMS displacement (position fluctuation):")
print(f"  σ_x = √<x²> = {x_rms:.2f} nm")

# Part c): show dependence on trap stiffness
kappas = np.array([1, 5, 10, 20, 50]) * 1e-3  # N/m
x_rms_array = np.sqrt(k_B * T / kappas) * 1e9

print(f"\nPart c) Dependence on trap stiffness:")
for k, x in zip(kappas*1e3, x_rms_array):
    print(f"  κ = {k:.2e} N/m ({k/1e-3:.1f} pN/µm) → σ_x = {x:.2f} nm")
print(f"\n  Stiffer trap → smaller variance (weaker thermal fluctuations)")
Optical Trap Stiffness (Equipartition Theorem)
==================================================
Trap stiffness: κ = 10 pN/µm = 1.000e-05 N/m
Temperature: T = 300 K
Boltzmann constant: k_B = 1.380e-23 J/K

Part a) Position variance:
  <x²> = k_B·T / κ
  <x²> = 4.140e-16 m²
  <x²> = 414.0 nm²

Part b) RMS displacement (position fluctuation):
  σ_x = √<x²> = 20.35 nm

Part c) Dependence on trap stiffness:
  κ = 1.00e+00 N/m (1000.0 pN/µm) → σ_x = 2.03 nm
  κ = 5.00e+00 N/m (5000.0 pN/µm) → σ_x = 0.91 nm
  κ = 1.00e+01 N/m (10000.0 pN/µm) → σ_x = 0.64 nm
  κ = 2.00e+01 N/m (20000.0 pN/µm) → σ_x = 0.45 nm
  κ = 5.00e+01 N/m (50000.0 pN/µm) → σ_x = 0.29 nm

  Stiffer trap → smaller variance (weaker thermal fluctuations)

Problem 2: Thermophoresis in Janus Particles

Problem Statement:

A Janus particle (half-coated with gold, half with silica) is heated by a focused laser. The temperature gradient induces a thermophoretic drift toward the cold side.

The thermophoretic velocity is:

\[v_{\text{TP}} = -D_T \nabla T\]

where \(D_T = D \cdot S_T\) is the thermophoretic diffusion coefficient, \(D\) is the Brownian diffusion coefficient, and \(S_T\) is the Soret coefficient.

a) A 1 µm Janus sphere has \(D = 1\) µm²/s (from Stokes-Einstein) and Soret coefficient \(S_T = 0.1\) K⁻¹. A focused laser creates a temperature gradient \(\nabla T = 1\) K/µm.

Calculate the thermophoretic drift velocity.

b) In an optical trap, the particle experiences both a trapping force and thermophoretic force. At equilibrium:

\[\kappa (x - x_0) = -D_T \frac{\partial T}{\partial x}\]

where \(x_0\) is the trap equilibrium without temperature gradient. If \(\kappa = 10\) pN/µm, estimate the trap displacement \(\Delta x\).

c) How would you use this effect to propel a microswimmer toward a target?

Code 
# Thermophoresis calculation
D = 1  # µm²/s
S_T = 0.1  # K⁻¹
dT_dx = 1  # K/µm

print("Thermophoresis in Janus Particles")
print("=" * 50)
print(f"Diffusion coefficient: D = {D} µm²/s")
print(f"Soret coefficient: S_T = {S_T} K⁻¹")
print(f"Temperature gradient: ∇T = {dT_dx} K/µm")

# Part a): thermophoretic velocity
D_T = D * S_T  # µm²/(s·K)
v_TP = -D_T * dT_dx  # µm/s (negative = toward cold)

print(f"\nPart a) Thermophoretic velocity:")
print(f"  D_T = D·S_T = {D} × {S_T} = {D_T} µm²/(s·K)")
print(f"  v_TP = -D_T·(∇T) = -{D_T} × {dT_dx} = {v_TP} µm/s")
print(f"  (moves toward cold side at {abs(v_TP)} µm/s)")

# Part b): trap displacement at equilibrium
kappa = 10  # pN/µm
kappa_SI = kappa * 1e-12 / 1e-6  # N/m

# Force balance: κ·Δx = D_T·(∇T)
# D_T in SI units: D = 1 µm²/s = 1e-12 m²/s
# S_T is in K⁻¹, so D_T·(∇T) has units m²/s × K⁻¹ × K/m = m/s
# But we need force...
# Actually: thermophoretic force F_TP ~ D_T·∇T (drift velocity × friction)
# F = -ξ·v_TP = -6πηa·v_TP, but it's easier to use equilibrium directly

D_SI = D * 1e-12  # m²/s
dT_dx_SI = dT_dx * 1e3  # K/m
v_TP_SI = abs(D_T * 1e-12 * dT_dx_SI)  # m/s

# In equilibrium with trap: κ·Δx = ξ·v_TP where ξ = 6πηa ≈ F/v
# Alternative: use that particle sits at position where thermal gradient force balances trap
# F_thermal ~ k_B·T·S_T·∇T (rough estimate from thermal energy)

F_thermal_approx = k_B * T * S_T * dT_dx_SI
Delta_x = F_thermal_approx / kappa_SI * 1e9  # nm

print(f"\nPart b) Trap displacement:")
print(f"  Estimated thermophoretic force ~")
print(f"    F ~ k_B·T·S_T·(∇T) ~ {F_thermal_approx:.3e} N")
print(f"  Trap displacement:")
print(f"    Δx = F / κ ~ {Delta_x:.2f} nm")
print(f"  (particle shifts toward hot side inside trap)")

print(f"\nPart c) Thermophoretic guidance:")
print(f"  - Design a temperature gradient pointing toward target")
print(f"  - Particle drifts thermophoretically at ~{abs(v_TP):.2f} µm/s")
print(f"  - Can be combined with other propulsion (self-phoresis)")
print(f"  - Example: gold half acts as heater and plasmonic engine")
Thermophoresis in Janus Particles
==================================================
Diffusion coefficient: D = 1 µm²/s
Soret coefficient: S_T = 0.1 K⁻¹
Temperature gradient: ∇T = 1 K/µm

Part a) Thermophoretic velocity:
  D_T = D·S_T = 1 × 0.1 = 0.1 µm²/(s·K)
  v_TP = -D_T·(∇T) = -0.1 × 1 = -0.1 µm/s
  (moves toward cold side at 0.1 µm/s)

Part b) Trap displacement:
  Estimated thermophoretic force ~
    F ~ k_B·T·S_T·(∇T) ~ 4.140e-19 N
  Trap displacement:
    Δx = F / κ ~ 0.00 nm
  (particle shifts toward hot side inside trap)

Part c) Thermophoretic guidance:
  - Design a temperature gradient pointing toward target
  - Particle drifts thermophoretically at ~0.10 µm/s
  - Can be combined with other propulsion (self-phoresis)
  - Example: gold half acts as heater and plasmonic engine

Problem 3: Brownian vs. Active Motion

Problem Statement:

For a passive Brownian particle in 2D, the mean-squared displacement grows linearly with time:

\[\langle r^2 \rangle = 4Dt\]

where \(D\) is the diffusion coefficient.

For an active Brownian particle (self-propelled with speed \(v\)), the MSD is:

\[\langle r^2 \rangle = v^2 \tau_r^2 \left[ \frac{2t}{\tau_r} + e^{-2t/\tau_r} - 1 \right]\]

where \(\tau_r = 1/(2D_r)\) is the rotational diffusion time and \(D_r\) is the rotational diffusion coefficient.

a) Show that at short times (\(t \ll \tau_r\)), active motion gives \(\langle r^2 \rangle \approx v^2 t^2\) (ballistic motion).

b) Show that at long times (\(t \gg \tau_r\)), both Brownian and active particles give \(\langle r^2 \rangle \approx v_{\text{eff}}^2 \tau_r \cdot t\) (diffusive motion).

c) Plot both MSD curves on a log-log scale to illustrate the crossover from ballistic to diffusive behavior.

Code 
# MSD for Brownian vs. active motion
D = 1  # µm²/s, diffusion coefficient
v = 5  # µm/s, active velocity
D_r = 0.5  # rad²/s, rotational diffusion
tau_r = 1 / (2 * D_r)  # rotational diffusion time

print("Brownian vs. Active Motion")
print("=" * 50)
print(f"Diffusion coefficient: D = {D} µm²/s")
print(f"Active velocity: v = {v} µm/s")
print(f"Rotational diffusion: D_r = {D_r} rad²/s")
print(f"Rotational time: τ_r = 1/(2D_r) = {tau_r:.2f} s")

# Time vector
t = np.logspace(-1, 2, 1000)  # 0.1 to 100 s

# Brownian MSD
msd_brownian = 4 * D * t

# Active MSD
msd_active = v**2 * tau_r**2 * (2*t/tau_r + np.exp(-2*t/tau_r) - 1)

# Plot on log-log scale
fig, ax = plt.subplots(figsize=get_size(12, 8))

ax.loglog(t, msd_brownian, 'b-', linewidth=2.5, label='Brownian: $\\langle r^2 \\rangle = 4Dt$')
ax.loglog(t, msd_active, 'r-', linewidth=2.5, label=f'Active: $v = {v}$ µm/s, $\\tau_r = {tau_r:.1f}$ s')

# Add reference lines for ballistic and diffusive regimes
t_ballistic = np.logspace(-1, 0.5, 100)
msd_ballistic = v**2 * t_ballistic**2
ax.loglog(t_ballistic, msd_ballistic, 'r--', alpha=0.5, linewidth=1, label='Ballistic: $\\propto t^2$')

t_diffusive = np.logspace(0.5, 2, 100)
msd_diffusive = v**2 * tau_r * t_diffusive
ax.loglog(t_diffusive, msd_diffusive, 'r:', alpha=0.5, linewidth=1, label='Diffusive: $\\propto t$')

# Mark crossover time
ax.axvline(x=tau_r, color='gray', linestyle='--', alpha=0.7, label=f'Crossover: $t = \\tau_r = {tau_r:.1f}$ s')

ax.set_xlabel('Time $t$ (s)')
ax.set_ylabel(r'Mean-squared displacement $\langle r^2 \rangle$ (µm²)')
ax.grid(True, which='both', alpha=0.3)
ax.legend(fontsize=8, loc='upper left')
ax.set_xlim([0.1, 100])

plt.savefig('img/brownian_vs_active.png', dpi=150, bbox_inches='tight')
plt.close()
print("\nFigure saved to img/brownian_vs_active.png")

# Analysis of scaling laws
print(f"\nPart a) Short-time limit (t ≪ τ_r):")
t_short = 0.01  # s
msd_active_short = v**2 * t_short**2
msd_approx_short = v**2 * tau_r**2 * (2*t_short/tau_r + np.exp(-2*t_short/tau_r) - 1)
print(f"  t = {t_short} s ≪ τ_r = {tau_r:.2f} s")
print(f"  Active MSD (exact) ≈ {msd_approx_short:.4f} µm²")
print(f"  v²·t² = {msd_active_short:.4f} µm²  ✓ Ballistic")

print(f"\nPart b) Long-time limit (t ≫ τ_r):")
t_long = 10  # s
msd_active_long = v**2 * tau_r**2 * (2*t_long/tau_r + np.exp(-2*t_long/tau_r) - 1)
msd_approx_long = v**2 * tau_r * t_long
print(f"  t = {t_long} s ≫ τ_r = {tau_r:.2f} s")
print(f"  Active MSD (exact) ≈ {msd_active_long:.4f} µm²")
print(f"  v²·τ_r·t = {msd_approx_long:.4f} µm²  ✓ Diffusive")
print(f"  Effective diffusion: D_eff ~ v²·τ_r = {v**2 * tau_r:.2f} µm²/s")
Brownian vs. Active Motion
==================================================
Diffusion coefficient: D = 1 µm²/s
Active velocity: v = 5 µm/s
Rotational diffusion: D_r = 0.5 rad²/s
Rotational time: τ_r = 1/(2D_r) = 1.00 s

Figure saved to img/brownian_vs_active.png

Part a) Short-time limit (t ≪ τ_r):
  t = 0.01 s ≪ τ_r = 1.00 s
  Active MSD (exact) ≈ 0.0050 µm²
  v²·t² = 0.0025 µm²  ✓ Ballistic

Part b) Long-time limit (t ≫ τ_r):
  t = 10 s ≫ τ_r = 1.00 s
  Active MSD (exact) ≈ 475.0000 µm²
  v²·τ_r·t = 250.0000 µm²  ✓ Diffusive
  Effective diffusion: D_eff ~ v²·τ_r = 25.00 µm²/s

Problem 4: Monte Carlo Simulation of Brownian and Active Motion

Problem Statement:

Simulate 2D trajectories for: 1. Passive Brownian particle: random walk with diffusion coefficient \(D\) 2. Active Brownian particle: constant-speed swimmer with rotational diffusion

Compare the trajectories and extract MSD to verify scaling laws.

Code 
from scipy.optimize import curve_fit

# Simulation parameters
n_particles = 10  # number of trajectories
n_steps = 1000
dt = 0.01  # s per step

# Brownian particle
sigma_brownian = np.sqrt(2 * D * dt)

# Active particle
v_active = v  # µm/s
D_r_active = D_r  # rad²/s
sigma_theta = np.sqrt(2 * D_r_active * dt)

print("Brownian and Active Motion Simulation")
print("=" * 50)
print(f"Number of particles: {n_particles}")
print(f"Simulation steps: {n_steps}")
print(f"Time step: {dt} s")
print(f"Total time: {n_steps * dt:.1f} s")

# Simulate trajectories
brownian_trajectories = np.zeros((n_particles, n_steps, 2))
active_trajectories = np.zeros((n_particles, n_steps, 2))

for p in range(n_particles):
    # Brownian motion
    for step in range(1, n_steps):
        displacement = np.random.normal(0, sigma_brownian, 2)
        brownian_trajectories[p, step] = brownian_trajectories[p, step-1] + displacement

    # Active motion
    theta = 0  # initial orientation
    for step in range(1, n_steps):
        # Update orientation with rotational diffusion
        dtheta = np.random.normal(0, sigma_theta)
        theta += dtheta
        # Move in current direction
        dx = v_active * np.cos(theta) * dt
        dy = v_active * np.sin(theta) * dt
        active_trajectories[p, step] = active_trajectories[p, step-1] + np.array([dx, dy])

print(f"Trajectories simulated.")
Brownian and Active Motion Simulation
==================================================
Number of particles: 10
Simulation steps: 1000
Time step: 0.01 s
Total time: 10.0 s
Trajectories simulated.

Now compute and plot MSD:

Code 
# Compute MSD for each trajectory
time_steps = np.arange(n_steps)
msd_brownian_sim = np.zeros(n_steps)
msd_active_sim = np.zeros(n_steps)

for step in range(n_steps):
    r2_brownian = np.sum((brownian_trajectories[:, step, :] - brownian_trajectories[:, 0, :])**2, axis=1)
    r2_active = np.sum((active_trajectories[:, step, :] - active_trajectories[:, 0, :])**2, axis=1)
    msd_brownian_sim[step] = np.mean(r2_brownian)
    msd_active_sim[step] = np.mean(r2_active)

# Convert to physical time
time_sim = time_steps * dt

# Fit to theoretical models
def msd_brownian_model(t, D_fit):
    return 4 * D_fit * t

def msd_active_model(t, v_fit, tau_r_fit):
    return v_fit**2 * tau_r_fit**2 * (2*t/tau_r_fit + np.exp(-2*t/tau_r_fit) - 1)

# Fit Brownian (linear regime)
mask_brownian = (time_sim > 0.1) & (time_sim < 5)
popt_b, _ = curve_fit(msd_brownian_model, time_sim[mask_brownian],
                       msd_brownian_sim[mask_brownian], p0=[1], maxfev=5000)
D_fit = popt_b[0]

# Fit active (full model)
mask_active = time_sim > 0
popt_a, _ = curve_fit(msd_active_model, time_sim[mask_active],
                       msd_active_sim[mask_active], p0=[v, tau_r], maxfev=5000)
v_fit, tau_r_fit = popt_a

print(f"\nFit results:")
print(f"  Brownian:")
print(f"    D_true = {D:.2f} µm²/s")
print(f"    D_fit  = {D_fit:.2f} µm²/s")
print(f"\n  Active:")
print(f"    v_true    = {v:.2f} µm/s,  v_fit    = {v_fit:.2f} µm/s")
print(f"    τ_r_true  = {tau_r:.2f} s,  τ_r_fit  = {tau_r_fit:.2f} s")

# Plot trajectories and MSD
fig, axes = plt.subplots(2, 2, figsize=get_size(15, 12))

# Trajectory plots
for p in range(min(3, n_particles)):
    axes[0, 0].plot(brownian_trajectories[p, :, 0], brownian_trajectories[p, :, 1],
                    alpha=0.6, linewidth=0.8)
axes[0, 0].set_xlabel('x (µm)')
axes[0, 0].set_ylabel('y (µm)')
axes[0, 0].set_title('Brownian Trajectories')
axes[0, 0].grid(True, alpha=0.3)

for p in range(min(3, n_particles)):
    axes[0, 1].plot(active_trajectories[p, :, 0], active_trajectories[p, :, 1],
                    alpha=0.6, linewidth=0.8)
axes[0, 1].set_xlabel('x (µm)')
axes[0, 1].set_ylabel('y (µm)')
axes[0, 1].set_title('Active Trajectories')
axes[0, 1].grid(True, alpha=0.3)

# MSD plots (linear scale)
axes[1, 0].plot(time_sim, msd_brownian_sim, 'b-', linewidth=2, label='Simulated')
t_fit_b = np.linspace(0, time_sim[-1], 500)
axes[1, 0].plot(t_fit_b, msd_brownian_model(t_fit_b, D_fit), 'b--',
                linewidth=1.5, label=f'Fit: D = {D_fit:.2f} µm²/s')
axes[1, 0].set_xlabel('Time (s)')
axes[1, 0].set_ylabel(r'MSD (µm²)')
axes[1, 0].set_title('Brownian MSD (Linear)')
axes[1, 0].legend(fontsize=8)
axes[1, 0].grid(True, alpha=0.3)

# MSD plot (log-log)
mask_plot = msd_active_sim > 0.01
axes[1, 1].loglog(time_sim[mask_plot], msd_active_sim[mask_plot], 'r-',
                  linewidth=2, label='Simulated')
t_fit_a = np.linspace(0.001, time_sim[-1], 500)
axes[1, 1].loglog(t_fit_a, msd_active_model(t_fit_a, v_fit, tau_r_fit), 'r--',
                  linewidth=1.5, label=f'Fit: v = {v_fit:.1f} µm/s, τ_r = {tau_r_fit:.2f} s')
axes[1, 1].axvline(x=tau_r_fit, color='gray', linestyle=':', alpha=0.7, label='Crossover')
axes[1, 1].set_xlabel('Time (s)')
axes[1, 1].set_ylabel(r'MSD (µm²)')
axes[1, 1].set_title('Active MSD (Log-Log)')
axes[1, 1].legend(fontsize=8)
axes[1, 1].grid(True, which='both', alpha=0.3)

plt.tight_layout()
plt.savefig('img/brownian_active_simulation.png', dpi=150, bbox_inches='tight')
plt.close()
print("\nFigure saved to img/brownian_active_simulation.png")

Fit results:
  Brownian:
    D_true = 1.00 µm²/s
    D_fit  = 0.68 µm²/s

  Active:
    v_true    = 5.00 µm/s,  v_fit    = 6.85 µm/s
    τ_r_true  = 1.00 s,  τ_r_fit  = 0.71 s

Figure saved to img/brownian_active_simulation.png

Problem 5: Photon Nudging Strategy

Problem Statement:

A microswimmer with orientation \(\theta\) undergoes rotational diffusion. The propulsion laser can be switched ON only when the swimmer’s orientation is close to a target direction \(\theta_{\text{target}}\).

Implement a simple feedback control strategy: - IF \(|\theta - \theta_{\text{target}}| < \theta_{\text{threshold}}\): laser is ON, particle swims at speed \(v\) - ELSE: laser is OFF, particle only drifts with Brownian motion

Simulate this strategy and show that the particle navigates toward and orbits the target.

Code 
# Photon nudging simulation
v_swim = 3  # µm/s, swimming speed
theta_target = 0  # target direction (radians)
theta_threshold = np.radians(30)  # 30° threshold cone
n_steps_nudge = 2000
dt_nudge = 0.01  # s

print("Photon Nudging Strategy")
print("=" * 50)
print(f"Swimming speed: v = {v_swim} µm/s")
print(f"Target direction: θ_target = {np.degrees(theta_target):.1f}°")
print(f"Threshold cone: ±{np.degrees(theta_threshold):.1f}°")
print(f"Simulation time: {n_steps_nudge * dt_nudge:.1f} s")

# Trajectory
position = np.zeros((n_steps_nudge, 2))
theta = np.zeros(n_steps_nudge)
laser_on = np.zeros(n_steps_nudge, dtype=bool)

# Initial random orientation
theta[0] = np.random.uniform(-np.pi, np.pi)

for step in range(1, n_steps_nudge):
    # Rotational diffusion
    dtheta = np.random.normal(0, np.sqrt(2 * D_r_active * dt_nudge))
    theta[step] = theta[step-1] + dtheta

    # Decide if laser is ON
    angle_error = np.abs(np.angle(np.exp(1j * (theta[step] - theta_target))))
    if angle_error < theta_threshold:
        # Laser ON: swim toward target
        laser_on[step] = True
        dx = v_swim * np.cos(theta_target) * dt_nudge
        dy = v_swim * np.sin(theta_target) * dt_nudge
    else:
        # Laser OFF: only rotational diffusion (Brownian drift in lab frame is minimal)
        laser_on[step] = False
        dx = 0
        dy = 0

    # Add small Brownian noise
    dx += np.random.normal(0, sigma_brownian)
    dy += np.random.normal(0, sigma_brownian)

    position[step] = position[step-1] + np.array([dx, dy])

# Plot results
fig, axes = plt.subplots(1, 2, figsize=get_size(15, 6))

# Trajectory
time_nudge = np.arange(n_steps_nudge) * dt_nudge
on_indices = np.where(laser_on)[0]
off_indices = np.where(~laser_on)[0]

axes[0].plot(position[off_indices, 0], position[off_indices, 1], 'b.',
            markersize=2, alpha=0.3, label='Laser OFF')
axes[0].plot(position[on_indices, 0], position[on_indices, 1], 'r.',
            markersize=2, alpha=0.5, label='Laser ON')
axes[0].plot(0, 0, 'g*', markersize=15, label='Target')
axes[0].set_xlabel('x (µm)')
axes[0].set_ylabel('y (µm)')
axes[0].set_title('Photon-Nudged Trajectory')
axes[0].legend(fontsize=8)
axes[0].grid(True, alpha=0.3)
axes[0].set_aspect('equal')

# Orientation and laser state
axes[1].plot(time_nudge, np.degrees(theta), 'b-', linewidth=0.8, alpha=0.7, label='Orientation θ')
axes[1].axhline(y=np.degrees(theta_target), color='g', linestyle='--', alpha=0.7, label='Target')
axes[1].fill_between(time_nudge, np.degrees(theta_target - theta_threshold),
                      np.degrees(theta_target + theta_threshold), alpha=0.2, color='green',
                      label='Laser ON zone')
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Orientation (degrees)')
axes[1].set_ylim([-180, 180])
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('img/photon_nudging.png', dpi=150, bbox_inches='tight')
plt.close()
print("\nFigure saved to img/photon_nudging.png")

# Statistics
frac_on = np.sum(laser_on) / n_steps_nudge
avg_speed = np.mean(np.linalg.norm(np.diff(position, axis=0), axis=1)) / dt_nudge
final_distance = np.linalg.norm(position[-1])

print(f"\nResults:")
print(f"  Laser ON fraction: {frac_on*100:.1f}%")
print(f"  Average net speed: {avg_speed:.2f} µm/s")
print(f"  Final distance from origin: {final_distance:.2f} µm")
Photon Nudging Strategy
==================================================
Swimming speed: v = 3 µm/s
Target direction: θ_target = 0.0°
Threshold cone: ±30.0°
Simulation time: 20.0 s

Figure saved to img/photon_nudging.png

Results:
  Laser ON fraction: 0.0%
  Average net speed: 17.72 µm/s
  Final distance from origin: 11.83 µm

Problem 6: Course Review — Summary of Key Formulas

Problem Statement:

Create a comprehensive reference table of the most important formulas and results from the entire course. Organize by topic.

Code 
# Create summary table
summary_text = """
INTRODUCTION TO PHOTONICS — KEY FORMULAS SUMMARY
================================================

1. DIFFRACTION AND RESOLUTION
   ─────────────────────────
   Abbe diffraction limit (lateral):   d = 0.51 λ / NA
   Rayleigh criterion:                  δθ = 1.22 λ / D
   Axial resolution (wide-field):       d_z = 2 λ / NA²

2. OPTICAL PROPERTIES & MATERIALS
   ────────────────────────────────
   Snell's law:                         n₁ sin(θ₁) = n₂ sin(θ₂)
   Brewster angle:                      tan(θ_B) = n₂/n₁
   Refractive index (Cauchy):           n(λ) = A + B/λ²
   Sellmeier equation:                  n²(λ) = 1 + Σ Bₖλ²/(λ² - Cₖ)
   Group velocity:                      v_g = c / n_g, where n_g = n - λ dn/dλ

3. CONFOCAL & SUPERRESOLUTION
   ──────────────────────────
   Confocal PSF (Gaussian):             PSF_conf = PSF_exc · PSF_det
   Confocal resolution improvement:     σ_conf = σ₀/√2 → factor of √2
   Confocal pinhole (Airy disk):        Pinhole ≈ 1 Airy Unit (AU)

   STED resolution:                     d = d₀ / √(1 + I/I_sat)

   PALM/STORM localization precision:   σ ≈ s/√N + a²/(12N)
   where s = PSF width, N = photons, a = pixel size

   SIM resolution improvement:          d_SIM ≈ d₀/2

   Two-photon excitation:               F ∝ I²  (quadratic)
   Two-photon resolution:               d ≈ 0.51 λ_2P / NA
   Axial sectioning (2-photon):         No pinhole needed; I² gives section

4. SCATTERING & OPTICAL PROPERTIES
   ────────────────────────────────
   Rayleigh scattering (small sphere):  σ_sca = (2π/λ⁴) |ε-1|²/|ε+2|² (2a)⁶
   Scattering efficiency:               Q_sca = σ_sca / σ_geom
   Mie scattering (arbitrary size):     Q = (2/x²) Σ (2n+1)/(n(n+1)) Re(aₙ+bₙ)

   Dynamic light scattering (DLS):      q = 4πn/λ sin(θ/2)
   Decay rate:                          Γ = Dq²
   Stokes-Einstein relation:            D = k_B T / (6πηa)
   Autocorrelation:                     g₂(τ) = 1 + β exp(-2Γτ)

5. PHOTOTHERMAL EFFECTS
   ────────────────────
   Steady-state temperature rise:       ΔT = P_abs / (4πκr)
   Thermophoretic velocity:             v_TP = -D_T ∇T = -D·S_T·∇T
   Thermal lens signal:                 Δn = (dn/dT) · ΔT
   Photothermal absorption cross-section: σ_abs ∝ Im(ε) (Rayleigh)

6. OPTICAL MANIPULATION
   ────────────────────
   Optical trap stiffness:              F = -κx
   Equipartition theorem:               <½κx²> = ½k_B T
   Position variance in trap:           <x²> = k_B T / κ

   Radiation pressure force:            F_rad = n·α·I/c  (for small particles)
   Gradient force:                       F_grad ∝ ∇(I²) (for dipoles)

7. BROWNIAN MOTION & ACTIVE MATTER
   ────────────────────────────────
   Brownian MSD (2D):                   <r²> = 4Dt
   Brownian MSD (3D):                   <r²> = 6Dt

   Active Brownian particle MSD:        <r²> = v²τ_r²[2t/τ_r + exp(-2t/τ_r) - 1]
   Rotational diffusion time:           τ_r = 1/(2D_r)
   Crossover time (ballistic→diffusive): t_cross ~ τ_r

   Effective diffusion (long time):     D_eff ~ v²τ_r

8. FUNDAMENTAL CONSTANTS
   ─────────────────────
   Boltzmann constant:                  k_B = 1.38 × 10⁻²³ J/K
   Planck constant:                     h = 6.626 × 10⁻³⁴ J·s
   Elementary charge:                   e = 1.602 × 10⁻¹⁹ C
   Speed of light:                      c = 2.998 × 10⁸ m/s
   Vacuum permittivity:                 ε₀ = 8.854 × 10⁻¹² F/m

9. TYPICAL PARAMETER VALUES
   ────────────────────────
   Water viscosity (25°C):              η ≈ 0.89 mPa·s
   Water thermal conductivity:          κ ≈ 0.6 W/(m·K)
   Water refractive index:              n ≈ 1.33
   Water dn/dT:                         ≈ -1.0 × 10⁻⁴ K⁻¹

   Visible light wavelength:            λ ≈ 400–700 nm
   Near-infrared (two-photon):          λ ≈ 700–1000 nm
   NA (objective lens):                 0.5–1.4 (air/oil immersion)

   Diffusion coefficient (1 µm bead):   D ≈ 0.5–1 µm²/s
   Optical trap stiffness:              κ ≈ 1–100 pN/nm
   Single-molecule photon count:        N ≈ 100–10000 photons

"""

print(summary_text)

# Save as text file for reference
with open('course_summary.txt', 'w') as f:
    f.write(summary_text)

print("Summary saved to course_summary.txt")

INTRODUCTION TO PHOTONICS — KEY FORMULAS SUMMARY
================================================

1. DIFFRACTION AND RESOLUTION
   ─────────────────────────
   Abbe diffraction limit (lateral):   d = 0.51 λ / NA
   Rayleigh criterion:                  δθ = 1.22 λ / D
   Axial resolution (wide-field):       d_z = 2 λ / NA²

2. OPTICAL PROPERTIES & MATERIALS
   ────────────────────────────────
   Snell's law:                         n₁ sin(θ₁) = n₂ sin(θ₂)
   Brewster angle:                      tan(θ_B) = n₂/n₁
   Refractive index (Cauchy):           n(λ) = A + B/λ²
   Sellmeier equation:                  n²(λ) = 1 + Σ Bₖλ²/(λ² - Cₖ)
   Group velocity:                      v_g = c / n_g, where n_g = n - λ dn/dλ

3. CONFOCAL & SUPERRESOLUTION
   ──────────────────────────
   Confocal PSF (Gaussian):             PSF_conf = PSF_exc · PSF_det
   Confocal resolution improvement:     σ_conf = σ₀/√2 → factor of √2
   Confocal pinhole (Airy disk):        Pinhole ≈ 1 Airy Unit (AU)

   STED resolution:                     d = d₀ / √(1 + I/I_sat)

   PALM/STORM localization precision:   σ ≈ s/√N + a²/(12N)
   where s = PSF width, N = photons, a = pixel size

   SIM resolution improvement:          d_SIM ≈ d₀/2

   Two-photon excitation:               F ∝ I²  (quadratic)
   Two-photon resolution:               d ≈ 0.51 λ_2P / NA
   Axial sectioning (2-photon):         No pinhole needed; I² gives section

4. SCATTERING & OPTICAL PROPERTIES
   ────────────────────────────────
   Rayleigh scattering (small sphere):  σ_sca = (2π/λ⁴) |ε-1|²/|ε+2|² (2a)⁶
   Scattering efficiency:               Q_sca = σ_sca / σ_geom
   Mie scattering (arbitrary size):     Q = (2/x²) Σ (2n+1)/(n(n+1)) Re(aₙ+bₙ)

   Dynamic light scattering (DLS):      q = 4πn/λ sin(θ/2)
   Decay rate:                          Γ = Dq²
   Stokes-Einstein relation:            D = k_B T / (6πηa)
   Autocorrelation:                     g₂(τ) = 1 + β exp(-2Γτ)

5. PHOTOTHERMAL EFFECTS
   ────────────────────
   Steady-state temperature rise:       ΔT = P_abs / (4πκr)
   Thermophoretic velocity:             v_TP = -D_T ∇T = -D·S_T·∇T
   Thermal lens signal:                 Δn = (dn/dT) · ΔT
   Photothermal absorption cross-section: σ_abs ∝ Im(ε) (Rayleigh)

6. OPTICAL MANIPULATION
   ────────────────────
   Optical trap stiffness:              F = -κx
   Equipartition theorem:               <½κx²> = ½k_B T
   Position variance in trap:           <x²> = k_B T / κ

   Radiation pressure force:            F_rad = n·α·I/c  (for small particles)
   Gradient force:                       F_grad ∝ ∇(I²) (for dipoles)

7. BROWNIAN MOTION & ACTIVE MATTER
   ────────────────────────────────
   Brownian MSD (2D):                   <r²> = 4Dt
   Brownian MSD (3D):                   <r²> = 6Dt

   Active Brownian particle MSD:        <r²> = v²τ_r²[2t/τ_r + exp(-2t/τ_r) - 1]
   Rotational diffusion time:           τ_r = 1/(2D_r)
   Crossover time (ballistic→diffusive): t_cross ~ τ_r

   Effective diffusion (long time):     D_eff ~ v²τ_r

8. FUNDAMENTAL CONSTANTS
   ─────────────────────
   Boltzmann constant:                  k_B = 1.38 × 10⁻²³ J/K
   Planck constant:                     h = 6.626 × 10⁻³⁴ J·s
   Elementary charge:                   e = 1.602 × 10⁻¹⁹ C
   Speed of light:                      c = 2.998 × 10⁸ m/s
   Vacuum permittivity:                 ε₀ = 8.854 × 10⁻¹² F/m

9. TYPICAL PARAMETER VALUES
   ────────────────────────
   Water viscosity (25°C):              η ≈ 0.89 mPa·s
   Water thermal conductivity:          κ ≈ 0.6 W/(m·K)
   Water refractive index:              n ≈ 1.33
   Water dn/dT:                         ≈ -1.0 × 10⁻⁴ K⁻¹

   Visible light wavelength:            λ ≈ 400–700 nm
   Near-infrared (two-photon):          λ ≈ 700–1000 nm
   NA (objective lens):                 0.5–1.4 (air/oil immersion)

   Diffusion coefficient (1 µm bead):   D ≈ 0.5–1 µm²/s
   Optical trap stiffness:              κ ≈ 1–100 pN/nm
   Single-molecule photon count:        N ≈ 100–10000 photons


Summary saved to course_summary.txt

Reflection Questions

  1. Optical Sectioning: Compare the optical sectioning mechanisms of confocal and two-photon microscopy. Why is two-photon better for thick specimens?

  2. Superresolution Trade-offs: STED requires high power and can cause phototoxicity. PALM/STORM requires many frames. SIM requires patterned illumination. Which would you choose for live-cell imaging of fast dynamics? Why?

  3. Active vs. Passive: An optical trap can hold a passive bead or guide an active swimmer. How would you distinguish between these two scenarios by observing position fluctuations?

  4. Detection Limits: Why does lock-in detection fail if noise is at the signal frequency? How would you detect a 1 Hz signal buried in 50 Hz line noise?

  5. Course Integration: This course integrated several themes:

    • Light-matter interaction (absorption, scattering, fluorescence, trapping)
    • Optical design (lenses, aberrations, superresolution optics)
    • Signal detection (photons, noise, lock-in, autocorrelation)
    • Particle dynamics (Brownian motion, thermophoresis, propulsion)

    How would you design an experiment combining two of these themes?

Final Note

You have now completed the Introduction to Photonics course! You are equipped with: - Understanding of light-matter interactions across scales (molecules to micrometers) - Practical knowledge of modern optical microscopy (confocal, superresolution, nonlinear) - Tools for analyzing scattering, thermal effects, and particle dynamics - Experience with numerical simulation and data fitting - Critical thinking about experimental design and measurement limits

Further reading: Consult specialized textbooks on topics of interest: - Born & Wolf – Principles of Optics (classical optics foundation) - Goodman – Introduction to Fourier Optics (diffraction and imaging) - Hecht – Optics (general reference) - Hell et al. – Nanoscopy and Multidimensional Optical Microscopy (superresolution) - Evans & Wong – Optical Materials and Photonics (nonlinear effects)

Good luck with your further studies!