Nonlinear Optics Second Harmonic Generation Shg Susceptibility Pockels Kerr Phase Matching

# Nonlinear Optics: SHG, Susceptibility Tensors, Pockels–Kerr Effects, and QPM Engineering

## 1. Introduction: Nonlinear Optical Response and High-Intensity Phenomena

When intense electromagnetic radiation propagates through a material, the optical response becomes nonlinear—higher-order polarization terms become significant. While linear optics (α¹) dominates for weak fields, nonlinear processes like second-harmonic generation (SHG), three-wave mixing, and four-wave mixing enable frequency conversion, optical parametric amplification, and the generation of new spectral components.

The nonlinear polarization is a power series in the electric field:

$$\mathbf{P}(t) = \epsilon_0 \left[\chi^{(1)} \mathbf{E}(t) + \chi^{(2)} \mathbf{E}(t)^2 + \chi^{(3)} \mathbf{E}(t)^3 + ...\right]$$

where:
- χ^(1): Linear susceptibility (first-order, responsible for refraction)
- χ^(2): Second-order nonlinear susceptibility (SHG, sum/difference frequency generation)
- χ^(3): Third-order susceptibility (Kerr effect, self-phase modulation, four-wave mixing)

The second-order susceptibility χ^(2) vanishes by symmetry in centrosymmetric media (isotropic liquids, cubic crystals) but is large in non-centrosymmetric crystals (GaAs, LiNbO₃, KDP). The third-order χ^(3) survives in all media, but is typically much smaller than χ^(2) where the latter exists.

The efficiency of nonlinear frequency conversion depends critically on phase matching: the wave vectors of interacting waves must satisfy k₁ + k₂ = k₃, preventing destructive interference over propagation distances.

## 2. Nonlinear Susceptibility Tensors and Symmetry

The nonlinear polarization component P_i^(n) depends on products of field components:

$$P_i^{(2)} = \epsilon_0 \chi_{ijk}^{(2)} E_j E_k$$
$$P_i^{(3)} = \epsilon_0 \chi_{ijkl}^{(3)} E_j E_k E_l$$

For second-order (SHG), the susceptibility tensor χ^(2)_{ijk} has 27 components, but crystal symmetry reduces this. For example:

  • Zincblende (GaAs): 3 independent components (χ_xyz = χ_xzy = χ_zxy)
  • Wurtzite (GaN): 5 independent components
  • Trigonal (LiNbO₃): 4 independent components

The third-order susceptibility χ^(3) has 81 components, further reduced by symmetry. For isotropic media:

$$\chi_{ijkl}^{(3)} = \chi_e \left(\delta_{ij} \delta_{kl} + \delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}\right) + \chi_o (\delta_{ij} \delta_{kl})$$

where χ_e (electronic) and χ_o (orientational) are scalar components.

The frequency dependence χ^(n)(ω₁, ω₂, ..., ω_n) is crucial. Near resonances, χ^(n) can be enhanced by factors of 10³–10⁶. For example, resonant SHG in semiconductors near bandgap shows enormous enhancement.

## 3. Second-Harmonic Generation: Coupled Wave Equations

Consider a pump wave at frequency ω₁ with amplitude A₁(z) generating second-harmonic at ω₂ = 2ω₁. The coupled-wave equations (in the slowly-varying envelope approximation) are:

$$\frac{\partial A_1}{\partial z} + \frac{\alpha_1}{2} A_1 = i \frac{\omega_1 d_{\text{eff}}}{c n_1} A_1^* A_2 e^{i\Delta k z}$$
$$\frac{\partial A_2}{\partial z} + \frac{\alpha_2}{2} A_2 = i \frac{\omega_2 d_{\text{eff}}}{2 c n_2} A_1^2 e^{i\Delta k z}$$

where:
- d_eff: Effective nonlinear coefficient (tensor contraction)
- Δk = k₂ - 2k₁ = (n₂ω₂ - 2n₁ω₁)/c: Phase mismatch
- α_i: Absorption coefficient at frequency ω_i
- n_i: Refractive index at ω_i

### Phase Matching (Δk = 0)

If Δk = 0 (perfect phase matching), the SHG intensity grows quadratically with propagation distance:

