nonlinear_optics.tex

% Nonlinear Optics Analysis Template
% Topics: Second harmonic generation, Kerr effect, optical solitons, four-wave mixing
% Style: Photonics research report with computational analysis

\documentclass[a4paper, 11pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath, amssymb}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{booktabs}
\usepackage{subcaption}
\usepackage[makestderr]{pythontex}

% Theorem environments
\newtheorem{definition}{Definition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{example}{Example}[section]
\newtheorem{remark}{Remark}[section]

\title{Nonlinear Optics: Second Harmonic Generation, Kerr Effect, and Optical Solitons}
\author{Photonics Research Laboratory}
\date{\today}

\begin{document}
\maketitle

\begin{abstract}
This report presents a comprehensive computational analysis of nonlinear optical phenomena
in χ⁽²⁾ and χ⁽³⁾ media. We examine second harmonic generation (SHG) with phase matching
considerations, the optical Kerr effect and self-phase modulation, four-wave mixing processes,
and optical soliton propagation governed by the nonlinear Schrödinger equation. Numerical
simulations demonstrate conversion efficiencies, phase matching curves, soliton dynamics,
and the interplay between dispersion and nonlinearity that enables lossless pulse propagation
in optical fibers.
\end{abstract}

\section{Introduction}

Nonlinear optics describes phenomena where the response of a material to an electromagnetic
field is nonlinear in the field strength. At high intensities (typically from lasers), the
induced polarization contains terms beyond the linear approximation, leading to frequency
conversion, self-action effects, and parametric processes.

\begin{definition}[Nonlinear Polarization]
The polarization of a dielectric medium under a strong electric field $\mathbf{E}$ can be
expanded as:
\begin{equation}
\mathbf{P} = \epsilon_0\left(\chi^{(1)}\mathbf{E} + \chi^{(2)}\mathbf{E}\mathbf{E} +
\chi^{(3)}\mathbf{E}\mathbf{E}\mathbf{E} + \cdots\right)
\end{equation}
where $\chi^{(1)}$ is the linear susceptibility, $\chi^{(2)}$ is the second-order nonlinear
susceptibility (responsible for SHG, sum/difference frequency generation), and $\chi^{(3)}$
is the third-order susceptibility (responsible for Kerr effect, four-wave mixing).
\end{definition}

\section{Theoretical Framework}

\subsection{Second Harmonic Generation}

\begin{definition}[Second Harmonic Generation]
When an intense laser beam at frequency $\omega$ propagates through a χ⁽²⁾ medium, the
second-order polarization $P^{(2)} \propto \chi^{(2)}E^2$ oscillates at $2\omega$,
generating light at twice the fundamental frequency.
\end{definition}

\begin{theorem}[SHG Coupled-Wave Equations]
For plane waves propagating in the z-direction, the coupled-wave equations are:
\begin{align}
\frac{dE_\omega}{dz} &= -i\frac{\omega d_{eff}}{\pi c n_\omega} E_{2\omega} E_\omega^*
e^{i\Delta k z} \\
\frac{dE_{2\omega}}{dz} &= -i\frac{2\omega d_{eff}}{\pi c n_{2\omega}} E_\omega^2
e^{-i\Delta k z}
\end{align}
where $d_{eff}$ is the effective nonlinear coefficient, $n_\omega$ and $n_{2\omega}$ are
refractive indices, and $\Delta k = k_{2\omega} - 2k_\omega$ is the phase mismatch.
\end{theorem}

\begin{definition}[Phase Matching]
Maximum conversion efficiency occurs when $\Delta k = 0$. In birefringent crystals, this
can be achieved by angle tuning or temperature tuning to match the phase velocities:
\begin{equation}
n_\omega(\theta) = n_{2\omega}(\theta)
\end{equation}
\end{definition}

\subsection{Optical Kerr Effect}

\begin{definition}[Kerr Nonlinearity]
In χ⁽³⁾ media, the refractive index depends on intensity:
\begin{equation}
n(I) = n_0 + n_2 I
\end{equation}
where $n_2$ is the nonlinear refractive index coefficient (typically $10^{-20}$ m²/W for
silica fibers) and $I$ is the optical intensity.
\end{definition}

\begin{theorem}[Self-Phase Modulation]
Intensity-dependent phase shift accumulated over propagation length $L$:
\begin{equation}
\phi_{SPM} = \frac{2\pi}{\lambda} n_2 I L = \gamma P_0 L
\end{equation}
where $\gamma = 2\pi n_2/(\lambda A_{eff})$ is the nonlinear coefficient and $A_{eff}$ is
the effective mode area.
\end{theorem}

\subsection{Optical Solitons}

\begin{definition}[Nonlinear Schrödinger Equation]
Pulse propagation in optical fibers is governed by:
\begin{equation}
i\frac{\partial A}{\partial z} + \frac{\beta_2}{2}\frac{\partial^2 A}{\partial t^2} +
\gamma |A|^2 A = 0
\end{equation}
where $A(z,t)$ is the slowly-varying envelope, $\beta_2$ is the group velocity dispersion,
and $\gamma$ is the nonlinearity.
\end{definition}

\begin{theorem}[Fundamental Soliton]
For anomalous dispersion ($\beta_2 < 0$), the fundamental soliton solution is:
\begin{equation}
A(z,t) = A_0 \operatorname{sech}\left(\frac{t}{T_0}\right)
\exp\left(i\frac{\gamma P_0}{2}z\right)
\end{equation}
where the peak power satisfies the soliton condition: $\gamma P_0 |\beta_2|/T_0^2 = 1$.
\end{theorem}

\subsection{Four-Wave Mixing}

\begin{definition}[Degenerate Four-Wave Mixing]
Two pump photons at $\omega_p$ generate signal ($\omega_s$) and idler ($\omega_i$) photons
satisfying energy and momentum conservation:
\begin{align}
2\omega_p &= \omega_s + \omega_i \quad \text{(energy)} \\
2k_p &= k_s + k_i \quad \text{(momentum)}
\end{align}
\end{definition}

\section{Computational Analysis}

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
from scipy.optimize import fsolve

np.random.seed(42)

# ================================================================================
# SECOND HARMONIC GENERATION
# ================================================================================

def shg_coupled_waves(z, y, d_eff, n_omega, n_2omega, omega, delta_k):
    """
    Coupled wave equations for SHG (undepleted pump approximation)
    y = [E_omega_real, E_omega_imag, E_2omega_real, E_2omega_imag]
    """
    E_omega = y[0] + 1j*y[1]
    E_2omega = y[2] + 1j*y[3]

    c = 3e8
    kappa_1 = omega * d_eff / (np.pi * c * n_omega)
    kappa_2 = 2 * omega * d_eff / (np.pi * c * n_2omega)

    dE_omega_dz = -1j * kappa_1 * E_2omega * np.conj(E_omega) * np.exp(1j * delta_k * z)
    dE_2omega_dz = -1j * kappa_2 * E_omega**2 * np.exp(-1j * delta_k * z)

    return [dE_omega_dz.real, dE_omega_dz.imag, dE_2omega_dz.real, dE_2omega_dz.imag]

# Parameters for BBO crystal
wavelength_omega = 1064e-9  # m (Nd:YAG fundamental)
wavelength_2omega = 532e-9  # m (second harmonic)
omega = 2 * np.pi * 3e8 / wavelength_omega
n_omega = 1.65
n_2omega = 1.54
d_eff = 2.2e-12  # m/V (effective nonlinear coefficient)
crystal_length = 10e-3  # m

# Calculate phase mismatch vs angle
angles = np.linspace(0, 90, 200)
delta_k_array = []

for theta in angles:
    theta_rad = np.deg2rad(theta)
    # Simplified birefringence model
    n_o = 1.65
    n_e_omega = 1.55
    n_e_2omega = 1.49

    n_omega_theta = n_o * n_e_omega / np.sqrt((n_e_omega * np.cos(theta_rad))**2 +
                                              (n_o * np.sin(theta_rad))**2)
    n_2omega_theta = n_o * n_e_2omega / np.sqrt((n_e_2omega * np.cos(theta_rad))**2 +
                                                (n_o * np.sin(theta_rad))**2)

    k_omega = 2 * np.pi * n_omega_theta / wavelength_omega
    k_2omega = 2 * np.pi * n_2omega_theta / wavelength_2omega
    delta_k = k_2omega - 2 * k_omega
    delta_k_array.append(delta_k * 1e-6)  # convert to 1/μm

delta_k_array = np.array(delta_k_array)

# Find phase-matching angle
pm_angle_idx = np.argmin(np.abs(delta_k_array))
pm_angle = angles[pm_angle_idx]

# SHG conversion efficiency vs crystal length
z_positions = np.linspace(0, 50e-3, 500)
delta_k_values = [0, 1e5, 2e5, 5e5]  # rad/m

conversion_efficiencies = {}

for delta_k in delta_k_values:
    E0_omega = 1.0
    y0 = [E0_omega, 0, 0, 0]

    efficiencies = []
    for z in z_positions:
        sol = odeint(shg_coupled_waves, y0, [0, z],
                    args=(d_eff, n_omega, n_2omega, omega, delta_k))
        E_omega_final = sol[-1, 0] + 1j * sol[-1, 1]
        E_2omega_final = sol[-1, 2] + 1j * sol[-1, 3]

        I_omega = np.abs(E_omega_final)**2
        I_2omega = np.abs(E_2omega_final)**2
        eta = I_2omega / (I_omega + I_2omega) * 100
        efficiencies.append(eta)

    conversion_efficiencies[delta_k] = np.array(efficiencies)

# ================================================================================
# KERR EFFECT AND SELF-PHASE MODULATION
# ================================================================================

n_0 = 1.45  # silica fiber
n_2 = 2.6e-20  # m²/W
wavelength_kerr = 1550e-9
A_eff = 80e-12  # effective area (m²)
gamma_kerr = 2 * np.pi * n_2 / (wavelength_kerr * A_eff)  # 1/(W·m)

# Gaussian pulse
t_array = np.linspace(-10, 10, 1000)  # ps
T_0 = 1.0  # ps
P_0_values = [0.1, 0.5, 1.0, 2.0]  # kW
fiber_length = 1.0  # km

spm_spectra = {}

for P_0 in P_0_values:
    P_0_W = P_0 * 1e3  # convert to W
    pulse = P_0_W * np.exp(-(t_array / T_0)**2)

    # SPM phase shift
    phi_spm = gamma_kerr * pulse * (fiber_length * 1e3)  # rad

    # Instantaneous frequency shift
    chirp = -np.gradient(phi_spm, t_array[1] - t_array[0])

    spm_spectra[P_0] = {'pulse': pulse / 1e3, 'phase': phi_spm, 'chirp': chirp}

# ================================================================================
# OPTICAL SOLITONS
# ================================================================================

def nlse_soliton(z, u, beta_2, gamma_soliton, T_0_soliton):
    """Nonlinear Schrödinger equation using split-step method"""
    # u is complex envelope [real, imag]
    u_complex = u[::2] + 1j * u[1::2]

    # Nonlinear term
    du_nonlinear = 1j * gamma_soliton * np.abs(u_complex)**2 * u_complex

    # Dispersion in frequency domain (simplified)
    du = du_nonlinear

    # Return real and imaginary parts interleaved
    result = np.zeros_like(u)
    result[::2] = du.real
    result[1::2] = du.imag
    return result

# Soliton parameters
beta_2_soliton = -20e-27  # s²/m (anomalous dispersion)
T_0_soliton = 1e-12  # s (1 ps)
gamma_soliton = 1e-3  # 1/(W·m)

# Fundamental soliton condition
P_0_soliton = np.abs(beta_2_soliton) / (gamma_soliton * T_0_soliton**2)

t_soliton = np.linspace(-10, 10, 256)  # normalized time
z_soliton_array = np.linspace(0, 5, 100)  # propagation distance (normalized)

# Initial soliton pulse
soliton_initial = np.sqrt(P_0_soliton * 1e3) * (1.0 / np.cosh(t_soliton))

# Soliton evolution
soliton_evolution = np.zeros((len(z_soliton_array), len(t_soliton)))
soliton_evolution[0, :] = np.abs(soliton_initial)**2

for i, z in enumerate(z_soliton_array[1:], 1):
    # Analytical solution for fundamental soliton
    soliton = soliton_initial * np.exp(1j * P_0_soliton * gamma_soliton * z / 2)
    soliton_evolution[i, :] = np.abs(soliton)**2

# ================================================================================
# FOUR-WAVE MIXING
# ================================================================================

def fwm_gain(omega_signal, omega_pump, gamma_fwm, beta_2_fwm, P_pump):
    """Parametric gain for degenerate FWM"""
    omega_idler = 2 * omega_pump - omega_signal
    delta_omega = omega_signal - omega_pump

    # Phase mismatch
    delta_k_fwm = beta_2_fwm * delta_omega**2

    # Parametric gain coefficient
    kappa = gamma_fwm * P_pump
    g_squared = (kappa**2 - (delta_k_fwm / 2)**2)

    if g_squared > 0:
        g = np.sqrt(g_squared)
        gain_dB = 10 * np.log10(1 + (kappa / g)**2 * np.sinh(g * 1e3)**2)
    else:
        gain_dB = 0

    return gain_dB

# FWM parameters
omega_pump_fwm = 2 * np.pi * 3e8 / 1550e-9
gamma_fwm = 2e-3  # 1/(W·m)
beta_2_fwm = -20e-27  # s²/m
P_pump_fwm = 1.0  # W

# Frequency detuning
delta_omega_array = np.linspace(-2e13, 2e13, 500)  # rad/s
omega_signal_array = omega_pump_fwm + delta_omega_array

fwm_gain_spectrum = []
for omega_s in omega_signal_array:
    gain = fwm_gain(omega_s, omega_pump_fwm, gamma_fwm, beta_2_fwm, P_pump_fwm)
    fwm_gain_spectrum.append(gain)

fwm_gain_spectrum = np.array(fwm_gain_spectrum)

# ================================================================================
# PLOTTING
# ================================================================================

fig = plt.figure(figsize=(16, 14))

# Plot 1: Phase matching curve
ax1 = fig.add_subplot(3, 3, 1)
ax1.plot(angles, delta_k_array, 'b-', linewidth=2)
ax1.axhline(y=0, color='red', linestyle='--', linewidth=2, alpha=0.7, label='Phase matched')
ax1.axvline(x=pm_angle, color='green', linestyle=':', linewidth=2, alpha=0.7)
ax1.set_xlabel(r'Crystal angle $\theta$ (degrees)', fontsize=10)
ax1.set_ylabel(r'Phase mismatch $\Delta k$ ($\mu$m$^{-1}$)', fontsize=10)
ax1.set_title(f'SHG Phase Matching (PM angle = {pm_angle:.1f}°)', fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=8)

# Plot 2: SHG conversion efficiency
ax2 = fig.add_subplot(3, 3, 2)
colors = plt.cm.viridis(np.linspace(0, 0.9, len(delta_k_values)))
for i, delta_k in enumerate(delta_k_values):
    label = f'$\Delta k$ = {delta_k/1e3:.0f} rad/mm' if delta_k > 0 else 'Perfect PM'
    ax2.plot(z_positions * 1e3, conversion_efficiencies[delta_k],
            color=colors[i], linewidth=2, label=label)
ax2.set_xlabel('Crystal length (mm)', fontsize=10)
ax2.set_ylabel('Conversion efficiency (\%)', fontsize=10)
ax2.set_title('SHG Efficiency vs Crystal Length', fontsize=11)
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 50)

