Chapter 1: Fundamentals of Spectroscopy

Introduction

Spectroscopy is the study of the interaction between electromagnetic radiation and matter. By analyzing how materials absorb, emit, or scatter light across different wavelengths, we can determine their chemical composition, electronic structure, molecular bonding, and physical properties. This chapter establishes the fundamental principles that underpin all spectroscopic techniques.

1. Light-Matter Interactions

1.1 The Nature of Electromagnetic Radiation

Light exhibits wave-particle duality, behaving both as electromagnetic waves and as discrete packets of energy called photons. The key relationships connecting wavelength, frequency, and energy are:

$$c = \lambda \nu$$

$$E = h\nu = \frac{hc}{\lambda}$$

where:

• $c = 2.998 \times 10^{8}$ m/s is the speed of light
• $\lambda$ is the wavelength
• $\nu$ is the frequency
• $h = 6.626 \times 10^{-34}$ J s is Planck's constant
• $E$ is the photon energy

In spectroscopy, we also frequently use the wavenumber $\tilde{\nu}$, defined as the number of wavelengths per unit length:

$$\tilde{\nu} = \frac{1}{\lambda} = \frac{\nu}{c}$$

Wavenumber is commonly expressed in cm$^{-1}$ and is directly proportional to energy, making it convenient for infrared spectroscopy.

"""
Spectroscopy Unit Converter
Convert between wavelength, frequency, wavenumber, and energy units.
"""
import numpy as np

# Physical constants
h = 6.62607015e-34   # Planck's constant (J s)
c = 2.99792458e8     # Speed of light (m/s)
eV = 1.602176634e-19 # Electron volt (J)

def wavelength_to_frequency(wavelength_nm):
    """Convert wavelength (nm) to frequency (Hz)."""
    wavelength_m = wavelength_nm * 1e-9
    return c / wavelength_m

def wavelength_to_wavenumber(wavelength_nm):
    """Convert wavelength (nm) to wavenumber (cm^-1)."""
    wavelength_cm = wavelength_nm * 1e-7
    return 1.0 / wavelength_cm

def wavelength_to_energy_eV(wavelength_nm):
    """Convert wavelength (nm) to energy (eV)."""
    wavelength_m = wavelength_nm * 1e-9
    energy_J = h * c / wavelength_m
    return energy_J / eV

def wavenumber_to_energy_eV(wavenumber_cm):
    """Convert wavenumber (cm^-1) to energy (eV)."""
    # E = hc * wavenumber (converting cm^-1 to m^-1)
    energy_J = h * c * wavenumber_cm * 100
    return energy_J / eV

# Example conversions
wavelengths = [250, 400, 550, 700, 1000, 10000]  # nm
print("Wavelength Conversion Table")
print("=" * 70)
print(f"{'Wavelength (nm)':<18} {'Frequency (THz)':<18} {'Wavenumber (cm-1)':<18} {'Energy (eV)':<12}")
print("-" * 70)

for wl in wavelengths:
    freq = wavelength_to_frequency(wl) / 1e12  # THz
    wn = wavelength_to_wavenumber(wl)
    energy = wavelength_to_energy_eV(wl)
    print(f"{wl:<18} {freq:<18.2f} {wn:<18.0f} {energy:<12.4f}")

# Output:
# Wavelength (nm)    Frequency (THz)    Wavenumber (cm-1)  Energy (eV)
# ----------------------------------------------------------------------
# 250                1199.17            40000              4.9594
# 400                749.48             25000              3.0996
# 550                545.08             18182              2.2543
# 700                428.27             14286              1.7712
# 1000               299.79             10000              1.2398
# 10000              29.98              1000               0.1240

1.2 Absorption, Emission, and Scattering

When light interacts with matter, three fundamental processes can occur:

Absorption occurs when the photon energy matches the energy difference between two quantum states of the material. The photon is absorbed, and the system transitions to an excited state:

$$\Delta E = E_{\text{final}} - E_{\text{initial}} = h\nu$$

Emission is the reverse process: an excited system releases energy as a photon when transitioning to a lower energy state. Types include:

• Fluorescence: Rapid emission (nanoseconds) from singlet excited states
• Phosphorescence: Slower emission (microseconds to seconds) from triplet states

Scattering involves the redirection of light by the material:

• Rayleigh scattering: Elastic scattering with no energy change (scattered photon has same energy as incident)
• Raman scattering: Inelastic scattering where energy is exchanged with molecular vibrations

"""
Visualize the electromagnetic spectrum and spectroscopic techniques.
"""
import numpy as np
import matplotlib.pyplot as plt

# Define spectral regions and their characteristics
regions = {
    'X-ray': {'range': (0.01, 10), 'unit': 'nm', 'technique': 'XPS, XRF',
              'transition': 'Core electrons'},
    'UV': {'range': (10, 400), 'unit': 'nm', 'technique': 'UV-Vis',
           'transition': 'Valence electrons'},
    'Visible': {'range': (400, 700), 'unit': 'nm', 'technique': 'UV-Vis',
                'transition': 'Valence electrons'},
    'Near-IR': {'range': (700, 2500), 'unit': 'nm', 'technique': 'NIR',
                'transition': 'Overtones/combinations'},
    'Mid-IR': {'range': (2.5, 25), 'unit': 'um', 'technique': 'FTIR',
               'transition': 'Molecular vibrations'},
    'Far-IR': {'range': (25, 1000), 'unit': 'um', 'technique': 'THz',
               'transition': 'Lattice vibrations'},
    'Microwave': {'range': (1, 100), 'unit': 'mm', 'technique': 'ESR',
                  'transition': 'Rotations, spin'}
}

# Create visualization of energy ranges
fig, ax = plt.subplots(figsize=(14, 6))

# Energy scale (eV) - logarithmic
energies = np.logspace(-4, 4, 1000)

