Advanced API

FFTContext{N, T<:Complex, F, I}

Thread-safe container for FFT plans. Plans can be shared across threads since fftw_execute is thread-safe. Create once and reuse.

Fields

• resolution: Grid resolution tuple
• fft_plan: Pre-computed forward FFT plan
• ifft_plan: Pre-computed inverse FFT plan

MatrixFreeWorkspace{N, T<:Complex}

Thread-local workspace arrays for matrix-free operations. Each thread should have its own workspace instance.

Fields

• work_real: Workspace array in real space
• work_fourier: Workspace array in Fourier space

MatrixFreeOperator{D<:Dimension, W<:WaveType, T<:Complex}

Matrix-free representation of the Hamiltonian operator H = RHS⁻¹ * LHS.

Fields

• solver: Reference to the solver (contains geometry, basis, material_arrays)
• k: Current wave vector
• ctx: FFT context (shared, thread-safe for execution)
• workspace: Work arrays (thread-local)

Thread Safety

The FFT plans in ctx can be safely shared across threads. However, each thread must have its own workspace to avoid data races. For parallel k-point loops:

ctx = FFTContext(solver)  # Create once
Threads.@threads for k in k_points
    ws = MatrixFreeWorkspace(ctx)  # One per thread
    op = MatrixFreeOperator(solver, k, ctx, ws)
    # ... use op ...
end

Grid Resolution Constraint

Important: The solver's grid resolution must be large enough to contain all plane wave indices. The required resolution in each dimension is:

resolution >= 2 * cutoff + 1

For example, with cutoff=7 (default), the resolution must be at least 15 in each dimension (e.g., 16×16 for 2D, 16×16×16 for 3D).

If the resolution is too small, plane wave coefficients outside the grid range will be silently dropped during FFT operations, leading to incorrect results.

apply_lhs!(y, op::MatrixFreeOperator, x)

Apply LHS operator: y = LHS * x For TE: LHS = Kx * ε⁻¹ * Kx + Ky * ε⁻¹ * Ky

apply_lhs!(y, op::MatrixFreeOperator{Dim1}, x)

Apply LHS operator for 1D: y = K * ε⁻¹ * K * x

apply_rhs!(y, op::MatrixFreeOperator, x)

Apply RHS operator: y = RHS * x

to_linear_map(op::MatrixFreeOperator)

Convert MatrixFreeOperator to a LinearMap for use with iterative solvers. Returns LHS as a LinearMap.

to_linear_map_rhs(op::MatrixFreeOperator)

Convert RHS to a LinearMap.

EffectiveHamiltonian{TL, TR, TF}

Effective Hamiltonian operator H = RHS⁻¹ * LHS for use with Krylov.jl.

Transforms the generalized eigenvalue problem: LHS * ψ = ω² * RHS * ψ into the standard form: H * ψ = ω² * ψ

For Green's function computation: G(z) = (z·I - H)⁻¹ = (z·I - RHS⁻¹·LHS)⁻¹

Fields

• LHS: Left-hand side matrix
• RHS: Right-hand side matrix
• RHS_factorization: LU factorization of RHS for efficient solves
• n: Problem dimension
• tmp: Temporary vector for intermediate computations

Usage with Krylov.jl

Use NegatedOperator(H) to get -H, then solve (σI + A)x = b where A = -H.

MatrixFreeEffectiveHamiltonian{D, W, T}

Matrix-free effective Hamiltonian operator H = RHS⁻¹ * LHS.

This is the fully matrix-free version of EffectiveHamiltonian, using O(N) memory instead of O(N²) for dense matrices.

For photonic crystals:

• TE: RHS = μ, so RHS⁻¹ = 1/μ (element-wise in real space)
• TM: RHS = ε, so RHS⁻¹ = 1/ε (element-wise in real space)

For phononic crystals:

• SH/PSV: RHS = ρ, so RHS⁻¹ = 1/ρ (element-wise in real space)