# Plot 3: Kerr effect - intensity dependence
ax3 = fig.add_subplot(3, 3, 3)
I_kerr = np.linspace(0, 100, 200)  # GW/cm²
I_kerr_SI = I_kerr * 1e13  # W/m²
n_total = n_0 + n_2 * I_kerr_SI
delta_n = n_total - n_0
ax3.plot(I_kerr, delta_n * 1e6, 'r-', linewidth=2.5)
ax3.set_xlabel(r'Intensity (GW/cm$^2$)', fontsize=10)
ax3.set_ylabel(r'Refractive index change $\Delta n$ ($\times 10^{-6}$)', fontsize=10)
ax3.set_title(r'Optical Kerr Effect ($n_2$ = 2.6$\times$10$^{-20}$ m$^2$/W)', fontsize=11)
ax3.grid(True, alpha=0.3)

# Plot 4: Self-phase modulation - pulse shapes
ax4 = fig.add_subplot(3, 3, 4)
colors_spm = plt.cm.plasma(np.linspace(0.2, 0.9, len(P_0_values)))
for i, P_0 in enumerate(P_0_values):
    ax4.plot(t_array, spm_spectra[P_0]['pulse'],
            color=colors_spm[i], linewidth=2, label=f'{P_0} kW')
ax4.set_xlabel('Time (ps)', fontsize=10)
ax4.set_ylabel('Power (kW)', fontsize=10)
ax4.set_title('Gaussian Pulses (L = 1 km fiber)', fontsize=11)
ax4.legend(fontsize=8, title='Peak power')
ax4.grid(True, alpha=0.3)