# Color mapping for visible spectrum
def wavelength_to_rgb(wavelength):
    """Convert wavelength (nm) to approximate RGB color."""
    if 380 <= wavelength < 440:
        r, g, b = (440 - wavelength) / 60, 0, 1
    elif 440 <= wavelength < 490:
        r, g, b = 0, (wavelength - 440) / 50, 1
    elif 490 <= wavelength < 510:
        r, g, b = 0, 1, (510 - wavelength) / 20
    elif 510 <= wavelength < 580:
        r, g, b = (wavelength - 510) / 70, 1, 0
    elif 580 <= wavelength < 645:
        r, g, b = 1, (645 - wavelength) / 65, 0
    elif 645 <= wavelength <= 780:
        r, g, b = 1, 0, 0
    else:
        r, g, b = 0.5, 0.5, 0.5
    return (r, g, b)

# Plot spectral regions as colored bars
y_positions = [1, 1, 1, 1, 1, 1, 1]
colors = ['#9b59b6', '#3498db', 'rainbow', '#e74c3c', '#f39c12', '#2ecc71', '#1abc9c']
labels = ['X-ray', 'UV', 'Visible', 'Near-IR', 'Mid-IR', 'Far-IR', 'Microwave']

# Energy ranges in eV
energy_ranges = [
    (100, 10000),   # X-ray
    (3.1, 100),     # UV
    (1.77, 3.1),    # Visible
    (0.5, 1.77),    # Near-IR
    (0.05, 0.5),    # Mid-IR
    (0.001, 0.05),  # Far-IR
    (0.00001, 0.001) # Microwave
]

for i, (label, (e_min, e_max)) in enumerate(zip(labels, energy_ranges)):
    if label == 'Visible':
        # Create rainbow gradient for visible region
        n_colors = 100
        visible_energies = np.linspace(1.77, 3.1, n_colors)
        for j in range(n_colors - 1):
            wavelength = 1239.8 / visible_energies[j]  # E(eV) to wavelength(nm)
            color = wavelength_to_rgb(wavelength)
            ax.axvspan(visible_energies[j], visible_energies[j+1],
                      ymin=0.3, ymax=0.7, color=color, alpha=0.8)
    else:
        color = colors[i]
        ax.axvspan(e_min, e_max, ymin=0.3, ymax=0.7, color=color, alpha=0.6, label=label)

ax.set_xscale('log')
ax.set_xlim(1e-5, 1e4)
ax.set_ylim(0, 2)
ax.set_xlabel('Energy (eV)', fontsize=12)
ax.set_title('Electromagnetic Spectrum and Spectroscopic Techniques', fontsize=14, fontweight='bold')

# Add technique labels
techniques = [
    (1000, 1.4, 'XPS\nXRF'),
    (10, 1.4, 'UV-Vis'),
    (2.4, 1.4, 'Visible'),
    (1, 1.4, 'NIR'),
    (0.1, 1.4, 'FTIR'),
    (0.01, 1.4, 'THz'),
    (0.0001, 1.4, 'ESR/NMR')
]

for x, y, text in techniques:
    ax.annotate(text, xy=(x, y), ha='center', fontsize=9, fontweight='bold')

# Add transition type labels
transitions = [
    (1000, 0.1, 'Core electrons'),
    (10, 0.1, 'Valence electrons'),
    (0.1, 0.1, 'Molecular vibrations'),
    (0.0001, 0.1, 'Rotations/Spin')
]

for x, y, text in transitions:
    ax.annotate(text, xy=(x, y), ha='center', fontsize=8, style='italic')

ax.set_yticks([])
plt.tight_layout()
plt.savefig('electromagnetic_spectrum.png', dpi=150, bbox_inches='tight')
plt.show()

print("Spectral regions and their characteristics saved.")

2. The Electromagnetic Spectrum in Spectroscopy

2.1 Spectral Regions and Molecular Transitions

Different regions of the electromagnetic spectrum probe different types of molecular and electronic transitions:

2.2 Energy Level Separation

The Born-Oppenheimer approximation allows us to separate molecular energy into electronic, vibrational, and rotational components:

$$E_{\text{total}} = E_{\text{electronic}} + E_{\text{vibrational}} + E_{\text{rotational}}$$

The typical energy scales are:

• $E_{\text{electronic}} \sim 1 - 10$ eV (UV-Vis region)
• $E_{\text{vibrational}} \sim 0.05 - 0.5$ eV (IR region)
• $E_{\text{rotational}} \sim 10^{-4} - 10^{-3}$ eV (microwave region)

"""
Visualize molecular energy levels and spectroscopic transitions.
"""
import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(12, 8))

# Ground state electronic level
ground_e = 0
# Excited state electronic level
excited_e = 4.0  # eV

# Vibrational levels (harmonic oscillator approximation)
# E_v = hv(v + 1/2), where v = 0, 1, 2, ...
vib_spacing = 0.15  # eV (typical for molecular vibrations)
n_vib_levels = 5

# Draw electronic states
ax.axhline(y=ground_e, xmin=0.1, xmax=0.4, color='blue', linewidth=3, label='Ground electronic state (S0)')
ax.axhline(y=excited_e, xmin=0.6, xmax=0.9, color='red', linewidth=3, label='Excited electronic state (S1)')

# Draw vibrational levels on ground state
for v in range(n_vib_levels):
    y = ground_e + vib_spacing * (v + 0.5)
    ax.axhline(y=y, xmin=0.12, xmax=0.38, color='blue', linewidth=1, alpha=0.7)
    ax.text(0.41, y, f'v={v}', fontsize=9, va='center')

# Draw vibrational levels on excited state
for v in range(n_vib_levels):
    y = excited_e + vib_spacing * (v + 0.5)
    ax.axhline(y=y, xmin=0.62, xmax=0.88, color='red', linewidth=1, alpha=0.7)
    ax.text(0.91, y, f"v'={v}", fontsize=9, va='center')

# Draw transitions
# Electronic transition (UV-Vis)
ax.annotate('', xy=(0.75, excited_e + 0.5*vib_spacing),
            xytext=(0.25, ground_e + 0.5*vib_spacing),
            arrowprops=dict(arrowstyle='->', color='purple', lw=2))
