Particle in a box model of Cyanine Dyes

Particle in a box model of Cyanine Dyes#

We discussed the particle in a box model as being particularly useful for thinking about electronic transitions in cyanine dyes. While this model would also suggest that we can abritrarily red-shift such dyes by increasing the number of conjugated pi units, and that while we do, the cyanine dyes would become better absorbers and emitters.
However, this is not actually true for reasons that are very interesting and somewhat complicated! Dr. Foley has a great podcast discussing these topics in more detail, if you are interested.

That said, the scaling relationships suggested by the particle in a box model are interesting in their own right, and useful in some limited contexts, so we try to illustrate them here!

First, let’s recall the expression for the energy eigenvalues of the particle in a box of length \(L\), \begin{equation} E_n = \frac{\hbar^2 \pi^2 n^2}{2 m L^2} \end{equation}

where following the discussion in this article, we will approximate the length of a cyanine dye as \begin{equation} L = a \left( k + 1 \right) + b \end{equation} where \(a\) = 2.49 Angstroms, \(b\) = 5.69 Angstroms, and \(k\) is the number of -CH=CH- units in the dye.

The relevance of this models comes from the fact that the the absorption behavior is dominated by transitions between the HOMO and LUMO orbitals in these conjugated dyes, and the particle in a box eigenfunctions / eigenvalues provide supprisingly reasonable models for these orbitals. As the number of -CH=CH- units (\(k\)) increases, the quantum number associated with the HOMO increases as follows: \begin{equation} n_{HOMO} = k + 3. \end{equation}

Let’s then explore how the HOMO-LUMO gap (which approximates the absorption energy) scales with the number of these units and the length of the box.

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, IntSlider
from IPython.display import display

# Constants in SI units
a = 2.49e-10  # Average -CH=CH- unit length
b = 5.69e-10  # Terminal moieties length
hbar = 1.054e-34  # Reduced Planck's constant
m = 9.109e-31  # Electron mass
joules_to_ev = 6.242e18  # Conversion factor

def particle_in_box_wavefunction(x, n, L):
    """
    Compute the normalized wavefunction for particle in a box
    
    Parameters:
    x: position array
    n: quantum number
    L: box length
    """
    return np.sqrt(2/L) * np.sin(n * np.pi * x / L)