$$I_2(z) = I_2(0) + \left|\frac{\omega_2 d_{\text{eff}}}{2 c n_2}\right|^2 I_1^2 z^2$$

For quasi-phase-matching (QPM) with Δk ≠ 0, the conversion oscillates with coherence length:

$$L_c = \frac{\pi}{|\Delta k|} = \frac{\pi c}{|n_2 \omega_2 - 2n_1 \omega_1|}$$

Over this length, SHG efficiency rises to maximum, then oscillates. By periodically reversing the sign of χ^(2) (via domain engineering), QPM keeps the process in-phase constructively at all z, achieving efficiencies comparable to perfect phase matching but over wider parameter ranges.

### SHG Conversion Efficiency

The conversion efficiency for collinear SHG is:

$$\eta = \frac{I_2(L)}{I_1(0)} = \left(\frac{\omega_2 d_{\text{eff}} L}{2 c n_1 n_2}\right)^2 I_1(0) \text{sinc}^2\left(\frac{\Delta k L}{2}\right)$$

where sinc(x) = sin(x)/x represents the phase-matching factor. At perfect phase matching (Δk = 0):

$$\eta_{\text{max}} = \left(\frac{\omega_2 d_{\text{eff}} L}{2 c n_1 n_2}\right)^2 I_1(0)$$

For LiNbO₃ with d_eff ~ 10 pm/V, L ~ 1 cm, and I₁ ~ 1 MW/cm², the conversion can exceed 50%.

## 4. Type I and Type II Phase Matching in Birefringent Crystals

In birefringent crystals (n_e ≠ n_o), phase matching is achieved by choosing the polarization directions of pump and signal:

Type I Phase Matching: ω₁(o) + ω₁(o) → ω₂(e)
- Ordinary-polarized pump generates extraordinary-polarized harmonic
- Matching condition: n_o(ω₁) + n_o(ω₁) = n_e(ω₂)

Type II Phase Matching: ω₁(o) + ω₁(e) → ω₂(e) [or permutations]
- Mixed ordinary/extraordinary pump generates second harmonic
- Matching condition: n_o(ω₁) + n_e(ω₁) = n_e(ω₂)

For a given material and pump wavelength λ₁, phase matching can be achieved by temperature tuning (dn/dT), angle tuning (varying the propagation angle inside the crystal), or wavelength tuning (using different pump wavelengths).

The phase-matching bandwidth Δλ is inversely related to crystal length L:

$$\Delta \lambda \approx \frac{\lambda^2}{\pi L |dn/d\lambda|}$$

For L = 5 mm at 1.06 μm, Δλ ~ 0.1 nm in LiNbO₃. Longer crystals give narrower bandwidth but higher conversion.

## 5. Sum and Difference Frequency Generation

When two different frequencies ω₁ and ω₂ are incident:

$$\mathbf{P}^{(2)} \propto \chi^{(2)} [E_1 e^{-i\omega_1 t} + E_2 e^{-i\omega_2 t}]^2$$

This produces:
- Sum frequency: ω_s = ω₁ + ω₂ (SFG)
- Difference frequency: ω_d = |ω₁ - ω₂| (DFG)
- Second harmonics: 2ω₁, 2ω₂

The efficiency of SFG is optimized when all three waves are phase-matched:

$$k_s = k_1 + k_2 \quad \Rightarrow \quad \Delta k = 0$$

Difference frequency generation (DFG) is the inverse of SFG: a strong pump at ω_s and a weak signal at ω₁ generate an idler at ω_d = ω_s - ω₁. This is the basis of optical parametric oscillators (OPOs), which can tune the signal wavelength continuously.

## 6. Third-Order Nonlinear Effects: Kerr Effect and Self-Phase Modulation

The Kerr effect (χ^(3) contribution) induces an intensity-dependent refractive index:

$$n(I) = n_0 + n_2 I$$

where n₂ is the nonlinear refractive-index coefficient. Typical values: n₂ ~ 10⁻¹⁴ cm²/W in glass, ~10⁻¹¹ cm²/W in semiconductors.

Self-phase modulation (SPM) occurs when an intense pulse propagates through a Kerr medium:

$$\phi_{\text{SPM}}(t) = n_2 I(t, z) \omega_0 L / c$$

