WIP: Implement two-photon absorption (TPA) support

.gitignore

@@ -10,3 +10,6 @@ _git2_*
 docs/build
 *.info
 .vscode/spellright.dict
+.DS_Store
+.gitignore
+deps/build.log

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
+
+## 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)

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)

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)