ax.text(0.5, 2.2, 'Electronic\ntransition\n(UV-Vis)', ha='center', fontsize=10, color='purple')

# Vibrational transition (IR)
ax.annotate('', xy=(0.25, ground_e + 1.5*vib_spacing),
            xytext=(0.25, ground_e + 0.5*vib_spacing),
            arrowprops=dict(arrowstyle='->', color='green', lw=2))
ax.text(0.08, ground_e + vib_spacing, 'Vibrational\n(IR)', ha='center', fontsize=9, color='green')

# Fluorescence (emission)
ax.annotate('', xy=(0.25, ground_e + 0.5*vib_spacing),
            xytext=(0.75, excited_e + 0.5*vib_spacing),
            arrowprops=dict(arrowstyle='->', color='orange', lw=2, linestyle='--'))
ax.text(0.5, 3.5, 'Fluorescence\n(emission)', ha='center', fontsize=10, color='orange')

# Add energy scale
ax.set_ylabel('Energy (eV)', fontsize=12)
ax.set_xlim(0, 1)
ax.set_ylim(-0.5, 5.5)
ax.set_xticks([])

# Labels
ax.text(0.25, -0.3, 'Ground State (S0)', ha='center', fontsize=11, fontweight='bold')
ax.text(0.75, -0.3, 'Excited State (S1)', ha='center', fontsize=11, fontweight='bold')

ax.set_title('Molecular Energy Levels and Spectroscopic Transitions', fontsize=14, fontweight='bold')
ax.legend(loc='upper left')

plt.tight_layout()
plt.savefig('energy_levels.png', dpi=150, bbox_inches='tight')
plt.show()

# Print energy scale comparison
print("\nEnergy Scale Comparison:")
print("=" * 50)
print(f"Electronic transition: ~{excited_e:.1f} eV ({1239.8/excited_e:.0f} nm)")
print(f"Vibrational transition: ~{vib_spacing:.2f} eV ({1239.8/vib_spacing:.0f} nm = {vib_spacing*8065:.0f} cm-1)")
print(f"Rotational transition: ~0.001 eV ({1239.8/0.001:.0f} um)")

3. Energy Quantization and Transitions

3.1 Quantum States and Energy Levels

In quantum mechanics, molecules can only exist in discrete energy states. The energy of a photon must exactly match the energy difference between two states for absorption or emission to occur:

$$\Delta E = E_f - E_i = h\nu$$

This quantization explains why spectra show discrete lines or bands rather than continuous absorption.

3.2 The Harmonic Oscillator Model

For molecular vibrations, the harmonic oscillator approximation gives quantized energy levels:

$$E_v = h\nu_0\left(v + \frac{1}{2}\right)$$

where $v = 0, 1, 2, ...$ is the vibrational quantum number and $\nu_0$ is the fundamental vibrational frequency. The vibrational frequency depends on the force constant $k$ and reduced mass $\mu$:

$$\nu_0 = \frac{1}{2\pi}\sqrt{\frac{k}{\mu}}$$

"""
Calculate vibrational frequencies for diatomic molecules using the harmonic oscillator model.
"""
import numpy as np

# Physical constants
h = 6.62607015e-34   # J s
c = 2.99792458e8     # m/s
amu = 1.66054e-27    # kg (atomic mass unit)

def calculate_vibrational_frequency(k, m1, m2):
    """
    Calculate vibrational frequency for a diatomic molecule.

    Parameters:
    -----------
    k : float
        Force constant (N/m)
    m1, m2 : float
        Atomic masses (amu)

    Returns:
    --------
    freq_hz : float
        Frequency in Hz
    wavenumber : float
        Wavenumber in cm^-1
    """
    # Convert masses to kg
    m1_kg = m1 * amu
    m2_kg = m2 * amu

    # Calculate reduced mass
    mu = (m1_kg * m2_kg) / (m1_kg + m2_kg)

    # Calculate frequency
    freq_hz = (1 / (2 * np.pi)) * np.sqrt(k / mu)
    wavenumber = freq_hz / (c * 100)  # Convert to cm^-1

    return freq_hz, wavenumber, mu

# Example molecules with typical force constants
molecules = [
    ('H-H', 575, 1.008, 1.008),      # H2
    ('H-Cl', 516, 1.008, 35.45),     # HCl
    ('C-O', 1902, 12.01, 16.00),     # CO
    ('N-N', 2297, 14.01, 14.01),     # N2
    ('O-O', 1177, 16.00, 16.00),     # O2
    ('C-H', 500, 12.01, 1.008),      # C-H stretch (typical)
    ('C=O', 1200, 12.01, 16.00),     # C=O stretch
]

print("Vibrational Frequencies of Diatomic Molecules")
print("=" * 70)
print(f"{'Molecule':<10} {'Force const (N/m)':<18} {'Reduced mass (amu)':<20} {'Wavenumber (cm-1)':<15}")
print("-" * 70)

for name, k, m1, m2 in molecules:
    freq, wn, mu = calculate_vibrational_frequency(k, m1, m2)
    mu_amu = mu / amu
    print(f"{name:<10} {k:<18} {mu_amu:<20.3f} {wn:<15.0f}")

# Demonstrate effect of isotope substitution
print("\n\nIsotope Effect: H-Cl vs D-Cl")
print("-" * 50)
_, wn_hcl, _ = calculate_vibrational_frequency(516, 1.008, 35.45)
_, wn_dcl, _ = calculate_vibrational_frequency(516, 2.014, 35.45)
print(f"H-Cl: {wn_hcl:.0f} cm-1")
print(f"D-Cl: {wn_dcl:.0f} cm-1")
print(f"Ratio: {wn_hcl/wn_dcl:.3f} (theoretical sqrt(2) = {np.sqrt(2):.3f})")

# Output shows isotope effect - heavier isotopes have lower vibrational frequencies

4. Selection Rules and Transition Probabilities