where L is the propagation distance. This nonlinear phase shift causes spectral broadening (self-steepening) and is crucial for supercontinuum generation and soliton formation in nonlinear fiber optics.

Cross-phase modulation (XPM): Two pulses at different wavelengths experience phase shifts from each other's intensity:

$$\phi_{\text{XPM}} = 2 n_2 I_{\text{other}} \omega_0 L / c$$

XPM is stronger than SPM by a factor of 2, enabling efficient wavelength conversion.

Four-wave mixing (FWM): When three waves ω₁, ω₂, ω₃ interact, a fourth wave at ω₄ = ω₁ + ω₂ - ω₃ is generated via χ^(3). In optical fiber:

$$\frac{\partial A_4}{\partial z} = i \gamma (\omega_4) A_1 A_2 A_3^* e^{i\Delta k_{\text{FWM}} z}$$

where γ is the FWM coefficient. FWM efficiency is maximum when degenerate FWM occurs: ω₁ = ω₃ (or ω₂ = ω₃), leading to parametric amplification.

## 7. Optical Parametric Oscillators and Amplifiers

An optical parametric oscillator (OPO) uses a nonlinear crystal inside an optical cavity. Three interacting waves:
- Pump (ω_p): supplied externally
- Signal (ω_s): cavity mode
- Idler (ω_i): complementary frequency (ω_p = ω_s + ω_i)

The cavity selects resonant modes at ω_s and ω_i, amplifying them via the pump. The signal wavelength tunes continuously by adjusting the phase-matching angle, temperature, or pump wavelength.

Tuning range: OPOs typically cover octaves of wavelength (e.g., 700–2500 nm by tuning LiNbO₃). Modern OPOs pumped by Nd:YAG or fiber lasers achieve efficiencies >50% and single-frequency output.

Optical parametric amplification (OPA): Without the cavity, a weak signal amplifies under pump-seeded nonlinear interaction. Gain coefficients can reach 1000× per cm with multi-pass geometries.

## 8. Integrated Photonics and Waveguide Nonlinearity

In integrated photonic waveguides (silicon nitride, LiNbO₃ on insulator), two effects enhance nonlinear efficiency:

1. Tight field confinement: Waveguide dimensions ~ λ create high optical intensities (I ~ GW/cm²) even for modest powers (P ~ 100 mW)
2. Long interaction length: Waveguides support SHG/FWM over centimeters or meters in very compact devices

Silicon photonics (Si waveguides, n ≈ 3.5) exhibit strong χ^(3) nonlinearity:
- Efficient FWM for parametric amplification
- SPM and XPM for soliton pulses
- Two-photon absorption (α_TPA ~ 0.5 cm/GW) limits efficiency at very high powers

Lithium niobate photonics (LiNbO₃ thin-film waveguides) combines χ^(2) and χ^(3):
- Ultra-high χ^(2) (d_eff ~ 100 pm/V) for SHG/DFG
- Electro-optic modulation (Pockels effect) with V_π ~ 2V
- Integrated OPOs on-chip
- Wavelength conversion across telecom C/L/O bands

Kerr nonlinearity in silicon nitride (Si₃N₄) enables:
- Broadband parametric amplification (>100 nm bandwidth)
- Soliton microcombs (optical frequency combs via Kerr nonlinearity in high-Q cavities)
- Integrated quantum photonic circuits

## 9. Numerical Solver: Coupled-Wave Equations for SHG

We implement a solver for SHG under phase mismatch:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp
from scipy.constants import c, pi

def shg_coupled_equations(y, z, params):
    """
    Coupled-wave equations for SHG:
    dA1/dz = i (ω₁ d_eff / (2cn₁)) |A1| A2 exp(i Δk z)
    dA2/dz = i (ω₂ d_eff / (2cn₂)) A1² exp(i Δk z)
    
    With depletion and phase mismatch.
    """
    A1_real, A1_imag, A2_real, A2_imag = y
    A1 = A1_real + 1j * A1_imag
    A2 = A2_real + 1j * A2_imag
    
    omega_1, omega_2, d_eff, n_1, n_2, delta_k, alpha_1, alpha_2 = params
    
    # Coupled-wave equations
    phase_mismatch = np.exp(1j * delta_k * z)
    
    dA1_dz = 1j * (omega_1 * d_eff / (2 * c * n_1)) * np.conj(A1) * A2 * phase_mismatch - (alpha_1/2) * A1
    dA2_dz = 1j * (omega_2 * d_eff / (2 * c * n_2)) * A1**2 * phase_mismatch - (alpha_2/2) * A2
    
    return [np.real(dA1_dz), np.imag(dA1_dz), np.real(dA2_dz), np.imag(dA2_dz)]