def update_cyanine_dye(k=3):
    """
    Update and display the cyanine dye properties for given k
    
    Parameters:
    k: number of -CH=CH- units
    """
    # Calculate box length
    L = a * (k + 1) + b
    
    # Quantum numbers
    n_homo = k + 3
    n_lumo = k + 4
    
    # Energy eigenvalues in Joules
    E_homo_J = (hbar**2 * np.pi**2 * n_homo**2) / (2 * m * L**2)
    E_lumo_J = (hbar**2 * np.pi**2 * n_lumo**2) / (2 * m * L**2)
    
    # Convert to eV
    E_homo = E_homo_J * joules_to_ev
    E_lumo = E_lumo_J * joules_to_ev
    
    # HOMO-LUMO gap
    E_gap = E_lumo - E_homo
    
    # Absorption wavelength in nm
    lambda_abs = 1240 / E_gap
    
    # Create position array for wavefunction plotting
    x = np.linspace(0, L, 1000)
    
    # Calculate wavefunctions
    psi_homo = particle_in_box_wavefunction(x, n_homo, L)
    psi_lumo = particle_in_box_wavefunction(x, n_lumo, L)
    
    # Create figure with subplots
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
    
    # Plot 1: Wavefunctions with energy levels
    ax1.axhline(y=E_homo, color='blue', linestyle='--', alpha=0.3, linewidth=1)
    ax1.axhline(y=E_lumo, color='red', linestyle='--', alpha=0.3, linewidth=1)
    
    # Scale wavefunctions for visualization and offset by their energies
    scale_factor = 0.3 * E_gap  # Scale relative to gap
    ax1.plot(x*1e10, E_homo + psi_homo * scale_factor, 'b-', linewidth=2, 
             label=f'HOMO (n={n_homo})')
    ax1.plot(x*1e10, E_lumo + psi_lumo * scale_factor, 'r-', linewidth=2, 
             label=f'LUMO (n={n_lumo})')
    
    # Add arrow showing transition
    mid_x = L * 1e10 / 2
    ax1.annotate('', xy=(mid_x, E_lumo), xytext=(mid_x, E_homo),
                arrowprops=dict(arrowstyle='<->', color='green', lw=2))
    ax1.text(mid_x + L*1e10*0.05, (E_homo + E_lumo)/2, 
             f'ΔE = {E_gap:.3f} eV\nλ = {lambda_abs:.1f} nm',
             fontsize=11, color='green', verticalalignment='center',
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    ax1.set_xlabel('Position (Ã…)', fontsize=12)
    ax1.set_ylabel('Energy (eV)', fontsize=12)
    ax1.set_title(f'HOMO and LUMO Wavefunctions\nk = {k} (-CH=CH- units), L = {L*1e10:.2f} Ã…', 
                  fontsize=13)
    ax1.legend(fontsize=11)
    ax1.grid(True, alpha=0.3)
    ax1.set_xlim(0, L*1e10)
    
    # Plot 2: Scaling relationships
    k_range = np.arange(0, 16)
    L_range = a * (k_range + 1) + b
    n_homo_range = k_range + 3
    n_lumo_range = k_range + 4
    
    E_homo_range = (hbar**2 * np.pi**2 * n_homo_range**2) / (2 * m * L_range**2) * joules_to_ev
    E_lumo_range = (hbar**2 * np.pi**2 * n_lumo_range**2) / (2 * m * L_range**2) * joules_to_ev
    E_gap_range = E_lumo_range - E_homo_range
    lambda_range = 1240 / E_gap_range
    
    ax2.plot(L_range*1e10, E_gap_range, 'go-', linewidth=2, markersize=6, 
             label='HOMO-LUMO Gap')
    ax2.plot(L*1e10, E_gap, 'ro', markersize=12, label=f'Current (k={k})', zorder=5)
    
    ax2.set_xlabel('Box Length (Ã…)', fontsize=12)
    ax2.set_ylabel('Absorption Energy (eV)', fontsize=12)
    ax2.set_title('Scaling of Absorption Energy with Box Length', fontsize=13)
    ax2.legend(fontsize=11)
    ax2.grid(True, alpha=0.3)
    
    # Add secondary y-axis for wavelength
    ax2_twin = ax2.twinx()
    ax2_twin.plot(L_range*1e10, lambda_range, 'mo-', linewidth=2, markersize=6, 
                  alpha=0.6, label='Wavelength')
    ax2_twin.plot(L*1e10, lambda_abs, 'co', markersize=12, alpha=0.8, zorder=5)
    ax2_twin.set_ylabel('Absorption Wavelength (nm)', fontsize=12, color='m')
    ax2_twin.tick_params(axis='y', labelcolor='m')
    
    plt.tight_layout()
    plt.show()
    
    # Print summary
    print(f"{'='*60}")
    print(f"Cyanine Dye Properties (k = {k})")
    print(f"{'='*60}")
    print(f"Box length (L):           {L*1e10:.2f} Ã…")
    print(f"HOMO quantum number:      n = {n_homo}")
    print(f"LUMO quantum number:      n = {n_lumo}")
    print(f"HOMO energy:              {E_homo:.3f} eV")
    print(f"LUMO energy:              {E_lumo:.3f} eV")
    print(f"Transition energy (ΔE):   {E_gap:.3f} eV")
    print(f"Absorption wavelength:    {lambda_abs:.1f} nm")
    print(f"{'='*60}")

# Create interactive slider
interact(update_cyanine_dye, 
         k=IntSlider(min=0, max=15, step=1, value=3, 
                     description='k (units):', 
                     continuous_update=False,
                     style={'description_width': 'initial'}))
<function __main__.update_cyanine_dye(k=3)>