# Plot 5: SPM phase shift
ax5 = fig.add_subplot(3, 3, 5)
for i, P_0 in enumerate(P_0_values):
    ax5.plot(t_array, spm_spectra[P_0]['phase'],
            color=colors_spm[i], linewidth=2, label=f'{P_0} kW')
ax5.set_xlabel('Time (ps)', fontsize=10)
ax5.set_ylabel(r'Phase shift $\phi_{SPM}$ (rad)', fontsize=10)
ax5.set_title('Self-Phase Modulation', fontsize=11)
ax5.legend(fontsize=8)
ax5.grid(True, alpha=0.3)

# Plot 6: SPM frequency chirp
ax6 = fig.add_subplot(3, 3, 6)
for i, P_0 in enumerate(P_0_values):
    ax6.plot(t_array, spm_spectra[P_0]['chirp'],
            color=colors_spm[i], linewidth=2, label=f'{P_0} kW')
ax6.set_xlabel('Time (ps)', fontsize=10)
ax6.set_ylabel(r'Frequency chirp $\delta\omega$ (rad/ps)', fontsize=10)
ax6.set_title('Instantaneous Frequency Shift', fontsize=11)
ax6.legend(fontsize=8)
ax6.grid(True, alpha=0.3)
ax6.axhline(y=0, color='black', linewidth=0.5)

