Third Harmonic Generation in 1D

This example demonstrates third harmonic generation (THG) in a one-dimensional Kerr nonlinear medium. A plane wave at frequency ω propagates through a material with χ⁽³⁾ nonlinearity and generates light at the third harmonic frequency 3ω.

Reference: This example is based on the Meep Third Harmonic Generation tutorial.

Physics Background

Kerr Nonlinearity

In a Kerr medium, the polarization includes a term proportional to the cube of the electric field:

$$ P = \epsilon_0 (\chi^1 E + \chi^3 E^3) $$

When a monochromatic wave at frequency ω enters such a medium, the E³ term generates new frequency components, including the third harmonic at 3ω.

Third Harmonic Efficiency

The third harmonic power scales as:

Quadratic with χ⁽³⁾ at low intensities
Saturates at high intensities due to phase mismatch

Simulation Overview

Setup the 1D Kerr Medium

Create a 1D simulation cell with a nonlinear medium characterized by χ⁽³⁾.

Launch Gaussian Pulse

A broadband Gaussian pulse centered at frequency ω is launched into the medium.

Monitor Transmitted Flux

Flux monitors capture both the fundamental (ω) and third harmonic (3ω) frequencies.

Analyze Conversion Efficiency

Compare transmitted power at different χ⁽³⁾ values to understand the nonlinear scaling.

Configuration Parameters

Parameter	Value	Description
`sz`	100	Cell size in z direction
`fcen`	1/3.0	Center frequency of source
`dpml`	1.0	PML thickness
`resolution`	25	Grid resolution
`chi3`	10⁻³ to 1	Kerr susceptibility range

Complete Code

"""
Third Harmonic Generation in 1D with OptixLog Integration

1D simulation of a plane wave propagating through a Kerr medium
and generating the third-harmonic frequency component.
"""

import os
import meep as mp
import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
import numpy as np
from typing import Tuple, List, Union
from optixlog import Optixlog


def third_harmonic_generation(
    k: float,
    amp: float = 1.0,
    nfreq: int = 10,
    flux_spectrum: bool = True,
    run=None,
    step_offset: int = 0
) -> Union[Tuple[List[float], List[float]], Tuple[float, float]]:
    """Computes the transmission spectrum of a plane wave propagating
       through a Kerr medium.

    Args:
        k: strength of Kerr susceptibility.
        amp: amplitude of the incident planewave.
        nfreq: number of frequencies in flux spectrum.
        flux_spectrum: compute the flux spectrum over broad bandwidth (True) or
                      just the two harmonic frequencies ω and 3ω (False).
        run: OptixLog run object for logging.
        step_offset: offset for step numbering in logging.

    Returns:
        The frequencies and transmitted flux over the flux spectrum or
        the transmitted flux at the harmonic frequencies ω and 3ω.
    """
    sz = 100  # size of cell in z direction
    fcen = 1 / 3.0  # center frequency of source
    df = fcen / 20.0  # frequency width of source
    dpml = 1.0  # PML thickness

    # Explicitly 1D simulation - faster execution
    dimensions = 1
    cell = mp.Vector3(0, 0, sz)
    pml_layers = [mp.PML(dpml)]
    resolution = 25

    # Kerr medium with χ(3) nonlinearity
    default_material = mp.Medium(index=1, chi3=k)

    sources = [
        mp.Source(
            mp.GaussianSource(fcen, fwidth=df),
            component=mp.Ex,
            center=mp.Vector3(0, 0, -0.5 * sz + dpml),
            amplitude=amp,
        )
    ]

    # Frequency range for flux calculation
    fmin = fcen / 2.0
    fmax = fcen * 4

    sim = mp.Simulation(
        cell_size=cell,
        sources=sources,
        boundary_layers=pml_layers,
        default_material=default_material,
        resolution=resolution,
        dimensions=dimensions,
    )

    mon_pt = mp.Vector3(0, 0, 0.5 * sz - dpml - 0.5)

    if flux_spectrum:
        trans = sim.add_flux(
            0.5 * (fmin + fmax),
            fmax - fmin,
            nfreq,
            mp.FluxRegion(mon_pt),
        )
    else:
        trans1 = sim.add_flux(fcen, 0, 1, mp.FluxRegion(mon_pt))
        trans3 = sim.add_flux(3 * fcen, 0, 1, mp.FluxRegion(mon_pt))

    sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mon_pt, 1e-6))

    if flux_spectrum:
        freqs = mp.get_flux_freqs(trans)
        trans_flux = mp.get_fluxes(trans)
        return freqs, trans_flux
    else:
        flux1 = mp.get_fluxes(trans1)[0]
        flux3 = mp.get_fluxes(trans3)[0]
        return flux1, flux3