def solve_shg(L, lambda_1=1.064e-6, d_eff=10e-12, n_1=2.0, n_2=2.1, 
              alpha_1=0, alpha_2=0, I_pump_0=1e11):
    """
    Solve SHG coupled-wave equations.
    
    Parameters:
    - L: Crystal length (m)
    - lambda_1: Pump wavelength (m)
    - d_eff: Effective nonlinear coefficient (m/V)
    - n_1, n_2: Refractive indices at ω₁, 2ω₁
    - alpha_i: Absorption coefficients
    - I_pump_0: Input pump intensity (W/m²)
    """
    omega_1 = 2 * pi * c / lambda_1
    omega_2 = 2 * omega_1
    
    # Phase mismatch
    delta_k_0 = omega_2 / c * n_2 - 2 * omega_1 / c * n_1
    
    # Convert intensity to amplitude
    # I = (1/2) ε₀ c n |E|²  →  |E| = √(2I / (ε₀ c n))
    epsilon_0 = 8.854e-12  # F/m
    E_pump_0 = np.sqrt(2 * I_pump_0 / (epsilon_0 * c * n_1))
    
    # Initial conditions
    A1_0 = E_pump_0 + 0j
    A2_0 = 0 + 0j
    y0 = [np.real(A1_0), np.imag(A1_0), np.real(A2_0), np.imag(A2_0)]
    
    # Propagation
    z_array = np.linspace(0, L, 500)
    params = (omega_1, omega_2, d_eff, n_1, n_2, delta_k_0, alpha_1, alpha_2)
    
    solution = odeint(shg_coupled_equations, y0, z_array, args=(params,), rtol=1e-8)
    
    A1_array = solution[:, 0] + 1j * solution[:, 1]
    A2_array = solution[:, 2] + 1j * solution[:, 3]
    
    # Convert back to intensity
    I_1 = 0.5 * epsilon_0 * c * n_1 * np.abs(A1_array)**2
    I_2 = 0.5 * epsilon_0 * c * n_2 * np.abs(A2_array)**2
    
    return z_array, I_1, I_2

def phase_matching_bandwidth(L, lambda_1, dn_dlambda=0.01e-6):
    """
    Phase-matching bandwidth for type I SHG.
    Δλ ≈ λ² / (π L |dn/dλ|)
    """
    delta_lambda = lambda_1**2 / (pi * L * np.abs(dn_dlambda))
    return delta_lambda

def kerr_spm_spectrum(L, P_peak=1, n_0=1.5, n_2=1e-20, lambda_0=1.064e-6):
    """
    Self-phase modulation spectral broadening.
    Δφ_SPM = n_2 I L / c
    Spectral broadening: Δω ≈ |dφ/dt| ~ n_2 P_peak L / c
    """
    I_peak = P_peak / 1e-12  # Convert W to peak intensity (rough)
    phi_spm = n_2 * I_peak * L / c
    
    # Spectral broadening estimate
    delta_omega = n_2 * P_peak * L / c
    delta_lambda_kerr = lambda_0**2 / (2*pi*c) * delta_omega
    
    return phi_spm, delta_lambda_kerr

# Main execution
print("=" * 70)
print("NONLINEAR OPTICS: SHG & PHASE MATCHING ENGINEERING")
print("=" * 70)

# Parameters (LiNbO₃ crystal, type I SHG)
L_crystal = 0.01  # 1 cm
lambda_pump = 1.064e-6  # 1.064 μm (Nd:YAG)
lambda_shg = 0.532e-6  # 0.532 μm (green)
d_eff_lnb = 27e-12  # pm/V for LiNbO₃ type I
n_1 = 2.23  # Refractive index at 1.064 μm (o-ray)
n_2 = 2.29  # Refractive index at 0.532 μm (e-ray)
I_pump = 1e11  # W/m² (reasonable for focused 1 W in 100 μm² spot)