4.1 Transition Dipole Moment

Not all transitions between energy levels are allowed. The probability of a transition depends on the transition dipole moment:

$$\boldsymbol{\mu}_{fi} = \langle \psi_f | \hat{\boldsymbol{\mu}} | \psi_i \rangle = \int \psi_f^* \hat{\boldsymbol{\mu}} \psi_i \, d\tau$$

where $\hat{\boldsymbol{\mu}}$ is the electric dipole operator. A transition is allowed if $\boldsymbol{\mu}_{fi} \neq 0$ and forbidden if $\boldsymbol{\mu}_{fi} = 0$.

4.2 Selection Rules for Different Spectroscopies

Selection rules specify which transitions are allowed based on symmetry and quantum mechanical considerations:

The complementarity of IR and Raman is a direct consequence of their different selection rules. For molecules with a center of symmetry, vibrations that are IR-active are Raman-inactive and vice versa (the rule of mutual exclusion).

"""
Analyze vibrational modes of CO2 and their IR/Raman activity based on selection rules.
"""
import numpy as np
import matplotlib.pyplot as plt

# CO2 is linear and centrosymmetric (point group D_infinity_h)
# It has 3N - 5 = 4 vibrational modes

vibrational_modes = {
    'Symmetric stretch': {
        'wavenumber': 1388,
        'symmetry': 'Sigma_g+',
        'dipole_change': False,  # No change in dipole moment
        'polarizability_change': True,  # Change in polarizability
        'IR_active': False,
        'Raman_active': True,
        'description': 'Both O atoms move away from C simultaneously'
    },
    'Asymmetric stretch': {
        'wavenumber': 2349,
        'symmetry': 'Sigma_u+',
        'dipole_change': True,  # Dipole moment changes
        'polarizability_change': False,
        'IR_active': True,
        'Raman_active': False,
        'description': 'One O moves toward C while other moves away'
    },
    'Bending (degenerate)': {
        'wavenumber': 667,
        'symmetry': 'Pi_u',
        'dipole_change': True,  # Dipole moment changes
        'polarizability_change': False,
        'IR_active': True,
        'Raman_active': False,
        'description': 'Molecule bends (2 degenerate modes)'
    }
}

print("CO2 Vibrational Modes and Selection Rules")
print("=" * 80)
print(f"{'Mode':<22} {'cm-1':<8} {'Symmetry':<12} {'IR Active':<12} {'Raman Active':<12}")
print("-" * 80)

for mode, props in vibrational_modes.items():
    ir = 'Yes' if props['IR_active'] else 'No'
    raman = 'Yes' if props['Raman_active'] else 'No'
    print(f"{mode:<22} {props['wavenumber']:<8} {props['symmetry']:<12} {ir:<12} {raman:<12}")

print("\nMutual Exclusion Rule:")
print("For centrosymmetric molecules like CO2, vibrations that are IR-active")
print("are Raman-inactive, and vice versa.")

# Simulate IR and Raman spectra
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

wavenumbers = np.linspace(400, 2600, 1000)

def gaussian_peak(x, center, amplitude, width):
    return amplitude * np.exp(-(x - center)**2 / (2 * width**2))

# IR spectrum
ir_spectrum = np.zeros_like(wavenumbers)
for mode, props in vibrational_modes.items():
    if props['IR_active']:
        ir_spectrum += gaussian_peak(wavenumbers, props['wavenumber'], 1.0, 30)

ax1.plot(wavenumbers, ir_spectrum, 'b-', linewidth=2)
ax1.fill_between(wavenumbers, ir_spectrum, alpha=0.3)
ax1.set_xlabel('Wavenumber (cm$^{-1}$)', fontsize=12)
ax1.set_ylabel('Absorbance', fontsize=12)
ax1.set_title('Simulated IR Spectrum of CO$_2$', fontsize=14, fontweight='bold')
ax1.invert_xaxis()

# Add labels
for mode, props in vibrational_modes.items():
    if props['IR_active']:
        ax1.annotate(mode, xy=(props['wavenumber'], 1.05), ha='center', fontsize=9)

# Raman spectrum
raman_spectrum = np.zeros_like(wavenumbers)
for mode, props in vibrational_modes.items():
    if props['Raman_active']:
        raman_spectrum += gaussian_peak(wavenumbers, props['wavenumber'], 1.0, 30)

ax2.plot(wavenumbers, raman_spectrum, 'r-', linewidth=2)
ax2.fill_between(wavenumbers, raman_spectrum, alpha=0.3, color='red')
ax2.set_xlabel('Raman Shift (cm$^{-1}$)', fontsize=12)
ax2.set_ylabel('Intensity', fontsize=12)
ax2.set_title('Simulated Raman Spectrum of CO$_2$', fontsize=14, fontweight='bold')
ax2.invert_xaxis()

for mode, props in vibrational_modes.items():
    if props['Raman_active']:
        ax2.annotate(mode, xy=(props['wavenumber'], 1.05), ha='center', fontsize=9)

plt.tight_layout()
plt.savefig('co2_spectra.png', dpi=150, bbox_inches='tight')
plt.show()

4.3 Fermi's Golden Rule

The transition probability per unit time from an initial state $|i\rangle$ to a final state $|f\rangle$ is given by Fermi's Golden Rule:

$$W_{i \to f} = \frac{2\pi}{\hbar} |\langle f | \hat{H}' | i \rangle|^2 \rho(E_f)$$

where $\hat{H}'$ is the perturbation Hamiltonian (describing the light-matter interaction) and $\rho(E_f)$ is the density of final states. For electric dipole transitions with light of electric field amplitude $E_0$:

$$W_{i \to f} \propto |\boldsymbol{\mu}_{fi}|^2 E_0^2$$

This shows that the absorption intensity is proportional to the square of the transition dipole moment.

5. Spectral Resolution and Instrumentation Basics

5.1 Spectral Resolution

Spectral resolution describes the ability to distinguish between two closely spaced spectral features. It is typically defined as:

$$R = \frac{\lambda}{\Delta\lambda} = \frac{\nu}{\Delta\nu} = \frac{\tilde{\nu}}{\Delta\tilde{\nu}}$$

where $\Delta\lambda$, $\Delta\nu$, or $\Delta\tilde{\nu}$ is the minimum resolvable separation. Higher resolution allows detection of finer spectral details but often requires longer measurement times.

5.2 Line Broadening Mechanisms

Spectral lines are never infinitely sharp. Several mechanisms contribute to line broadening:

• Natural linewidth (homogeneous): Due to the uncertainty principle, $\Delta E \cdot \Delta t \geq \hbar/2$, excited states with finite lifetimes have intrinsic energy uncertainty. This produces a Lorentzian lineshape.
• Doppler broadening (inhomogeneous): Due to thermal motion of molecules, different molecules see different frequencies. This produces a Gaussian lineshape.
• Collision broadening (homogeneous): Collisions interrupt the emission/absorption process, producing Lorentzian broadening.
• Instrumental broadening: Due to finite resolution of the spectrometer.

The convolution of Gaussian and Lorentzian lineshapes produces the Voigt profile, commonly observed in real spectra.

"""
Demonstrate Gaussian, Lorentzian, and Voigt lineshapes for spectral analysis.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import wofz  # Faddeeva function for Voigt profile

def gaussian(x, amplitude, center, sigma):
    """
    Gaussian lineshape.
    FWHM = 2.355 * sigma
    """
    return amplitude * np.exp(-(x - center)**2 / (2 * sigma**2))

def lorentzian(x, amplitude, center, gamma):
    """
    Lorentzian lineshape.
    FWHM = 2 * gamma
    """
    return amplitude * gamma**2 / ((x - center)**2 + gamma**2)

def voigt(x, amplitude, center, sigma, gamma):
    """
    Voigt profile - convolution of Gaussian and Lorentzian.
    Uses the Faddeeva function for efficient computation.
    """
    z = ((x - center) + 1j * gamma) / (sigma * np.sqrt(2))
    return amplitude * np.real(wofz(z)) / (sigma * np.sqrt(2 * np.pi))

# Create comparison plot
x = np.linspace(-5, 5, 500)
sigma = 1.0  # Gaussian width
gamma = 0.5  # Lorentzian width

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot lineshapes (normalized to unit height)
ax1 = axes[0]
y_gauss = gaussian(x, 1, 0, sigma)
y_lorentz = lorentzian(x, 1, 0, gamma) / lorentzian(0, 1, 0, gamma)  # Normalize
y_voigt = voigt(x, 1, 0, sigma, gamma)
y_voigt = y_voigt / np.max(y_voigt)  # Normalize

ax1.plot(x, y_gauss, 'b-', linewidth=2, label=f'Gaussian (sigma={sigma})')
ax1.plot(x, y_lorentz, 'r-', linewidth=2, label=f'Lorentzian (gamma={gamma})')
ax1.plot(x, y_voigt, 'g-', linewidth=2, label='Voigt')

ax1.set_xlabel('x - x0', fontsize=12)
ax1.set_ylabel('Normalized Intensity', fontsize=12)
ax1.set_title('Comparison of Lineshape Functions', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Log scale to show tail behavior
ax2 = axes[1]
ax2.semilogy(x, y_gauss, 'b-', linewidth=2, label='Gaussian')
ax2.semilogy(x, y_lorentz, 'r-', linewidth=2, label='Lorentzian')
ax2.semilogy(x, y_voigt, 'g-', linewidth=2, label='Voigt')

ax2.set_xlabel('x - x0', fontsize=12)
ax2.set_ylabel('Normalized Intensity (log scale)', fontsize=12)
ax2.set_title('Tail Behavior of Lineshapes', fontsize=14, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_ylim(1e-4, 2)

plt.tight_layout()
plt.savefig('lineshapes.png', dpi=150, bbox_inches='tight')
plt.show()

# FWHM calculations
print("\nFull Width at Half Maximum (FWHM):")
print("=" * 40)
print(f"Gaussian: FWHM = 2.355 * sigma = {2.355 * sigma:.3f}")
print(f"Lorentzian: FWHM = 2 * gamma = {2 * gamma:.3f}")
print(f"\nNote: Lorentzian has broader tails (slower decay)")
print("than Gaussian, which is important for baseline fitting.")

5.3 Basic Instrumentation Components

Most spectrometers share common components:

• Light source: Provides radiation across the spectral range of interest (tungsten lamp for visible, deuterium for UV, globar for IR)
• Wavelength selector: Isolates specific wavelengths (monochromator with diffraction grating, interferometer for FTIR, filters)
• Sample compartment: Holds the sample in the light path
• Detector: Converts light intensity to electrical signal (photomultiplier tube, CCD array, MCT detector for IR)

6. Practical Spectral Analysis with Python

"""
Complete spectral analysis pipeline: loading, preprocessing, peak finding, and fitting.
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks, savgol_filter
from scipy.optimize import curve_fit

# Generate synthetic UV-Vis spectrum with multiple peaks
np.random.seed(42)

def multi_gaussian(x, *params):
    """Sum of multiple Gaussian peaks."""
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        amplitude, center, sigma = params[i:i+3]
        y += amplitude * np.exp(-(x - center)**2 / (2 * sigma**2))
    return y

# Create wavelength axis
wavelengths = np.linspace(300, 700, 500)

# True spectrum parameters: (amplitude, center, sigma) for each peak
true_params = [
    0.8, 350, 15,   # Peak 1: UV absorption
    1.2, 420, 20,   # Peak 2: Main visible absorption
    0.5, 520, 25,   # Peak 3: Secondary absorption
    0.3, 600, 18,   # Peak 4: Weak absorption
]

# Generate spectrum with noise and baseline
true_spectrum = multi_gaussian(wavelengths, *true_params)
baseline = 0.1 + 0.0002 * (wavelengths - 500)**2  # Curved baseline
noise = np.random.normal(0, 0.02, len(wavelengths))
observed_spectrum = true_spectrum + baseline + noise

