diff --git a/.gitignore b/.gitignore index 9b1a5785..01bcf79f 100644 --- a/.gitignore +++ b/.gitignore @@ -10,3 +10,6 @@ _git2_* docs/build *.info .vscode/spellright.dict +.DS_Store +.gitignore +deps/build.log diff --git a/docs/src/model/nonlinear_responses.md b/docs/src/model/nonlinear_responses.md index 128ae3c0..558188d7 100644 --- a/docs/src/model/nonlinear_responses.md +++ b/docs/src/model/nonlinear_responses.md @@ -1,7 +1,73 @@ # Nonlinear responses ## Kerr effect +The instantaneous electronic Kerr effect arises from the real part of the third-order susceptibility $\chi^{(3)}$. + +In the envelope picture, the Kerr nonlinear polarisation is: + +$$P_\mathrm{Kerr}(t) = \frac{3}{4}\varepsilon_0 \gamma_3 \rho \lvert E \rvert^2 E$$ + +where $\gamma_3$ is the single-molecule third-order hyperpolarisability (as returned by `PhysData.γ3_gas`), $\rho$ is the number density (so that $\rho \gamma_3 = \chi^{(3)}$), and $E$ is the field envelope. This response is frequency-independent and operates in the time domain. ## Photoionisation & plasma -## Raman response \ No newline at end of file + +## Raman response + + +## Two-photon absorption (TPA) + +Two-photon absorption is the **imaginary part** of $\chi^{(3)}$, the absorptive counterpart to the Kerr effect (which is the real part). While the Kerr response is frequency-independent and can be evaluated in the time domain, $\mathrm{Im}(\chi^{(3)})(\omega)$ varies strongly across deep-UV bandwidths — by a factor of ~20 across 200–280 nm in fused silica — making a time-domain scalar response inadequate. + +### Physical model + +TPA is characterised by the two-photon absorption coefficient $\beta_2(\omega)$ (units: m/W), which is energetically allowed when the two-photon energy exceeds the material bandgap: + +$$\beta_2(\omega) = C \times (2\hbar\omega/\mathrm{eV} - E_g)^\alpha \times \Theta(2\hbar\omega - E_g)$$ + +where $E_g$ is the bandgap energy (eV), $C$ (m/W/eV$^\alpha$) and $\alpha$ are empirical fitting parameters, and $\Theta$ is the Heaviside step function. + +### Implementation as a spectral response + +TPA is implemented as a [`SpectralResponse`](@ref) — a frequency-domain nonlinear response that operates on $\mathcal{F}\{|E|^2 E\}(\omega)$. This is in contrast to time-domain responses (Kerr, Raman, Plasma) which operate on $E(t)$. + +The TPA contribution to the nonlinear polarisation is: + +$$P_\mathrm{TPA}(\omega) = \mathrm{coeff}(\omega) \times \rho \times \mathcal{F}\{|E|^2 E\}(\omega)$$ + +where the coefficient is derived by tracing the standard relationship $\beta_2 = 3\omega \, \mathrm{Im}(\chi^{(3)}) / (2 n_0^2 c^2)$ through Luna's propagation pipeline. In `TransModeAvg`, the field is normalised by $\mathrm{nlscale} = \sqrt{\varepsilon_0 c/2}$ before computing $|E|^2 E$, and the norm function applies $-i\omega^2 / (4 \cdot \mathrm{nlscale} \cdot c \cdot \beta)$. The resulting $\mathrm{nlscale}^4 = (\varepsilon_0 c/2)^2$ factor in the normalisation chain produces the $\varepsilon_0^2 c^2$ prefactor: + +$$\mathrm{coeff}(\omega) = \frac{-i \, \varepsilon_0^2 \, c^2 \, \beta_2(\omega)}{2\omega}$$ + +This coefficient enters the **same P_NL buffer** as the Kerr response, so the `Trans*` types' normalisation functions correctly convert both Kerr and TPA contributions to $\partial A/\partial z$. The $-i$ factor ensures TPA causes loss (negative real $\partial A / \partial z$) after the normalisation's own $-i$ factor. + +### Geometry-agnostic design + +`TPAResponse` is **geometry-agnostic**: it operates on 1D frequency vectors with callable signature `(out_ω, F_E2E_ω, ρ)`, regardless of whether the simulation uses mode-averaged (`TransModeAvg`), modal (`TransModal`), or 3D free-space (`TransFree`) propagation. All spatial integration, effective area, and field normalisation is handled by the `Trans*` types in `NonlinearRHS`, exactly as for time-domain responses. + +### Degenerate approximation + +The current implementation uses the **degenerate** TPA approximation, where all interacting photons have the same frequency. This still accounts for the overall frequency dependence of the TPA response (e.g. two photons at 200 nm, or two photons at 300 nm), but does not account for the frequency difference between photons in non-degenerate TPA interactions (e.g. one photon at 200 nm and another at 300 nm). This is appropriate for narrowband or moderate-bandwidth pulses. Non-degenerate TPA (photons at different frequencies) is not currently implemented. + +### Notes +- Only valid for **envelope** propagation (`EnvGrid`). Using with `RealGrid` (carrier-resolved) is not physically meaningful. +- Handles **scalar** (single-polarisation) fields. Vector TPA ($(|E_x|^2 + |E_y|^2)E$) is a future extension. +- The coefficient is derived for Luna's `TransModeAvg` pipeline. It works correctly across all `Trans*` types because each type's normalisation function compensates for its own field convention — the same mechanism that makes the Kerr response geometry-agnostic. + +### Usage + +```julia +using Luna + +# From material parameters (PhysData) +β₂_ω = PhysData.β₂_TPA.(grid.ω, :SiO2) +tpa = Nonlinear.TPAResponse(grid.ω, β₂_ω) + +# Include in responses tuple alongside Kerr +responses = (Nonlinear.Kerr_env(χ3), tpa) + +# Or via prop_gnlse interface +prop_gnlse(γ, flength, βs; ..., tpa=β₂_ω) +``` + +See also: [`PhysData.β₂_TPA`](@ref), [`Nonlinear.TPAResponse`](@ref), [`Nonlinear.SpectralResponse`](@ref) diff --git a/examples/low_level_interface/freespace/radial_env_tpa.jl b/examples/low_level_interface/freespace/radial_env_tpa.jl new file mode 100644 index 00000000..88627f09 --- /dev/null +++ b/examples/low_level_interface/freespace/radial_env_tpa.jl @@ -0,0 +1,84 @@ +using Luna +import Luna.PhysData: wlfreq +import FFTW +import Hankel + +# Deep-UV Gaussian beam in SiO₂ — envelope, radially symmetric +λ0 = 240e-9 +τfwhm = 2e-15 # broad bandwidth across TPA region +energy = 150e-9 +w0 = 15e-6 # beam waist +L = 10e-6 # 100 µm propagation (thin slab) +R = 250e-6 # Hankel aperture +N = 256 # radial points + +material = :MgF2 + +grid = Grid.EnvGrid(L, λ0, (160e-9, 600e-9), 1e-12) +q = Hankel.QDHT(R, N, dim=2) + +densityfun(z) = 1.0 # solid — χ³ already bulk value +nfun = PhysData.ref_index_fun(material, lookup=false) +normfun = NonlinearRHS.const_norm_radial(grid, q, nfun) +linop = LinearOps.make_const_linop(grid, q, nfun) + +χ3 = PhysData.χ3(material) +β₂_ω = PhysData.β₂_TPA.(grid.ω, material) +tpa = Nonlinear.TPAResponse(grid.ω, β₂_ω) + +inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, w0=w0) + +# With TPA +responses_tpa = (Nonlinear.Kerr_env(χ3), tpa) +Eω_tpa, transform_tpa, FT_tpa = Luna.setup(grid, q, densityfun, normfun, +responses_tpa, inputs) +output_tpa = Output.MemoryOutput(0, grid.zmax, 201) +Luna.run(Eω_tpa, grid, linop, transform_tpa, FT_tpa, output_tpa) + +# Without TPA +responses_notpa = (Nonlinear.Kerr_env(χ3),) +Eω_notpa, transform_notpa, FT_notpa = Luna.setup(grid, q, densityfun, normfun, +responses_notpa, inputs) +output_notpa = Output.MemoryOutput(0, grid.zmax, 201) +Luna.run(Eω_notpa, grid, linop, transform_notpa, FT_notpa, output_notpa) + +# Energy comparison (sum over all radial k-modes) +Eωk_tpa_in = output_tpa.data["Eω"][:, :, 1] +Eωk_tpa_out = output_tpa.data["Eω"][:, :, end] +Eωk_notpa_out = output_notpa.data["Eω"][:, :, end] + +energy_in = sum(abs2, Eωk_tpa_in) +energy_out_tpa = sum(abs2, Eωk_tpa_out) +energy_out_notpa = sum(abs2, Eωk_notpa_out) + +# With TPA, output energy should be less than without +@info "Energy comparison: TPA output = $energy_out_tpa, No TPA output = $energy_out_notpa" + +# Spectral suppression: TPA stronger at shorter wavelengths (higher ω) +# Sum over radial modes for spectral comparison +spec_tpa = dropdims(sum(abs2, Eωk_tpa_out, dims=2), dims=2) +spec_notpa = dropdims(sum(abs2, Eωk_notpa_out, dims=2), dims=2) + +ω_below = 2π * PhysData.c / 300e-9 # below SiO₂ TPA edge (~299 nm) +ω_above = 2π * PhysData.c / 220e-9 # well above TPA edge + +idx_below = argmin(abs.(grid.ω .- ω_below)) +idx_above = argmin(abs.(grid.ω .- ω_above)) + +ratio_below = spec_tpa[idx_below] / spec_notpa[idx_below] +ratio_above = spec_tpa[idx_above] / spec_notpa[idx_above] +# Above-edge suppression should be stronger (lower ratio) +@info "Spectral suppression: ratio above edge = $ratio_above, ratio below edge = $ratio_below" + +λ, Iλtpa = Processing.getIω(output_tpa, :λ, L) +λ, Iλnotpa = Processing.getIω(output_notpa, :λ, L) + +using PyPlot +figure() +plot(λ*1e9, Iλtpa[:,1,1], label="With TPA") +plot(λ*1e9, Iλnotpa[:,1,1], label="Without TPA") +xlabel("Wavelength (nm)") +ylabel("On-axis spectral Intensity (a.u.)") +legend() +title("Spectral suppression from TPA in SiO₂") +xlim(180, 400) diff --git a/examples/simple_interface/gnlse_tpa.jl b/examples/simple_interface/gnlse_tpa.jl new file mode 100644 index 00000000..18490d5d --- /dev/null +++ b/examples/simple_interface/gnlse_tpa.jl @@ -0,0 +1,36 @@ +# Quick example showing how to include TPA in a GNLSE simulation. +# Comparing with the same simulation without TPA to see the spectral suppression at shorter wavelengths. + +using Luna, PyPlot + +λ0 = 240e-9 +τfwhm = 5e-15 # short pulse = broad bandwidth +energy = 100e-9 +flength = 50e-6 +β2 = 1e-30 +βs = [0.0, 0.0, β2] +Aeff = pi*12e-6^2 +λlims = [160e-9, 400e-9] +trange = 1e-12 + +# Using material-based convention: n2, γ, and TPA auto-computed from :SiO2 +output_tpa = prop_gnlse(flength, βs; material=:SiO2, Aeff, λ0, τfwhm, energy, + pulseshape=:gauss, λlims, trange, + raman=false, shock=false, shotnoise=false, tpa=true) + +output_notpa = prop_gnlse(flength, βs; material=:SiO2, Aeff, λ0, τfwhm, energy, + pulseshape=:gauss, λlims, trange, + raman=false, shock=false, shotnoise=false) + +λ, Iλtpa = Processing.getIω(output_tpa, :λ, flength) +λ, Iλnotpa = Processing.getIω(output_notpa, :λ, flength) + +using PyPlot +figure() +plot(λ*1e9, Iλtpa, label="With TPA") +plot(λ*1e9, Iλnotpa, label="Without TPA") +xlabel("Wavelength (nm)") +ylabel("Spectral Intensity (a.u.)") +legend() +title("Spectral suppression from TPA in SiO₂") +xlim(180, 400) diff --git a/src/Interface.jl b/src/Interface.jl index 0518040d..dfbaf102 100644 --- a/src/Interface.jl +++ b/src/Interface.jl @@ -836,11 +836,17 @@ end """ prop_gnlse(γ, flength, βs; λ0, λlims, trange, kwargs...) + prop_gnlse(flength, βs; n2, Aeff, λ0, λlims, trange, kwargs...) + prop_gnlse(flength, βs; material, Aeff, λ0, λlims, trange, kwargs...) Simulate pulse propagation using the GNLSE. +The first form takes the nonlinear coefficient `γ` directly. The second and third forms +compute `γ` internally from `n2` (or a `material` lookup) and `Aeff`, which is required +for two-photon absorption (TPA) support. + # Mandatory arguments -- `γ::Number`: The nonlinear coefficient. +- `γ::Number`: The nonlinear coefficient (W⁻¹m⁻¹). Only for the first calling convention. - `flength::Number`: Length of the fibre. - `βs`: The Taylor expansion of the propagation constant about `λ0`. - `λ0`: (keyword argument) the reference wavelength for the simulation. For simple @@ -849,6 +855,14 @@ Simulate pulse propagation using the GNLSE. - `trange::Number`: The total width of the time grid. To make the number of samples a power of 2, the actual grid used is usually bigger. +# Nonlinear parameter specification (alternative to `γ`) +When `γ` is not given as a positional argument, the nonlinear coefficient is computed from +`n2` and `Aeff`. Exactly one of `n2` or `material` must be provided. +- `n2::Real`: Nonlinear refractive index (m²/W). +- `material::Symbol`: Material name (e.g. `:SiO2`), from which `n2` is looked up via + [`PhysData.n2`](@ref). +- `Aeff::Real`: Effective mode area (m²). Required when using `n2` or `material`. + # Grid options - `δt::Number`: Time step on the fine grid used for the nonlinear interaction. By default, this is determined by the wavelength grid. If `δt` is given **and smaller** than the @@ -889,6 +903,18 @@ Note that the current GNLSE model is single mode only. - `τ1`: the Raman oscillator period. - `τ2`: the Raman damping time. +# Two-photon absorption (TPA) +TPA is only available when using the `n2`/`material` + `Aeff` calling convention (not when +passing `γ` directly, since `Aeff` is needed to scale β₂ to mode-averaged units). +- `tpa`: TPA specification. Can be: + - `nothing` or `false`: no TPA (default). + - `true`: auto-compute β₂(ω) from `material` via [`PhysData.β₂_TPA`](@ref). + Requires `material` to be specified. + - `AbstractVector`: pre-computed β₂(ω) values (m/W) sampled on the simulation + frequency grid `grid.ω`. + - callable `tpa(ω) → β₂`: a function returning the TPA coefficient β₂ (m/W) + at angular frequency ω. + # Output options - `saveN::Integer`: Number of points along z at which to save the field. - `filepath`: If `nothing` (default), create a `MemoryOutput` to store the simulation results @@ -912,11 +938,16 @@ end """ prop_gnlse_args(γ, flength, βs; λ0, λlims, trange, kwargs...) -Prepare to simulate pulse propagation using the GNLSE. This -function takes the same arguments as `prop_gnlse` but instead or running the -simulation and returning the output, it returns the required arguments for `Luna.run`, -which is useful for repeated simulations in an indentical fibre with different initial -conditions. +Prepare to simulate pulse propagation using the GNLSE with a directly specified +nonlinear coefficient `γ` (W⁻¹m⁻¹). This function takes the same arguments as +`prop_gnlse` but instead of running the simulation and returning the output, it +returns the required arguments for `Luna.run`. + +!!! note + TPA (`tpa` kwarg) is not supported in this calling convention because `γ` does + not carry enough information to determine `Aeff`. Use + `prop_gnlse(flength, βs; n2, Aeff, tpa, ...)` or + `prop_gnlse(flength, βs; material, Aeff, tpa, ...)` instead. """ function prop_gnlse_args(γ, flength, βs; λ0, λlims, trange, δt=1, τfwhm=nothing, τw=nothing, ϕ=Float64[], @@ -926,6 +957,7 @@ function prop_gnlse_args(γ, flength, βs; λ0, λlims, trange, shotnoise=true, shock=true, loss=0.0, raman=true, fr=0.18, ramanmodel=:sdo, τ1=12.2e-15, τ2=32e-15, + tpa=nothing, saveN=201, filepath=nothing, scan=nothing, scanidx=nothing, filename=nothing) envelope = true @@ -960,6 +992,125 @@ function prop_gnlse_args(γ, flength, βs; λ0, λlims, trange, error("unrecognised value for ramanmodel") end end + if !isnothing(tpa) + error("TPA requires Aeff for correct scaling. Use prop_gnlse(flength, βs; n2, Aeff, tpa, ...) or prop_gnlse(flength, βs; material, Aeff, tpa, ...) instead.") + end + resp = Tuple(resp) + + inputs = makeinputs(mode_s, λ0, pulses, τfwhm, τw, ϕ, + power, energy, pulseshape, polarisation, propagator) + inputs = shotnoise_maybe(inputs, mode_s, shotnoise) + + norm! = NonlinearRHS.norm_mode_average_gnlse(grid, aeff; shock) + Eω, transform, FT = Luna.setup(grid, density, resp, inputs, βfun!, aeff, norm! = norm!) + stats = Stats.default(grid, Eω, mode_s, linop, transform) + output = makeoutput(grid, saveN, stats, filepath, scan, scanidx, filename) + + saveargs(output; γ, flength, βs, λlims, trange, envelope, thg, δt, + λ0, τfwhm, τw, ϕ, power, energy, pulseshape, polarisation, propagator, pulses, + shotnoise, shock, loss, raman, ramanmodel, fr, τ1, τ2, tpa, saveN, filepath, filename) + + return Eω, grid, linop, transform, FT, output +end + +""" + prop_gnlse_args(flength, βs; n2, Aeff, λ0, λlims, trange, kwargs...) + prop_gnlse_args(flength, βs; material, Aeff, λ0, λlims, trange, kwargs...) + +Prepare to simulate pulse propagation using the GNLSE with physically consistent +nonlinear parameters. The nonlinear coefficient `γ` is computed internally from +`n2` and `Aeff`, ensuring correct scaling for all nonlinear effects including TPA. + +# Nonlinear parameter specification (exactly one required) +- `n2::Real` : nonlinear refractive index (m²/W) +- `material::Symbol` : material name (e.g. `:SiO2`), from which `n2` is looked up + via [`PhysData.n2`](@ref) + +# Required keyword +- `Aeff::Real` : effective mode area (m²) + +# TPA keyword +- `tpa` : two-photon absorption specification: + - `nothing` or `false` : no TPA (default) + - `true` : auto-compute from `material` via [`PhysData.β₂_TPA`](@ref) + - `AbstractVector` : pre-computed β₂(ω) values (m/W) on `grid.ω` + - callable `tpa(ω) → β₂` : function returning β₂ (m/W) + + In all cases, β₂ is scaled by `1/Aeff` internally to convert from intensity-based + to mode-averaged units. + +All other keyword arguments are the same as for `prop_gnlse_args(γ, flength, βs; ...)`. +""" +function prop_gnlse_args(flength, βs; λ0, λlims, trange, + n2=nothing, Aeff=nothing, material::Union{Symbol,Nothing}=nothing, + δt=1, τfwhm=nothing, τw=nothing, ϕ=Float64[], + power=nothing, energy=nothing, + pulseshape=:gauss, propagator=nothing, + pulses=nothing, + shotnoise=true, shock=true, + loss=0.0, raman=true, fr=0.18, + ramanmodel=:sdo, τ1=12.2e-15, τ2=32e-15, + tpa=nothing, + saveN=201, filepath=nothing, + scan=nothing, scanidx=nothing, filename=nothing) + # Validate and resolve n2 + if !isnothing(material) + !isnothing(n2) && error("cannot specify both n2 and material") + n2 = PhysData.n2(material) + elseif isnothing(n2) + error("must specify either n2 or material") + end + isnothing(Aeff) && error("Aeff must be provided") + + envelope = true + thg = false + polarisation = :linear + grid = makegrid(flength, λ0, λlims, trange, envelope, thg, δt) + mode_s = SimpleFibre.SimpleMode(PhysData.wlfreq(λ0), βs; loss) + aeff = z -> 1.0 + density = z -> 1.0 + linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, mode_s, λ0) + + # Compute γ from n2 and Aeff + k0 = 2π/λ0 + γ = n2 * k0 / Aeff + + # Kerr: same derivation as γ-based method + # factor of 4/3 compensates for the factor of 3/4 in Nonlinear.jl, as + # n2 and γ are usually defined for the envelope case already + χ3 = 4/3 * (1 - fr) * n2 / Aeff * (PhysData.ε_0*PhysData.c) + resp = Any[Nonlinear.Kerr_env(χ3)] + if raman + # factor of 2 here compensates for factor 1/2 in Nonlinear.jl as fr is + # defined for the envelope case already + χ3R = 2 * fr * n2 / Aeff * (PhysData.ε_0*PhysData.c) + if ramanmodel == :SiO2 + push!(resp, Nonlinear.RamanPolarEnv(grid.to, Raman.raman_response(grid.to, :SiO2, + χ3R * PhysData.ε_0))) + elseif ramanmodel == :sdo + if isnothing(τ1) || isnothing(τ2) + error("for :sdo ramanmodel you must specify τ1 and τ2") + end + push!(resp, Nonlinear.RamanPolarEnv(grid.to, + Raman.CombinedRamanResponse(grid.to, + [Raman.RamanRespNormedSingleDampedOscillator(χ3R * PhysData.ε_0, 1/τ1, τ2)]))) + else + error("unrecognised value for ramanmodel") + end + end + + # Two-photon absorption: scale β₂ by 1/Aeff for mode-averaged propagation + if !isnothing(tpa) && tpa !== false + if tpa === true + isnothing(material) && error("tpa=true requires material to be specified") + β₂_ω = PhysData.β₂_TPA.(grid.ω, material) + elseif tpa isa AbstractVector + β₂_ω = tpa + else + β₂_ω = tpa.(grid.ω) + end + push!(resp, Nonlinear.TPAResponse(grid.ω, β₂_ω ./ Aeff)) + end resp = Tuple(resp) inputs = makeinputs(mode_s, λ0, pulses, τfwhm, τw, ϕ, @@ -973,7 +1124,8 @@ function prop_gnlse_args(γ, flength, βs; λ0, λlims, trange, saveargs(output; γ, flength, βs, λlims, trange, envelope, thg, δt, λ0, τfwhm, τw, ϕ, power, energy, pulseshape, polarisation, propagator, pulses, - shotnoise, shock, loss, raman, ramanmodel, fr, τ1, τ2, saveN, filepath, filename) + shotnoise, shock, loss, raman, ramanmodel, fr, τ1, τ2, tpa, n2, Aeff, material, + saveN, filepath, filename) return Eω, grid, linop, transform, FT, output end diff --git a/src/Luna.jl b/src/Luna.jl index 24de3188..76a5b0bb 100644 --- a/src/Luna.jl +++ b/src/Luna.jl @@ -357,7 +357,7 @@ save_modeinfo_maybe(output, t) = nothing function run(Eω, grid, linop, transform, FT, output; - min_dz=0, max_dz=grid.zmax/2, init_dz=1e-4, z0=0.0, + min_dz=0, max_dz=grid.zmax/2, init_dz=grid.zmax/4, z0=0.0, rtol=1e-6, atol=1e-10, safety=0.9, norm=RK45.weaknorm, status_period=1) diff --git a/src/Nonlinear.jl b/src/Nonlinear.jl index a90e7198..359e1532 100644 --- a/src/Nonlinear.jl +++ b/src/Nonlinear.jl @@ -1,6 +1,6 @@ module Nonlinear import Luna -import Luna.PhysData: ε_0, e_ratio +import Luna.PhysData: ε_0, c, e_ratio import Luna: Maths, Utils import FFTW import LinearAlgebra: mul!, ldiv! @@ -332,4 +332,143 @@ function (R::RamanPolar)(out, Et, ρ) end end +""" + SpectralResponse + +Abstract type for nonlinear responses that operate in the frequency domain +rather than the time domain. + +Unlike time-domain responses (Kerr, Raman, Plasma) with callable signature +`response!(out_t, E_t, ρ)`, spectral responses have signature: + + response!(out_ω, F_E2E_ω, ρ) + +where `F_E2E_ω` is ``\\mathcal{F}\\{|E|^2 E\\}(\\omega)`` and `out_ω` is the +frequency-domain nonlinear polarisation (same units as FFT of time-domain P_NL). + +Spectral responses are **geometry-agnostic**: they operate on 1D frequency +vectors regardless of whether the simulation is mode-averaged, modal, or +3D free-space. All spatial integration and field normalisation is handled +by the `Trans*` types in `NonlinearRHS`, exactly as for time-domain responses. +The spectral response contribution enters `P_NL(ω)` **before** the +normalisation step, so it goes through the same norm as Kerr. + +See also: [`TPAResponse`](@ref) +""" +abstract type SpectralResponse end + +""" + TPAResponse{cT} + +Two-photon absorption (TPA) response for envelope propagation. + +TPA is the imaginary part of the third-order susceptibility χ⁽³⁾. While the +real part (Kerr effect) is frequency-independent and handled in the time domain, +Im(χ⁽³⁾)(ω) varies strongly across deep-UV bandwidths and must be applied as +a frequency-domain multiplier on ``\\mathcal{F}\\{|E|^2 E\\}(\\omega)``: + + P_TPA(ω) = coeff(ω) × ρ × F{|E|²E}(ω) + +This produces output in the **same P_NL units** as the Kerr response, so the +`Trans*` types' normalisation functions correctly convert both Kerr and TPA +contributions to ∂A/∂z. This makes `TPAResponse` geometry-agnostic: it works +identically for mode-averaged (`TransModeAvg`), modal (`TransModal`), and +3D free-space (`TransFree`) propagation. + +The stored coefficient is: + + coeff(ω) = -i × ε₀² × c² × β₂(ω) / (2ω) + +derived by requiring consistency with the Kerr nonlinear polarisation +in Luna's propagation pipeline. In ``TransModeAvg``, the field is +divided by ``\\mathrm{nlscale} = \\sqrt{\\varepsilon_0 c/2}`` before +computing ``|E|^2 E``, and the norm function applies +``-i\\omega^2/(4 \\cdot \\mathrm{nlscale} \\cdot c \\cdot \\beta)``. +Tracing the standard relationship +``\\beta_2 = 3\\omega \\, \\mathrm{Im}(\\chi^{(3)}) / (2 n_0^2 c^2)`` +through yields the ``\\varepsilon_0^2 c^2`` prefactor +(from ``\\mathrm{nlscale}^4 = (\\varepsilon_0 c/2)^2``). +The ``-i`` factor ensures TPA causes loss (negative real +``\\partial A/\\partial z``) after the norm's own ``-i`` factor. + +# Construction +```julia +# From β₂(ω) values and frequency grid: +tpa = Nonlinear.TPAResponse(grid.ω, PhysData.β₂_TPA.(grid.ω, :SiO2)) + +# From a pre-computed coefficient array: +tpa = Nonlinear.TPAResponse(coeff_ω) +``` + +# Usage +Pass in the responses tuple alongside time-domain responses: +```julia +responses = (Nonlinear.Kerr_env(χ3), Nonlinear.TPAResponse(grid.ω, β₂_ω)) +``` +All `Trans*` types automatically route `SpectralResponse` subtypes to +the frequency-domain path. + +# Notes +- Only valid for envelope propagation (`EnvGrid`). Using with `RealGrid` + (carrier-resolved) is not physically meaningful. +- Handles scalar (single-polarisation) fields. Vector TPA + (``(|E_x|^2 + |E_y|^2) E``) is a future extension. +- The ``\\varepsilon_0^2 c^2`` prefactor arises from the ``\\mathrm{nlscale}`` + field normalisation in the ``Trans*`` types — it is not related to a + ``\\gamma_3 = \\chi^{(3)}/n_0^2`` convention. In Luna, ``\\gamma_3`` is + the single-molecule third-order hyperpolarisability, and + ``\\chi^{(3)} = \\rho \\, \\gamma_3``. + +See also: [`PhysData.β₂_TPA`](@ref), [`SpectralResponse`](@ref) +""" +struct TPAResponse{cT} <: SpectralResponse + coeff_ω::cT +end + +""" + TPAResponse(ω_grid, β₂_ω) → TPAResponse + +Construct a `TPAResponse` from an angular frequency grid `ω_grid` (rad/s) +and a vector of β₂ values `β₂_ω` (m/W). + +Converts using: + + coeff(ω) = -i × ε₀² × c² × β₂(ω) / (2ω) + +Sets coeff = 0 at ω = 0 to avoid the singularity (β₂ = 0 there anyway, +since no material has a TPA edge at DC). +""" +function TPAResponse(ω_grid::AbstractVector, β₂_ω::AbstractVector) + length(ω_grid) == length(β₂_ω) || throw( + DimensionMismatch("ω_grid and β₂_ω must have the same length")) + coeff = zeros(ComplexF64, length(ω_grid)) + for i in eachindex(ω_grid) + if ω_grid[i] != 0 + coeff[i] = -im * ε_0^2 * c^2 * β₂_ω[i] / (2 * ω_grid[i]) + end + end + return TPAResponse(coeff) +end + +""" + (tpa::TPAResponse)(out_ω, F_E2E_ω, ρ) + +Add TPA contribution to the frequency-domain nonlinear polarisation `out_ω`. + +- `out_ω`: frequency-domain P_NL buffer (1D), modified in-place +- `F_E2E_ω`: ``\\mathcal{F}\\{|E|^2 E\\}(\\omega)`` on the standard ω grid +- `ρ`: density (scalar), same as passed to time-domain responses + +Adds: `out_ω[i] += ρ × coeff[i] × F_E2E_ω[i]` + +This signature mirrors the time-domain `(out_t, E_t, ρ)` pattern but +operates in the frequency domain. The `Trans*` types handle all spatial +looping and normalisation. +""" +function (tpa::TPAResponse)(out_ω, F_E2E_ω, ρ) + @inbounds for i in eachindex(out_ω) + out_ω[i] += ρ * tpa.coeff_ω[i] * F_E2E_ω[i] + end +end + end diff --git a/src/NonlinearRHS.jl b/src/NonlinearRHS.jl index c020b333..74e03c88 100644 --- a/src/NonlinearRHS.jl +++ b/src/NonlinearRHS.jl @@ -7,6 +7,7 @@ import LinearAlgebra: mul!, ldiv! import NumericalIntegration: integrate, SimpsonEven import Luna: PhysData, Modes, Maths, Grid import Luna.PhysData: wlfreq +import Luna.Nonlinear: SpectralResponse """ to_time!(Ato, Aω, Aωo, IFTplan) @@ -140,12 +141,40 @@ function Et_to_Pt!(Pt, Et, responses, density, idcs) end end +""" + sp_Et_to_Pω!(Pω, F_E2E_ω, sp_resp, density) + +Apply all [`SpectralResponse`](@ref)s (e.g. TPA) to the frequency-domain +cubic envelope product ``\\mathcal{F}\\{|E|^2 E\\}(\\omega)``. Adds contributions +to the frequency-domain polarisation buffer `Pω` (same units as the Kerr +contribution), **before** normalisation. + +This is the spectral counterpart of [`Et_to_Pt!`](@ref). +""" +function sp_Et_to_Pω!(Pω, F_E2E_ω, sp_resp, density) + for resp in sp_resp + resp(Pω, F_E2E_ω, density) + end +end + +""" + sp_Et_to_Pω!(Pω, F_E2E_ω, sp_resp, density, idcs) + +Spatial-indexed version for 3D free-space (`TransFree`) and 2D column data. +Loops over spatial indices, applying spectral responses to each 1D frequency slice. +""" +function sp_Et_to_Pω!(Pω, F_E2E_ω, sp_resp, density, idcs) + for i in idcs + sp_Et_to_Pω!(view(Pω, :, i), view(F_E2E_ω, :, i), sp_resp, density) + end +end + """ TransModal Transform E(ω) -> Pₙₗ(ω) for modal fields. """ -mutable struct TransModal{tsT, lT, TT, FTT, rT, gT, dT, ddT, nT} +mutable struct TransModal{tsT, lT, TT, FTT, rT, gT, dT, ddT, nT, spT} ts::tsT full::Bool dimlimits::lT @@ -169,6 +198,7 @@ mutable struct TransModal{tsT, lT, TT, FTT, rT, gT, dT, ddT, nT} atol::Float64 mfcn::Int err::Array{ComplexF64,2} + sp_resp::spT # spectral responses (e.g. TPA), applied in ω domain before norm end function show(io::IO, t::TransModal) @@ -179,7 +209,12 @@ function show(io::IO, t::TransModal) samples = "time grid size: $(length(t.grid.t)) / $(length(t.grid.to))" resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n ") full = "full: $(t.full)" - out = join(["TransModal", modes, pol, grid, samples, full, resp], "\n ") + parts = ["TransModal", modes, pol, grid, samples, full, resp] + if !isempty(t.sp_resp) + sp = "spectral responses: "*join([string(typeof(ri)) for ri in t.sp_resp], "\n ") + push!(parts, sp) + end + out = join(parts, "\n ") print(io, out) end @@ -211,8 +246,12 @@ function TransModal(tT, grid, ts::Modes.ToSpace, FT, resp, densityfun, norm!; Prωo = Array{ComplexF64,2}(undef, length(grid.ωo), ts.npol) Prmω = Array{ComplexF64,2}(undef, length(grid.ω), ts.nmodes) IFT = inv(FT) + # Split responses into time-domain and spectral (frequency-domain) + td_resp = Tuple(r for r in resp if !(r isa SpectralResponse)) + sp_resp = Tuple(r for r in resp if r isa SpectralResponse) TransModal(ts, full, Modes.dimlimits(ts.ms[1]), Emω, Erω, Erωo, Er, Pr, Prω, Prωo, Prmω, - FT, resp, grid, densityfun, densityfun(0.0), norm!, 0, 0.0, rtol, atol, mfcn, similar(Prmω)) + FT, td_resp, grid, densityfun, densityfun(0.0), norm!, 0, 0.0, rtol, atol, mfcn, + similar(Prmω), sp_resp) end function TransModal(grid::Grid.RealGrid, args...; kwargs...) @@ -277,6 +316,17 @@ function Erω_to_Prω!(t, x) Et_to_Pt!(t.Pr, t.Er, t.resp, t.density) @. t.Pr *= t.grid.towin to_freq!(t.Prω, t.Prωo, t.Pr, t.FT) + # Spectral responses (e.g. TPA): add to Prω before ωwin and norm. + # Erω is stale (consumed by to_time! above), so reuse it for F{|E|²E}. + if !isempty(t.sp_resp) + @. t.Pr = abs2(t.Er) * t.Er + @. t.Pr *= t.grid.towin + to_freq!(t.Erω, t.Erωo, t.Pr, t.FT) + for col in axes(t.Prω, 2) + sp_Et_to_Pω!(view(t.Prω, :, col), view(t.Erω, :, col), + t.sp_resp, t.density) + end + end @. t.Prω *= t.grid.ωwin t.norm!(t.Prω) end @@ -320,7 +370,7 @@ end Transform E(ω) -> Pₙₗ(ω) for mode-averaged single-mode propagation. """ -struct TransModeAvg{TT, FTT, rT, gT, dT, nT, aT} +struct TransModeAvg{TT, FTT, rT, gT, dT, nT, aT, spT} Pto::Vector{TT} Eto::Vector{TT} Eωo::Vector{ComplexF64} @@ -331,13 +381,20 @@ struct TransModeAvg{TT, FTT, rT, gT, dT, nT, aT} densityfun::dT norm!::nT aeff::aT # function which returns effective area + sp_resp::spT # spectral responses (e.g. TPA), applied in ω domain before norm + A2A_ω::Vector{ComplexF64} # scratch buffer for F{|E|²E}(ω) end function show(io::IO, t::TransModeAvg) grid = "grid type: $(typeof(t.grid))" samples = "time grid size: $(length(t.grid.t)) / $(length(t.grid.to))" resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n ") - out = join(["TransModeAvg", grid, samples, resp], "\n ") + parts = ["TransModeAvg", grid, samples, resp] + if !isempty(t.sp_resp) + sp = "spectral responses: "*join([string(typeof(ri)) for ri in t.sp_resp], "\n ") + push!(parts, sp) + end + out = join(parts, "\n ") print(io, out) end @@ -346,7 +403,12 @@ function TransModeAvg(TT, grid, FT, resp, densityfun, norm!, aeff) Eto = zeros(TT, length(grid.to)) Pto = similar(Eto) Pωo = similar(Eωo) - TransModeAvg(Pto, Eto, Eωo, Pωo, FT, resp, grid, densityfun, norm!, aeff) + # Split responses into time-domain and spectral (frequency-domain) + td_resp = Tuple(r for r in resp if !(r isa SpectralResponse)) + sp_resp = Tuple(r for r in resp if r isa SpectralResponse) + A2A_ω = isempty(sp_resp) ? ComplexF64[] : zeros(ComplexF64, length(grid.ω)) + TransModeAvg(Pto, Eto, Eωo, Pωo, FT, td_resp, grid, densityfun, norm!, aeff, + sp_resp, A2A_ω) end function TransModeAvg(grid::Grid.RealGrid, FT, resp, densityfun, norm!, aeff) @@ -365,6 +427,14 @@ function (t::TransModeAvg)(nl, Eω, z) Et_to_Pt!(t.Pto, t.Eto, t.resp, t.densityfun(z)) @. t.Pto *= t.grid.towin to_freq!(nl, t.Pωo, t.Pto, t.FT) + # Spectral responses (e.g. TPA): add to nl before norm, so they go through + # the same normalisation pipeline as time-domain responses (Kerr, Raman). + if !isempty(t.sp_resp) + @. t.Pto = abs2(t.Eto) * t.Eto + @. t.Pto *= t.grid.towin + to_freq!(t.A2A_ω, t.Pωo, t.Pto, t.FT) + sp_Et_to_Pω!(nl, t.A2A_ω, t.sp_resp, t.densityfun(z)) + end t.norm!(nl, z) for i in eachindex(nl) !t.grid.sidx[i] && continue @@ -403,7 +473,7 @@ end Transform E(ω) -> Pₙₗ(ω) for radially symetric free-space propagation """ -struct TransRadial{TT, HTT, FTT, nT, rT, gT, dT, iT} +struct TransRadial{TT, HTT, FTT, nT, rT, gT, dT, iT, spT} QDHT::HTT # Hankel transform (space to k-space) FT::FTT # Fourier transform (time to frequency) normfun::nT # Function which returns normalisation factor @@ -415,6 +485,8 @@ struct TransRadial{TT, HTT, FTT, nT, rT, gT, dT, iT} Eωo::Array{ComplexF64,2} # Buffer array for field on oversampled frequency grid Pωo::Array{ComplexF64,2} # Buffer array for NL polarisation on oversampled frequency grid idcs::iT # CartesianIndices for Et_to_Pt! to iterate over + sp_resp::spT # spectral responses (e.g. TPA), applied in ω domain before norm + A2A_kω::Array{ComplexF64,2} # scratch buffer for F{|E|²E}(ω,kr) end function show(io::IO, t::TransRadial) @@ -423,7 +495,12 @@ function show(io::IO, t::TransRadial) resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n ") nr = "radial points: $(t.QDHT.N)" R = "aperture: $(t.QDHT.R)" - out = join(["TransRadial", grid, samples, nr, R, resp], "\n ") + parts = ["TransRadial", grid, samples, nr, R, resp] + if !isempty(t.sp_resp) + sp = "spectral responses: "*join([string(typeof(ri)) for ri in t.sp_resp], "\n ") + push!(parts, sp) + end + out = join(parts, "\n ") print(io, out) end @@ -433,7 +510,14 @@ function TransRadial(TT, grid, HT, FT, responses, densityfun, normfun) Pto = similar(Eto) Pωo = similar(Eωo) idcs = CartesianIndices(size(Pto)[2:end]) - TransRadial(HT, FT, normfun, responses, grid, densityfun, Pto, Eto, Eωo, Pωo, idcs) + # Split responses into time-domain and spectral (frequency-domain) + td_resp = Tuple(r for r in responses if !(r isa SpectralResponse)) + sp_resp = Tuple(r for r in responses if r isa SpectralResponse) + A2A_kω = isempty(sp_resp) ? + zeros(ComplexF64, (0, 0)) : + zeros(ComplexF64, (length(grid.ω), HT.N)) + TransRadial(HT, FT, normfun, td_resp, grid, densityfun, Pto, Eto, Eωo, Pωo, idcs, + sp_resp, A2A_kω) end """ @@ -470,6 +554,14 @@ function (t::TransRadial)(nl, Eω, z) @. t.Pto *= t.grid.towin # apodisation mul!(t.Pto, t.QDHT, t.Pto) # transform r -> k to_freq!(nl, t.Pωo, t.Pto, t.FT) # transform t -> ω + # Spectral responses (e.g. TPA): Eto still holds (t, r) field from above + if !isempty(t.sp_resp) + @. t.Pto = abs2(t.Eto) * t.Eto # |E|²E in (t, r) — reuse Pto buffer + @. t.Pto *= t.grid.towin + mul!(t.Pto, t.QDHT, t.Pto) # transform r -> k + to_freq!(t.A2A_kω, t.Pωo, t.Pto, t.FT) # transform t -> ω + sp_Et_to_Pω!(nl, t.A2A_kω, t.sp_resp, t.densityfun(z), t.idcs) + end nl .*= t.grid.ωwin .* (-im.*t.grid.ω)./(2 .* t.normfun(z)) end @@ -523,7 +615,7 @@ function norm_radial(grid, q, nfun) return norm end -mutable struct TransFree{TT, FTT, nT, rT, gT, xygT, dT, iT} +mutable struct TransFree{TT, FTT, nT, rT, gT, xygT, dT, iT, spT} FT::FTT # 3D Fourier transform (space to k-space and time to frequency) normfun::nT # Function which returns normalisation factor resp::rT # nonlinear responses (tuple of callables) @@ -536,6 +628,8 @@ mutable struct TransFree{TT, FTT, nT, rT, gT, xygT, dT, iT} Pωo::Array{ComplexF64, 3} # buffer for oversampled frequency-domain NL polarisation scale::Float64 # scale factor to be applied during oversampling idcs::iT # iterating over these slices Eto/Pto into Vectors, one at each position + sp_resp::spT # spectral responses (e.g. TPA), applied in ω domain before norm + A2A_kω::Array{ComplexF64, 3} # scratch buffer for F{|E|²E}(ω,ky,kx) end function show(io::IO, t::TransFree) @@ -544,7 +638,12 @@ function show(io::IO, t::TransFree) resp = "responses: "*join([string(typeof(ri)) for ri in t.resp], "\n ") y = "y grid: $(minimum(t.xygrid.y)) to $(maximum(t.xygrid.y)), N=$(length(t.xygrid.y))" x = "x grid: $(minimum(t.xygrid.x)) to $(maximum(t.xygrid.x)), N=$(length(t.xygrid.x))" - out = join(["TransFree", grid, samples, y, x, resp], "\n ") + parts = ["TransFree", grid, samples, y, x, resp] + if !isempty(t.sp_resp) + sp = "spectral responses: "*join([string(typeof(ri)) for ri in t.sp_resp], "\n ") + push!(parts, sp) + end + out = join(parts, "\n ") print(io, out) end @@ -556,8 +655,14 @@ function TransFree(TT, scale, grid, xygrid, FT, responses, densityfun, normfun) Pto = similar(Eto) Pωo = similar(Eωo) idcs = CartesianIndices((Ny, Nx)) - TransFree(FT, normfun, responses, grid, xygrid, densityfun, - Pto, Eto, Eωo, Pωo, scale, idcs) + # Split responses into time-domain and spectral (frequency-domain) + td_resp = Tuple(r for r in responses if !(r isa SpectralResponse)) + sp_resp = Tuple(r for r in responses if r isa SpectralResponse) + A2A_kω = isempty(sp_resp) ? + zeros(ComplexF64, (0, 0, 0)) : + zeros(ComplexF64, (length(grid.ω), Ny, Nx)) + TransFree(FT, normfun, td_resp, grid, xygrid, densityfun, + Pto, Eto, Eωo, Pωo, scale, idcs, sp_resp, A2A_kω) end """ @@ -601,6 +706,15 @@ function (t::TransFree)(nl, Eωk, z) @. t.Pto *= t.grid.towin # apodisation mul!(t.Pωo, t.FT, t.Pto) # transform (t, y, x) -> (ω, ky, kx) copy_scale!(nl, t.Pωo, length(t.grid.ω), 1/t.scale) + # Spectral responses (e.g. TPA): add to nl before inline norm, so they go + # through the same normalisation as time-domain responses. + if !isempty(t.sp_resp) + @. t.Pto = abs2(t.Eto) * t.Eto + @. t.Pto *= t.grid.towin + mul!(t.Pωo, t.FT, t.Pto) # transform (t, y, x) -> (ω, ky, kx) + copy_scale!(t.A2A_kω, t.Pωo, length(t.grid.ω), 1/t.scale) + sp_Et_to_Pω!(nl, t.A2A_kω, t.sp_resp, t.densityfun(z), t.idcs) + end nl .*= t.grid.ωwin .* (-im.*t.grid.ω)./(2 .* t.normfun(z)) end diff --git a/src/PhysData.jl b/src/PhysData.jl index f05372f4..9d1d9ca4 100644 --- a/src/PhysData.jl +++ b/src/PhysData.jl @@ -686,6 +686,85 @@ function n2_glass(material::Symbol; λ=nothing) end end +""" + TPA_params + +Material parameters for the empirical two-photon absorption coefficient model: + + β₂(ω) = C × (2ℏω/eV − E_g)^α [m/W] + +for 2ℏω > E_g (above the two-photon absorption edge), zero otherwise. + +### Fields per entry +- `E_g`: optical bandgap in eV +- `C`: pre-factor in m/W/eV^α +- `α`: power-law exponent + +### Fitting procedure +Fit `log(β₂)` vs `log(2ℏω/eV − E_g)` to measured data. + +### Units +- Input ω: rad/s +- Output β₂: m/W (SI) +- E_g in eV; C in m/W/eV^α + +### References +- Patankar et al., Appl. Opt. 56, 8309 (2017) (https://doi.org/10.1364/AO.56.008309) +- Nikogosyan et al., Opt. Commun. 200, 109 (2001) +- Van Stryland et al., Opt. Lett. 10, 490 (1985) +- Taylor A J, Gibson R B and Roberts J P 1988 Opt. Lett. 13 814 +- Kittelmann O and Ringling J 1994 Opt. Lett. 19 2053 +""" +const TPA_params = Dict( + :SiO2 => (E_g = 8.3, C = 1e-13, α = 3.5), # Nikogosyan, Taylor, Patankar, Kittelmann + :MgF2 => (E_g = 10.8, C = 8.4e-14, α = 4.35), # Kittelmann, Patankar +) + +""" + β₂_TPA(ω, material::Symbol; E_g=nothing, C=nothing, α=nothing) → Float64 + +Return the degenerate two-photon absorption coefficient β₂ (m/W) at +angular frequency `ω` (rad/s) for the given solid `material`. + +Uses the empirical power-law model: + + β₂(ω) = C × (2ℏω/eV − E_g)^α × Θ(2ℏω − E_g) + +where E_g is the material bandgap in eV, C is the pre-factor (m/W/eV^α), +α is the power-law exponent, and Θ is the Heaviside step function. +Returns 0.0 for frequencies below the two-photon edge (2ℏω ≤ E_g). + +# Arguments +- `ω`: angular frequency in rad/s (must be positive) +- `material`: one of `:SiO2`, `:MgF2` + +# Keyword arguments (override material defaults) +- `E_g`: bandgap energy in eV +- `C`: pre-factor in m/W/eV^α +- `α`: power-law exponent (dimensionless) +""" +function β₂_TPA(ω::Real, material::Symbol; + E_g=nothing, C=nothing, α=nothing) + haskey(TPA_params, material) || error( + "Unknown material for TPA: $material. Supported: $(keys(TPA_params))") + p = TPA_params[material] + E_g_eV = isnothing(E_g) ? p.E_g : E_g + C_val = isnothing(C) ? p.C : C + α_val = isnothing(α) ? p.α : α + + # Two-photon energy in eV + two_photon_eV = 2 * ħ * abs(ω) / electron + + excess_eV = two_photon_eV - E_g_eV + excess_eV <= 0 && return 0.0 + + isnan(C_val) && error( + "TPA parameters for :$material have not been fitted yet. " * + "See TPA_params docstring for fitting procedure.") + + return C_val * excess_eV^α_val +end + """ density(material::Symbol, P=1.0, T=roomtemp) diff --git a/test/runtests.jl b/test/runtests.jl index a2552f24..d6f237bb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -157,4 +157,9 @@ end include(joinpath(testdir, "test_gnlse.jl")) end +@testset "TPA" begin + @info("================= test_tpa.jl") + include(joinpath(testdir, "test_tpa.jl")) +end + end \ No newline at end of file diff --git a/test/test_tpa.jl b/test/test_tpa.jl new file mode 100644 index 00000000..10010964 --- /dev/null +++ b/test/test_tpa.jl @@ -0,0 +1,231 @@ +using Luna +import Test: @test, @testset, @test_throws +import Luna: PhysData, Nonlinear, NonlinearRHS, Grid, Processing +import Logging + +logger = Logging.SimpleLogger(stdout, Logging.Warn) +old_logger = Logging.global_logger(logger) + +@testset "β₂_TPA coefficient values" begin + # 300 nm → two-photon energy = 2 * 1240/300 = 8.27 eV < E_g = 8.3 eV → below edge + ω_300 = 2π * PhysData.c / 300e-9 + @test PhysData.β₂_TPA(ω_300, :SiO2) == 0.0 + + # 264 nm → 2ℏω = 2 * 1240/264 ≈ 9.39 eV > 8.3 eV → nonzero + ω_264 = 2π * PhysData.c / 264e-9 + β₂_264 = PhysData.β₂_TPA(ω_264, :SiO2) + @test β₂_264 > 0 + @test 1e-16 < β₂_264 < 1e-12 # order of magnitude check + + # 213 nm → 2ℏω ≈ 11.6 eV → much stronger + ω_213 = 2π * PhysData.c / 213e-9 + β₂_213 = PhysData.β₂_TPA(ω_213, :SiO2) + @test β₂_213 > β₂_264 # monotonicity +end + +@testset "TPAResponse coefficient conversion" begin + Nω = 100 + ω_grid = range(0.5e15, 2.0e15, length=Nω) + # Synthetic β₂ values: nonzero above some threshold + β₂_ω = [ω > 1.0e15 ? 1e-12 : 0.0 for ω in ω_grid] + + tpa = Nonlinear.TPAResponse(collect(ω_grid), β₂_ω) + + # coeff should be purely imaginary (negative imaginary part) where β₂ > 0 + for i in eachindex(ω_grid) + if β₂_ω[i] > 0 + @test real(tpa.coeff_ω[i]) ≈ 0.0 atol=1e-40 + @test imag(tpa.coeff_ω[i]) < 0 # -i gives negative imaginary + else + @test tpa.coeff_ω[i] == 0.0 + end + end + + # Check magnitude: coeff = -i * ε₀² * c² * β₂ / (2ω) + idx = findfirst(x -> x > 0, β₂_ω) + expected = -im * PhysData.ε_0^2 * PhysData.c^2 * β₂_ω[idx] / (2 * ω_grid[idx]) + @test tpa.coeff_ω[idx] ≈ expected + + # Callable: adds to output buffer + out_ω = zeros(ComplexF64, Nω) + F_E2E = ones(ComplexF64, Nω) + ρ = 1.0 + tpa(out_ω, F_E2E, ρ) + @test all(out_ω .== tpa.coeff_ω) + + # DimensionMismatch + @test_throws DimensionMismatch Nonlinear.TPAResponse([1.0, 2.0], [1.0]) +end + +@testset "TPA energy loss — mode-averaged" begin + # Propagate a pulse through a thin medium with TPA only (n2≈0). + # TPA should cause energy loss. + λ0 = 240e-9 + τfwhm = 5e-15 + energy = 20e-9 # 1 nJ + flength = 100e-6 # 100 µm + + β2 = 1e-30 # small normal dispersion + βs = [0.0, 0.0, β2] + Aeff = π * (12e-6)^2 # realistic effective area + n2_val = 1e-25 # very small n2 (effectively no Kerr) + λlims = [160e-9, 600e-9] + trange = 1e-12 + + # Synthetic TPA: constant β₂_TPA = 1e-11 m/W across band (strong enough to see effect) + β₂_tpa_val = 1e-11 # ~1 cm/GW + tpa_func = ω -> β₂_tpa_val + + # With TPA + output_tpa = prop_gnlse(flength, βs; n2=n2_val, Aeff, λ0, τfwhm, energy, + pulseshape=:gauss, λlims, trange, + raman=false, shock=false, shotnoise=false, tpa=tpa_func) + + # Without TPA + output_notpa = prop_gnlse(flength, βs; n2=n2_val, Aeff, λ0, τfwhm, energy, + pulseshape=:gauss, λlims, trange, + raman=false, shock=false, shotnoise=false) + + # Energy comparison + Eω_in_tpa = output_tpa["Eω"][:, 1] + Eω_out_tpa = output_tpa["Eω"][:, end] + energy_in = sum(abs2.(Eω_in_tpa)) + energy_out_tpa = sum(abs2.(Eω_out_tpa)) + + Eω_out_notpa = output_notpa["Eω"][:, end] + energy_out_notpa = sum(abs2.(Eω_out_notpa)) + + # Without TPA, energy should be conserved (no loss, no Raman) + @test isapprox(energy_out_notpa, sum(abs2.(output_notpa["Eω"][:, 1])), rtol=1e-13) + + # With TPA, output energy should be less than input + @test energy_out_tpa < energy_in + # And less than without TPA + @test energy_out_tpa < energy_out_notpa + + # Fractional loss should be positive and reasonable (not 100%) + frac_loss = 1 - energy_out_tpa / energy_in + @test 0.7 < frac_loss < 1 +end + +@testset "TPA spectral suppression" begin + # Broadband pulse: TPA with frequency-dependent β₂ should suppress + # short-wavelength components more than long-wavelength ones. + λ0 = 240e-9 + τfwhm = 5e-15 # short pulse = broad bandwidth + energy = 20e-9 + flength = 100e-6 + β2 = 1e-30 + βs = [0.0, 0.0, β2] + Aeff = π * (12e-6)^2 + n2_val = PhysData.n2(:SiO2) + λlims = [160e-9, 600e-9] + trange = 1e-12 + + tpa_func = ω -> PhysData.β₂_TPA(ω, :SiO2) + + output_tpa = prop_gnlse(flength, βs; n2=n2_val, Aeff, λ0, τfwhm, energy, + pulseshape=:gauss, λlims, trange, + raman=false, shock=false, shotnoise=false, tpa=tpa_func) + + output_notpa = prop_gnlse(flength, βs; n2=n2_val, Aeff, λ0, τfwhm, energy, + pulseshape=:gauss, λlims, trange, + raman=false, shock=false, shotnoise=false) + + ω = output_tpa["grid"]["ω"] + spec_tpa = abs2.(output_tpa["Eω"][:, end]) + spec_notpa = abs2.(output_notpa["Eω"][:, end]) + + # Find spectral ratio (TPA / no-TPA) at two frequencies: + ω_below = 2π * PhysData.c / 250e-9 # less TPA + ω_above = 2π * PhysData.c / 220e-9 # more TPA + + idx_below = argmin(abs.(ω .- ω_below)) + idx_above = argmin(abs.(ω .- ω_above)) + + # Below edge: spectrum should be essentially unchanged + if spec_notpa[idx_below] > 1e-30 # only test if there's signal + ratio_below = spec_tpa[idx_below] / spec_notpa[idx_below] + ratio_above = spec_tpa[idx_above] / spec_notpa[idx_above] + # Above-edge suppression should be stronger (lower ratio) + @test ratio_above < ratio_below + end +end + +@testset "TPA spectral suppression — radial free-space" begin + import Luna: Grid, Nonlinear, NonlinearRHS, LinearOps, Output, Fields + import Luna.PhysData: wlfreq + import FFTW + import Hankel + + # Deep-UV Gaussian beam in SiO₂ — envelope, radially symmetric + λ0 = 240e-9 + τfwhm = 2e-15 # broad bandwidth across TPA region + energy = 150e-9 + w0 = 15e-6 # beam waist + L = 10e-6 # 100 µm propagation (thin slab) + R = 250e-6 # Hankel aperture + N = 256 # radial points + + grid = Grid.EnvGrid(L, λ0, (160e-9, 600e-9), 1e-12) + q = Hankel.QDHT(R, N, dim=2) + + densityfun(z) = 1.0 # solid — χ³ already bulk value + nfun = PhysData.ref_index_fun(:SiO2, lookup=false) + normfun = NonlinearRHS.const_norm_radial(grid, q, nfun) + linop = LinearOps.make_const_linop(grid, q, nfun) + + χ3_SiO2 = PhysData.χ3(:SiO2) + β₂_ω = PhysData.β₂_TPA.(grid.ω, :SiO2) + tpa = Nonlinear.TPAResponse(grid.ω, β₂_ω) + + inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τfwhm, energy=energy, w0=w0) + + # With TPA + responses_tpa = (Nonlinear.Kerr_env(χ3_SiO2), tpa) + Eω_tpa, transform_tpa, FT_tpa = Luna.setup(grid, q, densityfun, normfun, + responses_tpa, inputs) + output_tpa = Output.MemoryOutput(0, grid.zmax, 201) + Luna.run(Eω_tpa, grid, linop, transform_tpa, FT_tpa, output_tpa) + + # Without TPA + responses_notpa = (Nonlinear.Kerr_env(χ3_SiO2),) + Eω_notpa, transform_notpa, FT_notpa = Luna.setup(grid, q, densityfun, normfun, + responses_notpa, inputs) + output_notpa = Output.MemoryOutput(0, grid.zmax, 201) + Luna.run(Eω_notpa, grid, linop, transform_notpa, FT_notpa, output_notpa) + + # Energy comparison (sum over all radial k-modes) + Eωk_tpa_in = output_tpa.data["Eω"][:, :, 1] + Eωk_tpa_out = output_tpa.data["Eω"][:, :, end] + Eωk_notpa_out = output_notpa.data["Eω"][:, :, end] + + energy_in = sum(abs2, Eωk_tpa_in) + energy_out_tpa = sum(abs2, Eωk_tpa_out) + energy_out_notpa = sum(abs2, Eωk_notpa_out) + + # With TPA, output energy should be less than without + @test energy_out_tpa < energy_out_notpa + + # Spectral suppression: TPA stronger at shorter wavelengths (higher ω) + # Sum over radial modes for spectral comparison + spec_tpa = dropdims(sum(abs2, Eωk_tpa_out, dims=2), dims=2) + spec_notpa = dropdims(sum(abs2, Eωk_notpa_out, dims=2), dims=2) + + ω_below = 2π * PhysData.c / 300e-9 # below SiO₂ TPA edge (~299 nm) + ω_above = 2π * PhysData.c / 220e-9 # well above TPA edge + + idx_below = argmin(abs.(grid.ω .- ω_below)) + idx_above = argmin(abs.(grid.ω .- ω_above)) + + if spec_notpa[idx_below] > 1e-30 && spec_notpa[idx_above] > 1e-30 + ratio_below = spec_tpa[idx_below] / spec_notpa[idx_below] + ratio_above = spec_tpa[idx_above] / spec_notpa[idx_above] + # Above-edge suppression should be stronger (lower ratio) + @test ratio_above < ratio_below + @test ratio_above < 0.6 # significant suppression above edge + end +end + +## +Logging.global_logger(old_logger)