print(f"\nCrystal Parameters (LiNbO₃, Type I SHG):")
print(f"  Pump wavelength: {lambda_pump*1e6:.3f} μm ({c/lambda_pump/1e12:.0f} THz)")
print(f"  SHG wavelength: {lambda_shg*1e6:.3f} μm ({c/lambda_shg/1e12:.0f} THz)")
print(f"  d_eff: {d_eff_lnb*1e12:.1f} pm/V")
print(f"  Crystal length: {L_crystal*100:.1f} cm")
print(f"  Refractive indices: n_1 = {n_1:.3f}, n_2 = {n_2:.3f}")

# Phase mismatch calculation
omega_1 = 2*pi*c / lambda_pump
omega_2 = 2*omega_1
delta_k = omega_2/c*n_2 - 2*omega_1/c*n_1
L_coherence = pi / np.abs(delta_k)

print(f"\nPhase Matching Analysis:")
print(f"  Phase mismatch: Δk = {delta_k:.1e} m^-1")
print(f"  Coherence length: L_c = {L_coherence*1e6:.3f} μm")
print(f"  Crystal length / L_c = {L_crystal/L_coherence:.2f}")
print("  → Process is {} phase-matched".format(
    "nearly" if L_crystal < 2*L_coherence else "NOT"
))

# Solve SHG with phase mismatch
z_array, I_1, I_2 = solve_shg(L_crystal, lambda_pump, d_eff_lnb, n_1, n_2,
                               alpha_1=0, alpha_2=0, I_pump_0=I_pump)

conversion_eff = I_2[-1] / I_1[0] if I_1[0] > 0 else 0
print(f"\nSHG Conversion (with phase mismatch):")
print(f"  Input pump intensity: {I_pump:.2e} W/m²")
print(f"  Output SHG intensity: {I_2[-1]:.2e} W/m²")
print(f"  Conversion efficiency: {conversion_eff*100:.4f} %")

# Phase-matching bandwidth
pm_bandwidth = phase_matching_bandwidth(L_crystal, lambda_pump, dn_dlambda=0.015e-6)
print(f"\nPhase-Matching Bandwidth:")
print(f"  FWHM bandwidth: Δλ ≈ {pm_bandwidth*1e12:.3f} pm")

# Kerr effect estimate
phi_spm, delta_lambda_kerr = kerr_spm_spectrum(0.01, P_peak=1, n_2=1e-20, lambda_0=lambda_pump)
print(f"\nKerr (χ³) Nonlinearity (for comparison):")
print(f"  SPM phase shift: {phi_spm:.3f} rad/cm")
print(f"  Spectral broadening: {delta_lambda_kerr*1e12:.3f} pm")

# Plotting
fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# Panel 1: SHG intensity evolution
ax = axes[0, 0]
z_mm = z_array * 1000  # Convert to mm
I_1_mw = I_1 / 1e6  # Convert to MW/m²
I_2_mw = I_2 / 1e6
ax.plot(z_mm, I_1_mw, 'b-', linewidth=2.5, label='Pump (ω₁)', marker='o', markersize=3, markevery=50)
ax.plot(z_mm, I_2_mw, 'r-', linewidth=2.5, label='SHG (2ω₁)', marker='s', markersize=3, markevery=50)
ax.fill_between(z_mm, I_1_mw, alpha=0.2, color='blue')
ax.fill_between(z_mm, I_2_mw, alpha=0.2, color='red')
ax.set_xlabel('Crystal Depth z (mm)', fontsize=11)
ax.set_ylabel('Intensity (MW/m²)', fontsize=11)
ax.set_title('SHG Intensity Evolution with Phase Mismatch', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=10)

# Panel 2: Phase-matching sinc function
ax = axes[0, 1]
delta_k_range = np.linspace(-2/L_coherence, 2/L_coherence, 200)
eta_pm = np.sinc(delta_k_range * L_crystal / (2*pi))**2
ax.plot(delta_k_range/1e3, eta_pm*100, 'purple', linewidth=2.5)
ax.axvline(x=delta_k/1e3, color='green', linestyle='--', linewidth=2, label=f'Actual Δk = {delta_k:.1e} m⁻¹')
ax.fill_between(delta_k_range/1e3, eta_pm*100, alpha=0.2, color='purple')
ax.set_xlabel('Phase Mismatch Δk (10³ m⁻¹)', fontsize=11)
ax.set_ylabel('Conversion Efficiency η (%)', fontsize=11)
ax.set_title('SHG Efficiency vs Phase Mismatch', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=10)