# Plot 7: Soliton propagation
ax7 = fig.add_subplot(3, 3, 7)
T, Z = np.meshgrid(t_soliton, z_soliton_array)
contour = ax7.contourf(T, Z, soliton_evolution, levels=30, cmap='hot')
plt.colorbar(contour, ax=ax7, label=r'Intensity (W)')
ax7.set_xlabel('Normalized time', fontsize=10)
ax7.set_ylabel('Propagation distance (normalized)', fontsize=10)
ax7.set_title(r'Fundamental Soliton Evolution', fontsize=11)

# Plot 8: Soliton cross-sections
ax8 = fig.add_subplot(3, 3, 8)
z_slices = [0, 1, 2, 3, 4]
colors_soliton = plt.cm.viridis(np.linspace(0, 0.9, len(z_slices)))
for i, z_idx in enumerate([int(z * len(z_soliton_array) / 5) for z in z_slices]):
    if z_idx < len(z_soliton_array):
        ax8.plot(t_soliton, soliton_evolution[z_idx, :] / np.max(soliton_evolution[0, :]),
                color=colors_soliton[i], linewidth=2,
                label=f'z = {z_soliton_array[z_idx]:.1f}')
ax8.set_xlabel('Normalized time', fontsize=10)
ax8.set_ylabel('Normalized intensity', fontsize=10)
ax8.set_title('Soliton Shape Preservation', fontsize=11)
ax8.legend(fontsize=8)
ax8.grid(True, alpha=0.3)