Grid Resolution Constraint

Important: The solver's grid resolution must be large enough to contain all plane wave indices. The required resolution in each dimension is:

resolution >= 2 * max_index + 1

where max_index is the maximum absolute value of plane wave indices (typically equal to the cutoff parameter). For example, with cutoff=7 (indices -7 to 7), the resolution must be at least 15 in each dimension.

If the resolution is too small, Fourier coefficients outside the grid will be silently dropped, leading to incorrect results (errors can be ~100%).

Usage

op = MatrixFreeOperator(solver, k)
H = MatrixFreeEffectiveHamiltonian(op)
A = NegatedOperator(H)  # A = -H for RSCG
x, stats = rscg(A, b, shifts)

RHS⁻¹ methods

• :approximate: Element-wise 1/ε in real space (fast, approximate for inhomogeneous media)
• :cg: Iterative CG to solve RHS * y = x (slower, exact)

NegatedOperator{T}

Wrapper that negates an operator: A = -H.

This transforms the Green's function problem: G(z) = (z·I - H)⁻¹ → solve (z·I + A)x = b where A = -H

Krylov.jl's cg_lanczos_shift solves (A + σI)x = b, so we use: A = -H, σ = z → (z·I - H)x = b

Usage

H = EffectiveHamiltonian(LHS, RHS)
A = NegatedOperator(H)
# Now use A with Krylov.cg_lanczos_shift

GFMethod

Abstract type for Green's function computation methods. Used by compute_greens_function, compute_dos, and compute_ldos.

DirectGF()

Dense direct solve method (LU factorization). Most accurate but O(N²) memory and O(N³) per frequency.

RSKGF(; atol=1e-10, rtol=1e-10, itmax=0, verbose=0)

Dense RSK method using ReducedShiftedKrylov.jl. O(N²) memory, efficient for many frequencies.

MatrixFreeGF(; rhs_inv_method=:approximate, atol=1e-10, rtol=1e-10, itmax=0, verbose=0)

Matrix-free RSCG method. O(N) memory, O(N log N) per iteration.

Arguments

• rhs_inv_method: Method for RHS⁻¹ application
  • :approximate: Element-wise 1/ε (fast, approximate for inhomogeneous media)
  • :cg: Iterative CG (slower, exact)

Recommended Usage

• Peak/resonance detection: accurate even without full RSCG convergence
• Large-scale calculations where Dense methods are memory-prohibitive
• For accurate absolute values, use DirectGF() instead

RHSInvMethod

Abstract type for RHS⁻¹ application methods in matrix-free effective Hamiltonian.

Concrete subtypes:

• ApproximateRHSInv: Fast element-wise 1/ε in real space (approximate)
• CGRHSInv: Iterative CG to solve RHS * y = x (exact)

See also: MatrixFreeEffectiveHamiltonian, MatrixFreeGF

ApproximateRHSInv <: RHSInvMethod

Approximate method for RHS⁻¹ application: element-wise 1/ε (or 1/ρ) in real space.

Fast but inaccurate for high-contrast media (e.g., Si/Air with ε ratio ≈ 12:1).

Example

method = MatrixFreeGF(rhs_inv_method=ApproximateRHSInv())

CGRHSInv <: RHSInvMethod

CG-based method for RHS⁻¹ application: solve RHS * y = x iteratively.

Accurate but slower than approximate method. Recommended for high-contrast media.

Fields

• atol::Float64: Absolute tolerance (default: 1e-10)
• rtol::Float64: Relative tolerance (default: 1e-10)
• maxiter::Int: Maximum iterations (default: 100)

Example

# Default parameters
method = MatrixFreeGF(rhs_inv_method=CGRHSInv())

# Custom parameters for high-contrast media
method = MatrixFreeGF(rhs_inv_method=CGRHSInv(atol=1e-12, rtol=1e-12, maxiter=200))

compute_greens_function(solver::Solver, k, ω_values, source; η=1e-3)