# Panel 3: Nonlinear coefficient χ² vs wavelength (example)
ax = axes[1, 0]
lambda_range = np.linspace(0.4e-6, 2e-6, 100)
# Artificial d_eff vs wavelength (illustrative)
d_eff_range = 30 * np.exp(-((lambda_range - 1.064e-6)/0.3e-6)**2)
ax.plot(lambda_range*1e6, d_eff_range, 'b-', linewidth=2.5)
ax.fill_between(lambda_range*1e6, d_eff_range, alpha=0.2, color='blue')
ax.axvline(x=lambda_pump*1e6, color='red', linestyle='--', linewidth=2, label='Pump wavelength')
ax.set_xlabel('Wavelength (μm)', fontsize=11)
ax.set_ylabel('Effective Coefficient d_eff (pm/V)', fontsize=11)
ax.set_title('Wavelength-Dependent χ² Susceptibility', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=10)

# Panel 4: Conversion efficiency vs crystal length
ax = axes[1, 1]
L_range = np.linspace(0, 5e-2, 50)
efficiency_range = []

for L_test in L_range:
    if L_test == 0:
        efficiency_range.append(0)
    else:
        _, I1_test, I2_test = solve_shg(L_test, lambda_pump, d_eff_lnb, n_1, n_2)
        eff = I2_test[-1] / I1_test[0] if I1_test[0] > 0 else 0
        efficiency_range.append(eff * 100)

ax.semilogy(L_range*100, np.array(efficiency_range)+1e-6, 'g-', linewidth=2.5, marker='o', markersize=4)
ax.fill_between(L_range*100, np.array(efficiency_range)+1e-6, alpha=0.2, color='green')
ax.axvline(x=L_crystal*100, color='red', linestyle='--', linewidth=2, label=f'Current length: {L_crystal*100:.1f} cm')
ax.set_xlabel('Crystal Length (cm)', fontsize=11)
ax.set_ylabel('Conversion Efficiency (%)', fontsize=11)
ax.set_title('SHG Efficiency vs Crystal Length', fontsize=12, fontweight='bold')
ax.grid(True, alpha=0.3, which='both')
ax.legend(fontsize=10)

plt.tight_layout()
plt.savefig('nonlinear_optics_shg.png', dpi=150, bbox_inches='tight')
print(f"\nFigure saved: nonlinear_optics_shg.png")
plt.close()

## 10. Periodically Poled Structures and Quasi-Phase Matching

Periodically Poled LiNbO₃ (PPLN) uses ferroelectric domain engineering to reverse χ^(2) periodically. With period Λ, the domains are flipped such that:

$$\chi^{(2)}(z) = \chi^{(2)} \times \text{sign}[\cos(2\pi z/\Lambda)]$$

This generates a spatial Fourier component at the QPM wavelength:

$$\Delta k_{\text{eff}} = \Delta k_0 + \frac{2\pi m}{\Lambda}$$

where m is an integer. By choosing Λ appropriately, ΔK_eff can be tuned to zero, achieving phase matching over much broader bandwidths than birefringent phase matching.

PPLN enables:
- Broadband OPO: Tunable signal across octaves
- Multi-wavelength SHG: Different pump wavelengths phase-matched simultaneously
- Wavelength conversion: Arbitrary input-to-output frequency pairs

## 11. Applications in Integrated Photonics

Recent advances use thin-film LiNbO₃ waveguides (100–500 nm height) to achieve:
- Ultra-compact SHG (100 μm² footprint)
- On-chip OPOs with cavity Q > 10⁶
- Efficient frequency combs
- Quantum photonic gates (parametric entanglement)

## 12. Four-Wave Mixing and Parametric Amplification

When four waves ω₁, ω₂, ω₃, ω₄ interact via χ^(3), energy can transfer from pump to signal through idler generation:

$$\frac{\partial A_s}{\partial z} = i \gamma A_p A_i^* e^{i\Delta k z}$$

Parametric gain reaches maximum when degenerate (degenerate four-wave mixing, DFWM):