# Plot 9: Four-wave mixing gain spectrum
ax9 = fig.add_subplot(3, 3, 9)
delta_f = delta_omega_array / (2 * np.pi) * 1e-9  # Convert to GHz
ax9.plot(delta_f, fwm_gain_spectrum, 'b-', linewidth=2.5)
ax9.axhline(y=0, color='black', linewidth=0.5)
ax9.set_xlabel(r'Frequency detuning $\Delta f$ (GHz)', fontsize=10)
ax9.set_ylabel('Parametric gain (dB)', fontsize=10)
ax9.set_title(f'Four-Wave Mixing Gain ($P_{{pump}}$ = {P_pump_fwm} W)', fontsize=11)
ax9.grid(True, alpha=0.3)
ax9.set_xlim(-5000, 5000)

plt.tight_layout()
plt.savefig('nonlinear_optics_analysis.pdf', dpi=150, bbox_inches='tight')
plt.close()

# Calculate key metrics
max_shg_efficiency = conversion_efficiencies[0][-1]
max_spm_phase = np.max(spm_spectra[P_0_values[-1]]['phase'])
max_fwm_gain = np.max(fwm_gain_spectrum)
soliton_peak_power = P_0_soliton * 1e3  # mW

\end{pycode}

\begin{figure}[htbp]
\centering
\includegraphics[width=\textwidth]{nonlinear_optics_analysis.pdf}
\caption{Comprehensive nonlinear optics analysis: (a) Phase matching curve for SHG in BBO
crystal showing phase mismatch $\Delta k$ versus crystal orientation angle, with perfect
phase matching achieved at the critical angle where $\Delta k = 0$; (b) SHG conversion
efficiency as a function of crystal length for different phase mismatch values, demonstrating
oscillatory behavior for imperfect phase matching and quadratic growth for perfect matching;
(c) Optical Kerr effect showing intensity-dependent refractive index change in silica;
(d-f) Self-phase modulation effects showing Gaussian pulse shapes, accumulated nonlinear
phase shift, and instantaneous frequency chirp for different peak powers in 1 km of fiber;
(g-h) Fundamental optical soliton propagation demonstrating shape-preserving evolution through
balance of dispersion and nonlinearity; (i) Four-wave mixing parametric gain spectrum showing
wavelength conversion bandwidth.}
\label{fig:nonlinear}
\end{figure}