def main():
    """Main simulation function with OptixLog integration"""
    
    # Initialize OptixLog client
    api_key = os.getenv("OPTIX_API_KEY", "your_api_key_here")
    
    client = Optixlog(api_key=api_key)
    project = client.project(name="MeepExamples", create_if_not_exists=True)
    
    run = project.run(
        name="third_harmonic_generation",
        config={
            "simulation_type": "third_harmonic_generation",
            "description": "1D plane wave through Kerr medium generating third harmonic",
            "framework": "meep",
            "cell_size_z": 100,
            "fcen": 1/3.0,
            "dpml": 1.0,
            "resolution": 25
        }
    )
    
    print(f"🚀 Starting Third Harmonic Generation simulation")
    print(f"📊 Logging to OptixLog run: {run.run_id}")
    
    # Log initial parameters
    run.log(step=0, 
            simulation_started=True,
            cell_size_z=100,
            center_frequency=1/3.0,
            pml_thickness=1.0,
            resolution=25)

    # Part 1: Plot transmitted power spectrum for several values of χ(3)
    print("\n📈 Part 1: Computing transmission spectra for varying χ(3)...")
    nfreq = 400
    logk = range(-3, 1)
    tflux = np.zeros((nfreq, len(logk)))
    
    for i, lk in enumerate(logk):
        chi3_value = 10**lk
        print(f"   Computing χ(3) = {chi3_value}")
        freqs, tflux[:, i] = third_harmonic_generation(
            chi3_value, nfreq=nfreq, flux_spectrum=True, run=run
        )
        
        # Log progress
        run.log(step=i+1, 
                chi3_value=chi3_value,
                log_chi3=lk,
                max_flux=float(np.max(tflux[:, i])),
                min_flux=float(np.min(tflux[:, i])))

    # Create and log transmission spectrum plot
    fig, ax = plt.subplots(figsize=(10, 7))
    ax.semilogy(freqs, tflux[:, 0], "bo-", markersize=2, label=r"$\chi^{(3)}$=0.001")
    ax.semilogy(freqs, tflux[:, 1], "ro-", markersize=2, label=r"$\chi^{(3)}$=0.01")
    ax.semilogy(freqs, tflux[:, 2], "go-", markersize=2, label=r"$\chi^{(3)}$=0.1")
    ax.semilogy(freqs, tflux[:, 3], "co-", markersize=2, label=r"$\chi^{(3)}$=1")
    ax.set_xlabel("Frequency", fontsize=12)
    ax.set_ylabel("Transmitted Power (a.u.)", fontsize=12)
    ax.set_title("Third Harmonic Generation: Transmitted Power Spectrum", fontsize=14)
    ax.set_xlim(0.2, 1.2)
    ax.set_ylim(1e-15, 1e2)
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    
    run.log_matplotlib("transmitted_power_vs_frequency", fig)
    plt.close(fig)

    # Part 2: Plot transmittance vs. χ(3) for frequencies ω and 3ω
    print("\n📈 Part 2: Computing transmittance vs χ(3) at harmonic frequencies...")
    logk = np.arange(-6.0, 0.2, 0.2)
    first_order = np.zeros(len(logk))
    third_order = np.zeros(len(logk))
    
    for i, lk in enumerate(logk):
        chi3_value = 10**lk
        first_order[i], third_order[i] = third_harmonic_generation(
            chi3_value, flux_spectrum=False, run=run
        )

    input_flux = first_order[0]
    
    # Log harmonic analysis results
    run.log(step=len(logk)+5,
            input_flux=float(input_flux),
            max_first_order=float(np.max(first_order)),
            max_third_order=float(np.max(third_order)),
            conversion_efficiency=float(np.max(third_order)/input_flux))

    # Create and log transmittance vs χ(3) plot
    fig, ax = plt.subplots(figsize=(10, 7))
    ax.loglog(10**logk, first_order / input_flux, "ro-", markersize=4, label=r"$\omega$ (fundamental)")
    ax.loglog(10**logk, third_order / input_flux, "bo-", markersize=4, label=r"$3\omega$ (third harmonic)")
    ax.loglog(10**logk, (10**logk) ** 2, "k--", linewidth=2, label="Quadratic reference")
    ax.set_xlabel(r"$\chi^{(3)}$ (Kerr susceptibility)", fontsize=12)
    ax.set_ylabel("Transmission / Incident Power", fontsize=12)
    ax.set_title("Third Harmonic Generation Efficiency", fontsize=14)
    ax.legend(fontsize=10)
    ax.grid(True, which="both", alpha=0.3)
    plt.tight_layout()
    
    run.log_matplotlib("transmittance_vs_chi3", fig)
    plt.close(fig)

    print(f"\n✅ Simulation complete!")
    print(f"🔗 View results at: https://optixlog.com")


if __name__ == "__main__":
    main()

Key Concepts

1D Simulation Setup

For 1D simulations, set dimensions=1 and use a cell with only one non-zero dimension:

dimensions = 1
cell = mp.Vector3(0, 0, sz)  # Only z-dimension

Kerr Medium Definition

Define a nonlinear medium using the chi3 parameter:

default_material = mp.Medium(index=1, chi3=k)

Flux Spectrum Analysis

Use add_flux with multiple frequency points to capture the spectrum:

trans = sim.add_flux(
    center_freq,
    freq_width,
    nfreq,  # Number of frequency points
    mp.FluxRegion(mon_pt),
)

Expected Results

Transmission Spectrum

The transmission spectrum shows:

Fundamental peak at ω ≈ 0.33
Third harmonic peak at 3ω ≈ 1.0
Peak height increases with larger χ⁽³⁾

Scaling Behavior

At low χ⁽³⁾ values:

Third harmonic power scales as (χ⁽³⁾)²
Matches the quadratic reference line

At high χ⁽³⁾ values:

Deviation from quadratic scaling
Saturation effects become important

OptixLog Integration

This example logs:

Metric	Description
`chi3_value`	Current χ⁽³⁾ being simulated
`max_flux`	Maximum flux in spectrum
`input_flux`	Input power for normalization
`conversion_efficiency`	Third harmonic / input ratio

Logged Plots

transmitted_power_vs_frequency — Spectrum for different χ⁽³⁾ values
transmittance_vs_chi3 — Scaling of harmonic efficiency

Running the Example

# Set your API key
export OPTIX_API_KEY="your_api_key_here"

# Run the simulation
python third-harmonic-1d.py

Third Harmonic Generation (1D)