$$G = \frac{\sinh^2(g L)}{(\Delta k L / 2)^2 + \sinh^2(g L)}$$

where g = (2γP)/(ω_0 c) is the parametric gain coefficient. For perfect phase matching and high pump power, gain can exceed 1000× over a 100m fiber.

## 13. Supercontinuum Generation and Spectral Broadening

In nonlinear fibers with carefully engineered dispersion, several nonlinear effects combine:
- SPM: Spectral broadening
- FWM: New frequency generation
- Stimulated Raman scattering (molecular vibrations)
- Modulation instability: Pump depletion

Together, these generate a supercontinuum: broadband white-light spanning 400–2000 nm from a single narrow-linewidth pump laser.

## 14. Nonlinear Optical Switching and Computing

High-speed optical logic gates exploit nonlinear optical effects:
- All-optical switching: Kerr nonlinearity gate a pump beam with a control beam
- Optical soliton logic: Temporal pulses as information carriers
- Photonic neural networks: Integrated nonlinear photonic circuits

## 15. Future Directions: Metasurfaces and Topological Nonlinearity

Nonlinear metasurfaces: 2D arrays of plasmonic or dielectric nanoresonators with enhanced χ^(2) and χ^(3), enabling ultrathin frequency converters.

Topological photonics: Topologically-protected edge states in photonic crystals support robust nonlinear interactions immune to disorder.

---

Computational Notes: The coupled-wave equations are integrated using adaptive ODE solvers with high accuracy (rtol=10^-8) to capture phase mismatch oscillations. Intensity is computed from field amplitude via ε₀cnE². Phase-matching bandwidth uses standard nonlinear optics formulas validated against experimental data for LiNbO₃ and other crystals. Figures display both depletion (pump intensity decrease) and SHG build-up (harmonic generation), showing characteristic conversion curves with phase-mismatch oscillations when Δk ≠ 0.

## 16. Coupled-Wave Formulation and SHG Limit

For three-wave mixing under slowly varying envelopes, one requested convention writes

$$ \frac{dA_1}{dz} = i\frac{\omega_1d_{\mathrm{eff}}}{n_1c} A_3A_2^*e^{-i\Delta kz}, $$
$$ \frac{dA_2}{dz} = i\frac{\omega_2d_{\mathrm{eff}}}{n_2c} A_1^2e^{i\Delta kz}. $$

For degenerate second-harmonic generation, A₁ is the fundamental envelope and A₂ is the second harmonic; the fundamental depletion equation contains A₁* A₂. Envelope normalization determines additional factors of two and impedance. A solver must conserve the appropriate Manley–Rowe photon flux when absorption is absent.

With undepleted pump and uniform crystal length L, integration of the generated field produces the normalized phase-matching envelope

$$ \eta= \operatorname{sinc}^2\left(\frac{\Delta kL}{2}\right), \qquad \Delta k=k_2-2k_1. $$

Here sinc(x) = sin(x)/x. Absolute conversion additionally contains pump intensity, d_eff², L², refractive indices, frequency, mode overlap, and Fresnel or propagation losses.

## 17. Susceptibility Tensor, Pockels Effect, and Kerr Response

The nonlinear polarization expansion is

$$ P_i= \epsilon_0\left[ \chi_{ij}^{(1)}E_j+ \chi_{ijk}^{(2)}E_jE_k+ \chi_{ijkl}^{(3)}E_jE_kE_l+\cdots \right]. $$

Crystal symmetry removes many tensor components. Contracted notation maps the last two symmetric indices of χ^(2) to a d matrix. The effective coefficient is the projection of that tensor onto the interacting polarizations and propagation direction.

The linear electro-optic or Pockels response is commonly written for the impermeability tensor:

$$ \Delta\left(\frac{1}{n^2}\right)_i = \sum_j r_{ij}E_j. $$

The r_ij matrix depends on crystal point group and axis convention. The optical Kerr effect produces an intensity-dependent refractive index

$$ \Delta n=n_2I. $$

The symbol n₂ here is a nonlinear refractive-index coefficient, not the second-harmonic refractive index. Stating units avoids ambiguity.

## 18. Quasi-Phase Matching and Duty Cycle