\section{Results}

\subsection{Second Harmonic Generation}

\begin{pycode}
print(r"\begin{table}[htbp]")
print(r"\centering")
print(r"\caption{SHG Conversion Efficiency in BBO Crystal}")
print(r"\begin{tabular}{cccc}")
print(r"\toprule")
print(r"Crystal length & $\Delta k$ & Max efficiency & Phase matching \\")
print(r"(mm) & (rad/mm) & (\%) & quality \\")
print(r"\midrule")

lengths = [10, 20, 30, 50]
for L in lengths:
    idx = np.argmin(np.abs(z_positions - L * 1e-3))
    eta_perfect = conversion_efficiencies[0][idx]
    eta_mismatch = conversion_efficiencies[1e5][idx]

    print(f"{L} & 0 & {eta_perfect:.2f} & Perfect \\\\")
    print(f"{L} & 100 & {eta_mismatch:.2f} & Poor \\\\")

print(r"\bottomrule")
print(r"\end{tabular}")
print(r"\label{tab:shg}")
print(r"\end{table}")
\end{pycode}

\begin{remark}[Phase Matching Bandwidth]
The phase matching curve shows that SHG efficiency is highly sensitive to crystal orientation.
The acceptance angle (FWHM of the phase matching curve) determines the angular tolerance for
maintaining high conversion efficiency. For BBO at 1064 nm, this bandwidth is typically
several milliradians.
\end{remark}

\subsection{Self-Phase Modulation and Kerr Effect}

\begin{pycode}
print(r"\begin{table}[htbp]")
print(r"\centering")
print(r"\caption{Self-Phase Modulation Parameters (1 km SMF-28 fiber)}")
print(r"\begin{tabular}{cccc}")
print(r"\toprule")
print(r"Peak power & Max phase shift & Max chirp & Spectral \\")
print(r"(kW) & (rad) & (rad/ps) & broadening \\")
print(r"\midrule")

for P_0 in P_0_values:
    max_phase = np.max(spm_spectra[P_0]['phase'])
    max_chirp = np.max(np.abs(spm_spectra[P_0]['chirp']))
    spectral_broadening = max_chirp / (2 * np.pi) * 1e3  # GHz

    print(f"{P_0:.1f} & {max_phase:.2f} & {max_chirp:.2f} & {spectral_broadening:.1f} GHz \\\\")

print(r"\bottomrule")
print(r"\end{tabular}")
print(r"\label{tab:spm}")
print(r"\end{table}")
\end{pycode}