Compute Green's function G(ω) = (ω² + iη - H)⁻¹ * source for multiple frequencies.

Uses direct solve (dense method). For large systems, use matrix-free iterative methods.

Arguments

• solver: The PhoXonic solver
• k: Wave vector
• ω_values: Frequencies at which to compute G
• source: Source vector in plane wave basis
• η: Broadening parameter (small positive imaginary part)

Returns

• Vector of Green's function values G(ω) * source for each ω

compute_greens_function(solver, k, ω_values, source, method::GFMethod; η=1e-3)

Compute Green's function G(ω) * source using the specified method.

Methods

• DirectGF(): Dense direct solve (most accurate, O(N²) memory)
• RSKGF(): Dense RSK using ReducedShiftedKrylov.jl
• MatrixFreeGF(): Matrix-free RSCG (O(N) memory, best for large systems)

Examples

# Direct solve (default)
G = compute_greens_function(solver, k, ω_values, source; η=1e-2)

# Matrix-free for large systems
G = compute_greens_function(solver, k, ω_values, source, MatrixFreeGF(); η=1e-2)

compute_dos(solver::Solver, ω_values, k_points; η=1e-3)

Compute density of states (DOS) by summing over k-points.

DOS(ω) = -1/π Im[Tr G(ω)] summed over k-points

Arguments

• solver: The PhoXonic solver
• ω_values: Frequencies at which to compute DOS
• k_points: List of k-points to sum over (should cover Brillouin zone)
• η: Broadening parameter

Returns

• Vector of DOS values at each frequency

compute_dos(solver, ω_values, k_points, method::GFMethod; η=1e-3, n_random=10)

Compute density of states (DOS) using the specified method.

Methods

• DirectGF(): Dense direct solve (most accurate)
• RSKGF(): Dense RSK with stochastic trace estimation
• MatrixFreeGF(): Matrix-free RSCG with stochastic trace estimation

Examples

# Direct solve (default)
dos = compute_dos(solver, ω_values, k_points; η=1e-2)

# Matrix-free for large systems
dos = compute_dos(solver, ω_values, k_points, MatrixFreeGF(); η=1e-2)

compute_ldos(solver::Solver, position, ω_values, k_points; η=1e-3)

Compute local density of states (LDOS) at a given position.

LDOS(r, ω) = -1/π Im[G(r, r, ω)]

Arguments

• solver: The PhoXonic solver
• position: Position in real space (Vec2 or [x, y])
• ω_values: Frequencies at which to compute LDOS
• k_points: List of k-points to sum over
• η: Broadening parameter

Returns

• Vector of LDOS values at each frequency

compute_ldos(solver, position, ω_values, k_points, method::GFMethod; η=1e-3)

Compute local density of states (LDOS) using the specified method.

Methods

• DirectGF(): Dense direct solve (most accurate, O(N²) memory)
• RSKGF(): Dense RSK using ReducedShiftedKrylov.jl
• MatrixFreeGF(): Matrix-free RSCG (O(N) memory, best for large systems)

Examples

# Direct solve (default when method not specified)
ldos = compute_ldos(solver, pos, ω_values, k_points; η=1e-2)

# Matrix-free for large systems
ldos = compute_ldos(solver, pos, ω_values, k_points, MatrixFreeGF(); η=1e-2)

# Matrix-free with exact RHS⁻¹
ldos = compute_ldos(solver, pos, ω_values, k_points, MatrixFreeGF(rhs_inv_method=:cg); η=1e-2)

# Dense RSK
ldos = compute_ldos(solver, pos, ω_values, k_points, RSKGF(); η=1e-2)

See Also