# Step 1: Baseline correction using polynomial fit
def baseline_correction(x, y, degree=2, n_iterations=3):
    """
    Iterative polynomial baseline correction.
    Excludes peaks progressively for better baseline estimation.
    """
    baseline_points = np.ones(len(y), dtype=bool)

    for _ in range(n_iterations):
        coeffs = np.polyfit(x[baseline_points], y[baseline_points], degree)
        baseline = np.polyval(coeffs, x)
        residual = y - baseline
        # Exclude points significantly above baseline
        threshold = np.std(residual[baseline_points])
        baseline_points = residual < 2 * threshold

    return baseline, y - baseline

estimated_baseline, corrected_spectrum = baseline_correction(
    wavelengths, observed_spectrum, degree=2
)

# Step 2: Smoothing (optional)
smoothed_spectrum = savgol_filter(corrected_spectrum, window_length=11, polyorder=3)

# Step 3: Peak detection
peaks, properties = find_peaks(smoothed_spectrum, height=0.1, prominence=0.05, width=3)

print("Automatic Peak Detection Results:")
print("=" * 50)
print(f"{'Peak':<8} {'Wavelength (nm)':<18} {'Height':<12}")
print("-" * 50)
for i, peak in enumerate(peaks):
    print(f"{i+1:<8} {wavelengths[peak]:<18.1f} {smoothed_spectrum[peak]:<12.3f}")

# Step 4: Multi-peak fitting
# Initial guess from detected peaks
initial_guess = []
for peak in peaks:
    initial_guess.extend([smoothed_spectrum[peak], wavelengths[peak], 15])

try:
    popt, pcov = curve_fit(multi_gaussian, wavelengths, smoothed_spectrum,
                           p0=initial_guess, maxfev=5000)
    fitted_spectrum = multi_gaussian(wavelengths, *popt)

    print("\nFitted Peak Parameters:")
    print("=" * 60)
    print(f"{'Peak':<8} {'Center (nm)':<15} {'Amplitude':<12} {'FWHM (nm)':<12}")
    print("-" * 60)
    for i in range(0, len(popt), 3):
        center = popt[i+1]
        amplitude = popt[i]
        fwhm = 2.355 * abs(popt[i+2])
        print(f"{i//3 + 1:<8} {center:<15.1f} {amplitude:<12.3f} {fwhm:<12.1f}")

    fit_success = True