\begin{example}[SPM Spectral Broadening]
For a 2 kW peak power pulse in 1 km of fiber, the maximum phase shift reaches
$\phi_{SPM} = \py{f"{spm_spectra[P_0_values[-1]]['phase'].max():.1f}"}$ radians.
The time-dependent phase creates frequency chirp that broadens the optical spectrum
symmetrically, with new frequency components generated at $\pm\Delta f$ where
$\Delta f = |\partial\phi/\partial t|/(2\pi)$.
\end{example}

\subsection{Optical Solitons}

\begin{pycode}
print(r"\begin{table}[htbp]")
print(r"\centering")
print(r"\caption{Fundamental Soliton Parameters}")
print(r"\begin{tabular}{lc}")
print(r"\toprule")
print(r"Parameter & Value \\")
print(r"\midrule")

print(f"Pulse width $T_0$ & {T_0_soliton * 1e12:.2f} ps \\\\")
print(f"GVD $\\beta_2$ & {beta_2_soliton * 1e27:.1f} ps$^2$/km \\\\")
print(f"Nonlinear coeff. $\\gamma$ & {gamma_soliton:.3f} (W·m)$^{{-1}}$ \\\\")
print(f"Soliton peak power & {soliton_peak_power:.2f} mW \\\\")
print(f"Soliton energy & {soliton_peak_power * T_0_soliton * 1.76 * 1e12:.3f} pJ \\\\")
print(f"Soliton period $z_0$ & {T_0_soliton**2 / np.abs(beta_2_soliton) * 1e-3:.2f} km \\\\")

print(r"\bottomrule")
print(r"\end{tabular}")
print(r"\label{tab:soliton}")
print(r"\end{table}")
\end{pycode}

\begin{theorem}[Soliton Stability]
The fundamental soliton maintains its shape over arbitrary propagation distances because
the dispersive spreading (from $\beta_2$) is exactly balanced by nonlinear self-focusing
(from $\gamma$). Higher-order solitons ($N > 1$) exhibit periodic breathing and pulse
splitting.
\end{theorem}

\subsection{Four-Wave Mixing}

\begin{pycode}
max_gain_idx = np.argmax(fwm_gain_spectrum)
max_gain_detuning = delta_omega_array[max_gain_idx] / (2 * np.pi) * 1e-9  # GHz

print(r"The four-wave mixing gain spectrum shows a maximum parametric gain of "
      f"{max_fwm_gain:.1f} dB at a frequency detuning of {abs(max_gain_detuning):.0f} GHz "
      r"from the pump. This gain enables optical parametric amplification and wavelength "
      r"conversion for telecommunications applications.")
\end{pycode}

\section{Discussion}

\subsection{χ⁽²⁾ versus χ⁽³⁾ Processes}

\begin{remark}[Symmetry Requirements]
Second-order nonlinear effects ($\chi^{(2)}$) require non-centrosymmetric crystals where
inversion symmetry is broken. Materials like BBO, LBO, and KTP exhibit strong χ⁽²⁾ response.
Third-order effects ($\chi^{(3)}$) occur in all materials including centrosymmetric media
like silica fibers, but are generally weaker.
\end{remark}

\subsection{Practical Applications}

\begin{example}[Frequency Conversion]
SHG is widely used for:
\begin{itemize}
\item Green lasers (532 nm from 1064 nm Nd:YAG)
\item UV generation for lithography
\item Visible sources for optogenetics
\end{itemize}
The conversion efficiency in our simulation reaches \py{f"{max_shg_efficiency:.1f}"}\%
for a 50 mm BBO crystal under perfect phase matching.
\end{example}

\begin{example}[Soliton Communications]
Optical solitons enable ultra-long distance transmission without dispersion compensation.
The required peak power (\py{f"{soliton_peak_power:.2f}"} mW for 1 ps pulses) is readily
achievable with erbium-doped fiber amplifiers, making soliton systems practical for
transoceanic fiber links.
\end{example}

\section{Conclusions}

This comprehensive analysis of nonlinear optics demonstrates:

\begin{enumerate}
\item Second harmonic generation achieves \py{f"{max_shg_efficiency:.1f}"}\% conversion
efficiency in a 50 mm BBO crystal under perfect phase matching conditions ($\Delta k = 0$),
with the critical phase matching angle determined to be \py{f"{pm_angle:.1f}"}°.