• compute_dos: Compute density of states (trace of Green's function)
• compute_greens_function: Compute Green's function directly

compute_dos_stochastic(solver::Solver, ω_values, k_points; η=1e-3, n_random=10)

Compute DOS using stochastic trace estimation.

This approximates Tr[G] by averaging over random vectors: Tr[G] ≈ (1/nrandom) Σi <vi | G | vi>

More efficient than exact trace for large systems.

Arguments

• solver: The PhoXonic solver
• ω_values: Frequencies at which to compute DOS
• k_points: List of k-points
• η: Broadening parameter
• n_random: Number of random vectors for stochastic averaging

compute_zak_phase(solver::Solver{Dim1}, bands; n_k=50)

Compute the Zak phase for a 1D system.

The Zak phase is the Berry phase acquired by a Bloch state when k traverses the full Brillouin zone from 0 to 2π/a. For systems with inversion symmetry, the Zak phase is quantized to 0 or π.

Arguments

• solver: 1D solver (Photonic1D or Longitudinal1D)
• bands: Which bands to compute (e.g., 1:3)
• n_k: Number of k-points for integration (default: 50)

Returns

• ZakPhaseResult: Contains phases for each band

Example

lat = lattice_1d(1.0)
geo = Geometry(lat, Dielectric(1.0), [(Segment(0.2, 0.8), Dielectric(4.0))])
solver = Solver(Photonic1D(), geo, 64; cutoff=10)
result = compute_zak_phase(solver, 1:2)
println(result.phases)  # Zak phases for bands 1 and 2

ZakPhaseResult

Result of 1D Zak phase calculation.

Fields

• phases::Vector{Float64}: Zak phase for each band, in [-π, π]
• bands::UnitRange{Int}: Which bands were computed

Note

Per-band Zak phases have an inherent ±π gauge ambiguity due to endpoint gauge dependence in the truncated PWE basis. The multi-band total sum(phases) is gauge-invariant and reliable.

compute_wilson_spectrum(solver::Solver{Dim2}, bands; n_k_path=21, n_k_loop=50, loop_direction=:b2)

Compute the Wilson loop spectrum for a 2D system.

The Wilson spectrum shows the eigenvalues of the Wilson loop operator as a function of momentum along a path in the Brillouin zone. The Wilson loop is computed along the perpendicular direction at each k-point.

Arguments

• solver: 2D solver (TEWave, TMWave, SHWave, or PSVWave)
• bands: Which bands to include in the Wilson loop (e.g., 1:2)
• n_k_path: Number of k-points along the scanning path (default: 21)
• n_k_loop: Number of k-points for Wilson loop integration (default: 50)
• loop_direction: Direction of Wilson loop (:b1 or :b2, default: :b2)

Returns

• WilsonSpectrumResult: Contains kvalues, phases matrix, bands, and loopdirection

Example

lat = square_lattice(1.0)
geo = Geometry(lat, Dielectric(1.0), [(Circle([0.0, 0.0], 0.3), Dielectric(9.0))])
solver = Solver(TMWave(), geo, (32, 32); cutoff=5)
result = compute_wilson_spectrum(solver, 1:2; n_k_path=21, n_k_loop=50)

WilsonSpectrumResult

Result of 2D Wilson loop spectrum calculation.

Fields

• k_values::Vector{Float64}: k-values along the scanning direction
• phases::Matrix{Float64}: Wilson phases (nkpath × n_bands)
• bands::UnitRange{Int}: Which bands were computed
• loop_direction::Symbol: Direction of Wilson loop integration (:b1 or :b2)

winding_number(result::WilsonSpectrumResult, band_index::Int)

Calculate the winding number of a Wilson phase band.

The winding number counts how many times the Wilson phase winds around the Brillouin zone as k varies along the scanning path. A non-zero winding number indicates non-trivial topology.

Arguments

• result: Wilson spectrum result from compute_wilson_spectrum
• band_index: Which band (1 to n_bands) to compute winding for

Returns

• Int: The winding number (positive for upward winding, negative for downward)

Example

result = compute_wilson_spectrum(solver, 1:2)
w = winding_number(result, 1)  # Winding of first Wilson band