except RuntimeError:
    print("Fitting did not converge")
    fit_success = False

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Original spectrum with baseline
ax1 = axes[0, 0]
ax1.plot(wavelengths, observed_spectrum, 'k-', linewidth=1, alpha=0.7, label='Observed')
ax1.plot(wavelengths, baseline, 'r--', linewidth=2, label='True baseline')
ax1.plot(wavelengths, estimated_baseline, 'b--', linewidth=2, label='Estimated baseline')
ax1.set_xlabel('Wavelength (nm)', fontsize=11)
ax1.set_ylabel('Absorbance', fontsize=11)
ax1.set_title('Step 1: Baseline Correction', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Corrected and smoothed spectrum
ax2 = axes[0, 1]
ax2.plot(wavelengths, corrected_spectrum, 'gray', linewidth=1, alpha=0.5, label='Corrected')
ax2.plot(wavelengths, smoothed_spectrum, 'b-', linewidth=2, label='Smoothed')
ax2.plot(wavelengths[peaks], smoothed_spectrum[peaks], 'ro', markersize=10, label='Detected peaks')
ax2.set_xlabel('Wavelength (nm)', fontsize=11)
ax2.set_ylabel('Absorbance', fontsize=11)
ax2.set_title('Step 2-3: Smoothing and Peak Detection', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Fitted spectrum
ax3 = axes[1, 0]
ax3.plot(wavelengths, smoothed_spectrum, 'ko', markersize=2, alpha=0.5, label='Data')
if fit_success:
    ax3.plot(wavelengths, fitted_spectrum, 'r-', linewidth=2, label='Fitted sum')
    # Plot individual peaks
    colors = ['#3498db', '#2ecc71', '#e74c3c', '#9b59b6', '#f39c12']
    for i in range(0, len(popt), 3):
        peak_params = popt[i:i+3]
        individual_peak = multi_gaussian(wavelengths, *peak_params)
        ax3.fill_between(wavelengths, individual_peak, alpha=0.3,
                        color=colors[i//3 % len(colors)],
                        label=f'Peak {i//3 + 1}')
ax3.set_xlabel('Wavelength (nm)', fontsize=11)
ax3.set_ylabel('Absorbance', fontsize=11)
ax3.set_title('Step 4: Multi-Peak Fitting', fontsize=12, fontweight='bold')
ax3.legend(loc='upper right')
ax3.grid(True, alpha=0.3)

# Residuals
ax4 = axes[1, 1]
if fit_success:
    residuals = smoothed_spectrum - fitted_spectrum
    ax4.plot(wavelengths, residuals, 'g-', linewidth=1)
    ax4.axhline(y=0, color='k', linestyle='--', linewidth=1)
    ax4.fill_between(wavelengths, residuals, alpha=0.3, color='green')
    ax4.set_xlabel('Wavelength (nm)', fontsize=11)
    ax4.set_ylabel('Residual', fontsize=11)
    ax4.set_title(f'Fitting Residuals (RMS = {np.sqrt(np.mean(residuals**2)):.4f})',
                  fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('spectral_analysis_pipeline.png', dpi=150, bbox_inches='tight')
plt.show()

Exercises

h = 6.626e-34  # J s
c = 2.998e8    # m/s
eV_to_J = 1.602e-19

# (a) Green laser at 532 nm
wl_a = 532e-9  # m
E_a = h * c / wl_a / eV_to_J  # eV
wn_a = 1 / (wl_a * 100)  # cm^-1
print(f"(a) 532 nm laser: {E_a:.3f} eV, {wn_a:.0f} cm^-1")

# (b) IR band at 1720 cm^-1
wn_b = 1720  # cm^-1
wl_b = 1 / wn_b  # cm
wl_b_um = wl_b * 1e4  # um
E_b = h * c * wn_b * 100 / eV_to_J  # eV
print(f"(b) 1720 cm^-1: {wl_b_um:.2f} um, {E_b:.4f} eV")

# (c) Al K-alpha at 1486.6 eV
E_c = 1486.6 * eV_to_J  # J
wl_c = h * c / E_c  # m
wl_c_nm = wl_c * 1e9  # nm
print(f"(c) 1486.6 eV: {wl_c_nm:.4f} nm")

# Output:
# (a) 532 nm laser: 2.331 eV, 18797 cm^-1
# (b) 1720 cm^-1: 5.81 um, 0.2132 eV
# (c) 1486.6 eV: 0.8340 nm

# Beer-Lambert Law: A = epsilon * c * l
A = 0.850
epsilon = 2340  # L mol^-1 cm^-1
l = 1.00  # cm
M = 158.03  # g/mol

c_mol_L = A / (epsilon * l)
c_mg_L = c_mol_L * M * 1000  # Convert to mg/L

print(f"Concentration: {c_mol_L:.4e} mol/L")
print(f"Concentration: {c_mg_L:.2f} mg/L")

# Output:
# Concentration: 3.6325e-04 mol/L
# Concentration: 57.41 mg/L

import numpy as np

# Vibrational frequency: nu = (1/2pi) * sqrt(k/mu)
# For the same force constant: nu1/nu2 = sqrt(mu2/mu1)

# Masses (amu)
m_H = 1.008
m_D = 2.014
m_C = 12.01

# Reduced masses for C-H and C-D
mu_CH = (m_C * m_H) / (m_C + m_H)
mu_CD = (m_C * m_D) / (m_C + m_D)

# Frequency ratio
ratio = np.sqrt(mu_CD / mu_CH)

# Predicted C-D frequency
nu_CH = 3019  # cm^-1
nu_CD = nu_CH / ratio

print(f"Reduced mass C-H: {mu_CH:.4f} amu")
print(f"Reduced mass C-D: {mu_CD:.4f} amu")
print(f"Frequency ratio (nu_CH/nu_CD): {ratio:.4f}")
print(f"Predicted C-D stretch: {nu_CD:.0f} cm^-1")
print(f"(Experimental value: ~2256 cm^-1)")

# Output:
# Reduced mass C-H: 0.9299 amu
# Reduced mass C-D: 1.7247 amu
# Frequency ratio: 1.3620
# Predicted C-D stretch: 2216 cm^-1

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial

def franck_condon_spectrum(E_elec, vib_freq_eV, S_factor, n_levels=10,
                           broadening=0.05, E_range=None):
    """
    Simulate absorption spectrum with vibrational fine structure.

    Parameters:
    -----------
    E_elec : float
        Electronic transition energy (eV)
    vib_freq_eV : float
        Vibrational frequency in eV
    S_factor : float
        Huang-Rhys factor (dimensionless displacement squared / 2)
    n_levels : int
        Number of vibrational levels to include
    broadening : float
        Gaussian broadening width (eV)
    E_range : tuple
        Energy range for spectrum (E_min, E_max)

    Returns:
    --------
    energies : array
        Energy axis
    spectrum : array
        Absorption intensity
    """
    # Calculate Franck-Condon factors for v=0 -> v' transitions
    fc_factors = []
    transition_energies = []

    for v_prime in range(n_levels):
        # FC factor for harmonic oscillators with v=0
        # |<0|v'>|^2 = exp(-S) * S^v' / v'!
        fc = np.exp(-S_factor) * (S_factor ** v_prime) / factorial(v_prime)
        fc_factors.append(fc)

        # Transition energy: E_elec + (v' + 1/2 - 1/2) * hv = E_elec + v' * hv
        E_trans = E_elec + v_prime * vib_freq_eV
        transition_energies.append(E_trans)

    fc_factors = np.array(fc_factors)
    transition_energies = np.array(transition_energies)

    # Create energy axis
    if E_range is None:
        E_min = E_elec - 0.5
        E_max = E_elec + (n_levels + 1) * vib_freq_eV
        E_range = (E_min, E_max)

    energies = np.linspace(E_range[0], E_range[1], 1000)

    # Build spectrum as sum of Gaussians
    spectrum = np.zeros_like(energies)
    for E_trans, fc in zip(transition_energies, fc_factors):
        spectrum += fc * np.exp(-(energies - E_trans)**2 / (2 * broadening**2))

    # Normalize
    spectrum /= np.max(spectrum)

    return energies, spectrum, transition_energies, fc_factors

# Example: Simulate spectrum for a typical organic chromophore
E_electronic = 3.0  # eV (visible region)
vib_frequency = 0.15  # eV (~1200 cm^-1, C-C stretch)
huang_rhys = 1.5  # Moderate displacement

energies, spectrum, E_trans, fc = franck_condon_spectrum(
    E_electronic, vib_frequency, huang_rhys, n_levels=8, broadening=0.04
)

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Spectrum
ax1.plot(energies, spectrum, 'b-', linewidth=2)
ax1.fill_between(energies, spectrum, alpha=0.3)
ax1.set_xlabel('Energy (eV)', fontsize=12)
ax1.set_ylabel('Normalized Absorbance', fontsize=12)
ax1.set_title('Simulated Absorption Spectrum with Vibronic Structure',
              fontsize=12, fontweight='bold')

# Mark individual transitions
for i, (E, f) in enumerate(zip(E_trans, fc)):
    if f > 0.01:
        ax1.axvline(x=E, color='red', linestyle='--', alpha=0.3)
        ax1.annotate(f"0-{i}", xy=(E, f), fontsize=8)

ax1.grid(True, alpha=0.3)

# Franck-Condon factors
ax2.bar(range(len(fc)), fc, color='#f093fb', edgecolor='#f5576c', linewidth=2)
ax2.set_xlabel("Final vibrational level (v')", fontsize=12)
ax2.set_ylabel('Franck-Condon Factor', fontsize=12)
ax2.set_title(f'Franck-Condon Factors (S = {huang_rhys})', fontsize=12, fontweight='bold')
ax2.set_xticks(range(len(fc)))
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('franck_condon_simulation.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"Maximum FC factor at v' = {np.argmax(fc)}")
print(f"Sum of FC factors: {np.sum(fc):.4f} (should be 1.0)")

# Generate test data
np.random.seed(123)
wavelengths = np.linspace(400, 600, 200)
true_spectrum = (
    0.8 * np.exp(-(wavelengths - 450)**2 / (2 * 20**2)) +
    1.0 * np.exp(-(wavelengths - 500)**2 / (2 * 25**2)) +
    0.5 * np.exp(-(wavelengths - 540)**2 / (2 * 15**2))
)
observed = true_spectrum + np.random.normal(0, 0.03, len(wavelengths))

import numpy as np
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
import matplotlib.pyplot as plt

def multi_gaussian(x, *params):
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        amp, center, sigma = params[i:i+3]
        y += amp * np.exp(-(x - center)**2 / (2 * sigma**2))
    return y

# Generate test data
np.random.seed(123)
wavelengths = np.linspace(400, 600, 200)
true_spectrum = (
    0.8 * np.exp(-(wavelengths - 450)**2 / (2 * 20**2)) +
    1.0 * np.exp(-(wavelengths - 500)**2 / (2 * 25**2)) +
    0.5 * np.exp(-(wavelengths - 540)**2 / (2 * 15**2))
)
observed = true_spectrum + np.random.normal(0, 0.03, len(wavelengths))

# Find peaks for initial guess
peaks, _ = find_peaks(observed, height=0.2, prominence=0.1)
initial_guess = []
for p in peaks:
    initial_guess.extend([observed[p], wavelengths[p], 20])

# Fit
popt, pcov = curve_fit(multi_gaussian, wavelengths, observed, p0=initial_guess)
perr = np.sqrt(np.diag(pcov))

# Results
print("Fitted Peak Parameters:")
print("=" * 70)
print(f"{'Peak':<8} {'Center (nm)':<15} {'Amplitude':<12} {'FWHM (nm)':<12}")
print("-" * 70)

for i in range(0, len(popt), 3):
    peak_num = i // 3 + 1
    center = popt[i + 1]
    amplitude = popt[i]
    fwhm = 2.355 * abs(popt[i + 2])
    print(f"{peak_num:<8} {center:<15.2f} {amplitude:<12.3f} {fwhm:<12.1f}")

# Plot
plt.figure(figsize=(10, 6))
plt.plot(wavelengths, observed, 'ko', markersize=3, alpha=0.5, label='Observed')
plt.plot(wavelengths, multi_gaussian(wavelengths, *popt), 'r-', linewidth=2, label='Fitted')

colors = ['blue', 'green', 'orange']
for i in range(0, len(popt), 3):
    individual = multi_gaussian(wavelengths, *popt[i:i+3])
    plt.fill_between(wavelengths, individual, alpha=0.3,
                     color=colors[i//3], label=f'Peak {i//3 + 1}')

plt.xlabel('Wavelength (nm)')
plt.ylabel('Absorbance')
plt.title('Multi-Peak Gaussian Fitting')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('exercise_6_solution.png', dpi=150)
plt.show()

Learning Objectives Check

After completing this chapter, you should be able to:

Fundamental Understanding

• Explain the wave-particle duality of light and convert between wavelength, frequency, wavenumber, and energy
• Describe the three fundamental light-matter interactions: absorption, emission, and scattering
• Identify which spectral region (UV, visible, IR, X-ray) is appropriate for different types of molecular transitions
• Explain the concept of quantized energy levels and their role in spectroscopy

Selection Rules

• Define the transition dipole moment and explain its role in determining transition probability
• State the key selection rules for electric dipole transitions, IR absorption, and Raman scattering
• Explain the complementarity of IR and Raman spectroscopy using the mutual exclusion rule
• Predict whether a given molecular vibration will be IR-active, Raman-active, or both

Practical Skills

• Perform unit conversions between spectroscopic quantities using Python
• Distinguish between Gaussian, Lorentzian, and Voigt lineshapes
• Apply baseline correction to spectroscopic data
• Detect peaks automatically and fit multi-component spectra

References

1. Hollas, J. M. (2004). Modern Spectroscopy (4th ed.). Wiley. - Comprehensive coverage of spectroscopic principles and selection rules.
2. Atkins, P., de Paula, J. (2014). Atkins' Physical Chemistry (10th ed.). Oxford University Press. - Quantum mechanical foundations of spectroscopy.
3. Banwell, C. N., McCash, E. M. (1994). Fundamentals of Molecular Spectroscopy (4th ed.). McGraw-Hill. - Clear introduction to molecular spectroscopy.
4. Harris, D. C., Bertolucci, M. D. (1989). Symmetry and Spectroscopy. Dover. - Detailed treatment of selection rules and group theory.
5. Stuart, B. (2004). Infrared Spectroscopy: Fundamentals and Applications. Wiley. - Practical guide to IR spectroscopy.
6. SciPy Documentation: Signal Processing (scipy.signal) and Optimization (scipy.optimize). https://docs.scipy.org/doc/scipy/

Next Steps

In this chapter, we have established the fundamental principles of spectroscopy: how light interacts with matter, the quantum mechanical basis for energy quantization, selection rules that govern transitions, and basic instrumentation concepts. We also developed practical Python skills for spectral analysis including unit conversion, lineshape functions, and peak fitting.

In Chapter 2, we will apply these fundamentals to UV-Vis spectroscopy, learning about electronic transitions, chromophores, the Beer-Lambert law for quantitative analysis, and practical applications in determining band gaps and concentrations.