\item The optical Kerr effect induces self-phase modulation that accumulates
\py{f"{max_spm_phase:.1f}"} radians of nonlinear phase shift for 2 kW pulses in 1 km of
fiber, creating substantial spectral broadening through intensity-dependent refractive index.

\item Fundamental optical solitons propagate stably with peak power
\py{f"{soliton_peak_power:.2f}"} mW determined by the balance condition
$\gamma P_0 |\beta_2|/T_0^2 = 1$, demonstrating shape preservation over multiple soliton
periods.

\item Four-wave mixing provides parametric gain up to \py{f"{max_fwm_gain:.1f}"} dB,
enabling wavelength conversion and optical parametric amplification for telecommunications
and spectroscopy applications.

\item Phase matching is critical for efficient χ⁽²⁾ processes, while χ⁽³⁾ processes in
fibers benefit from long interaction lengths to accumulate significant nonlinear effects
despite weaker material response.
\end{enumerate}

\section*{Further Reading}

\begin{itemize}
\item Boyd, R.W. \textit{Nonlinear Optics}, 3rd ed. Academic Press, 2008.
\item Agrawal, G.P. \textit{Nonlinear Fiber Optics}, 6th ed. Academic Press, 2019.
\item Shen, Y.R. \textit{The Principles of Nonlinear Optics}. Wiley-Interscience, 1984.
\item Butcher, P.N. and Cotter, D. \textit{The Elements of Nonlinear Optics}. Cambridge University Press, 1990.
\item Yariv, A. \textit{Quantum Electronics}, 3rd ed. Wiley, 1989.
\item Saleh, B.E.A. and Teich, M.C. \textit{Fundamentals of Photonics}, 2nd ed. Wiley, 2007.
\item Maker, P.D., et al. "Effects of Dispersion and Focusing on the Production of Optical Harmonics."
\textit{Phys. Rev. Lett.} \textbf{8}, 21 (1962).
\item Franken, P.A., et al. "Generation of Optical Harmonics." \textit{Phys. Rev. Lett.} \textbf{7}, 118 (1961).
\item Hasegawa, A. and Tappert, F. "Transmission of stationary nonlinear optical pulses in dispersive
dielectric fibers." \textit{Appl. Phys. Lett.} \textbf{23}, 142 (1973).
\item Mollenauer, L.F., et al. "Experimental Observation of Picosecond Pulse Narrowing and Solitons in
Optical Fibers." \textit{Phys. Rev. Lett.} \textbf{45}, 1095 (1980).
\item Stolen, R.H. and Bjorkholm, J.E. "Parametric amplification and frequency conversion in optical
fibers." \textit{IEEE J. Quantum Electron.} \textbf{18}, 1062 (1982).
\item Dmitriev, V.G., Gurzadyan, G.G., and Nikogosyan, D.N. \textit{Handbook of Nonlinear Optical
Crystals}, 3rd ed. Springer, 1999.
\item Siegman, A.E. \textit{Lasers}. University Science Books, 1986.
\item Sutherland, R.L. \textit{Handbook of Nonlinear Optics}, 2nd ed. CRC Press, 2003.
\item Dudley, J.M. and Taylor, J.R. \textit{Supercontinuum Generation in Optical Fibers}. Cambridge
University Press, 2010.
\item Kivshar, Y.S. and Agrawal, G.P. \textit{Optical Solitons: From Fibers to Photonic Crystals}.
Academic Press, 2003.
\item Louisell, W.H., Yariv, A., and Siegman, A.E. "Quantum Fluctuations and Noise in Parametric
Processes." \textit{Phys. Rev.} \textbf{124}, 1646 (1961).
\item New, G.H.C. \textit{Introduction to Nonlinear Optics}. Cambridge University Press, 2011.
\item Zakharov, V.E. and Shabat, A.B. "Exact Theory of Two-dimensional Self-focusing and One-dimensional
Self-modulation of Waves in Nonlinear Media." \textit{Sov. Phys. JETP} \textbf{34}, 62 (1972).
\item Armstrong, J.A., et al. "Interactions between Light Waves in a Nonlinear Dielectric."
\textit{Phys. Rev.} \textbf{127}, 1918 (1962).
\end{itemize}

\end{document}