Periodic inversion of the nonlinear coefficient supplies a reciprocal lattice vector. First-order quasi-phase matching uses

$$ \Lambda=\frac{2\pi}{\Delta k}. $$

For a binary grating with duty fraction D, the magnitude of the first Fourier coefficient scales as 2 sin(πD)/π. It is maximal at 50% duty cycle. Fabrication error, domain-wall roughness, chirped periods, temperature gradients, and wavelength dispersion broaden or reduce conversion.

Higher-order QPM uses odd spatial harmonics with smaller effective nonlinearity. QPM allows use of a large tensor coefficient without relying on birefringent phase matching, but it does not remove group-velocity mismatch or absorption.

## 19. Python Coupled-Wave, Phase-Mismatch, and Duty-Cycle Solver

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

L = 20e-3
kappa = 120.0
A10 = 1.0+0j
A20 = 0.0+0j

def propagate(delta_k, qpm=False, duty=0.5):
    def rhs(z, y):
        A1 = y[0]+1j*y[1]
        A2 = y[2]+1j*y[3]
        if qpm:
            period = 2*np.pi/max(abs(delta_k), 1e-12)
            phase = np.mod(z, period)/period
            sign = 1.0 if phase < duty else -1.0
        else:
            sign = 1.0
        dA1 = 1j*kappa*sign*np.conjugate(A1)*A2*np.exp(-1j*delta_k*z)
        dA2 = 0.5j*kappa*sign*A1**2*np.exp(1j*delta_k*z)
        return [dA1.real, dA1.imag, dA2.real, dA2.imag]
    z = np.linspace(0, L, 1400)
    sol = solve_ivp(rhs, (0, L), [A10.real, A10.imag, 0, 0],
                    t_eval=z, rtol=2e-8, atol=1e-10)
    A1 = sol.y[0]+1j*sol.y[1]
    A2 = sol.y[2]+1j*sol.y[3]
    return z, A1, A2

z, A1, A2 = propagate(0.0)
mismatch_L = np.linspace(-14*np.pi, 14*np.pi, 500)
undepleted = np.sinc((mismatch_L/2)/np.pi)**2

powers = []
for value in np.linspace(-8*np.pi, 8*np.pi, 160):
    _, _, out = propagate(value/L)
    powers.append(abs(out[-1])**2)

duty = np.linspace(0.05, 0.95, 300)
qpm_fourier = (2*np.sin(np.pi*duty)/np.pi)**2

fig, axes = plt.subplots(2, 2, figsize=(12, 9))
axes[0, 0].plot(z*1e3, abs(A1)**2, label="fundamental")
axes[0, 0].plot(z*1e3, abs(A2)**2, label="second harmonic")
axes[0, 0].set(xlabel="Propagation z (mm)", ylabel="Normalized power")
axes[0, 0].legend()

axes[0, 1].plot(mismatch_L, undepleted)
axes[0, 1].set(xlabel="Phase mismatch ΔkL", ylabel="sinc² efficiency")

axes[1, 0].plot(np.linspace(-8*np.pi, 8*np.pi, 160), powers)
axes[1, 0].set(xlabel="ΔkL", ylabel="Depleted-pump SH power")

axes[1, 1].plot(duty, qpm_fourier)
axes[1, 1].axvline(0.5, color="black", ls="--")
axes[1, 1].set(xlabel="Poling duty cycle", ylabel="Relative first-order QPM power")
for ax in axes.flat:
    ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()

The equations use normalized amplitudes to emphasize energy exchange. Quantitative power prediction requires a declared amplitude normalization, beam area or modal overlap, pulse duration, group velocities, absorption, focusing, and temperature-dependent Sellmeier dispersion.

## 20. Experimental Validation

A phase-matching result should state wavelength, polarization, crystal axes, temperature, length, aperture, focusing parameter, pulse bandwidth, duty cycle, and whether efficiency is internal or external. Measuring the full temperature or wavelength acceptance curve is stronger than quoting a single peak.

Pockels and Kerr measurements must separate instantaneous electronic response from thermal, photorefractive, free-carrier, and cascaded-χ² effects. Their time scales and power dependencies provide useful controls.

Take nonlinear optics second harmonic generation SHG susceptibility Pockels Kerr phase matching further

Ask the copilot about this term, or have our engineers assess it against your process.