The forward model

The forward model predicts the trace that a trial pulse would produce in a specified instrument. Choosing it is a physical decision: include propagation, beam geometry and collection when they materially affect the measurement. The dispersive tutorial illustrates the consequence of omitting propagation; validation examines more complete instrument simulations.

Start with the thin-medium model below, then add the effects justified by your experiment. ForwardModel is the NumPy core: it maps a complex spectrum \(\tilde E(\omega)\) to the trace for the thin and dispersive-slab models and provides the exact reverse-mode (Wirtinger) adjoint of that map, which COPRA and the analytic-gradient solvers descend on. It also evaluates the reduced smearing trace (a smearing= kernel), but without an adjoint. The thin medium is just the slab with \(\beta=0\), \(L=0\) and a single quadrature node (quadrature_nodes_weights() returns one node for zero thickness), so one code path and one adjoint serve both regimes. The JAX model (make_param_trace_fn()) implements the same models plus the spatial extensions below, differentiated by automatic differentiation; the two paths are cross-checked where they overlap.

The unknown field must have a declared reference plane. In the slab model it is the entrance field; the chromatic model further distinguishes an on-axis field from a transverse-integrated or effective field. See field reference conventions.

Thin medium

In the thin, non-dispersive limit the signal for one delay \(\tau\) is

\[\psi(\omega, \tau) \;=\; \mathcal{F}\!\left\{\, s\bigl(E(t),\, E(t-\tau)\bigr)\right\}(\omega), \qquad T(\omega, \tau) = \lvert \psi(\omega, \tau)\rvert^2,\]

where \(s(\cdot,\cdot)\) is the interaction signal, \(E(t)\) is the test field and \(E(t-\tau)\) the gate. Synthesising a full trace is a handful of FFTs per delay. This is the default — material=None.

import numpy as np, croak
g = croak.Grid(128, dt=0.3e-15)
ew = croak.gaussian_pulse(g, 2.5e-15)
delays = np.linspace(-18e-15, 18e-15, 87)

trace = croak.maketrace(g.omega, delays, ew, "shg")    # (Nomega, Ndelay)

maketrace() is a convenience wrapper; the underlying object is ForwardModel, whose trace() and signal_single() you can call directly.

Dispersive slab

To capture propagation inside the nonlinear medium, croak integrates the signal generation over the slab depth. A trial spectrum defines the input replicas at the entrance face; each propagates to depth \(z\) with the linear propagator \(e^{i\beta(\omega) z}\); the nonlinear signal is generated locally; and it propagates the remaining distance \(L-z\) to the exit face, accumulating its own phase. The exit signal is the coherent superposition over the thickness.

Discretised with Gauss–Legendre nodes \(z_q\) and weights \(w_q\) on \([0, L]\), the per-delay signal that ForwardModel computes is

(1)\[\psi(\omega,\tau) \;=\; \sum_q w_q\, e^{i\beta(\omega)(L-z_q)} \,\odot\, \mathcal{F}^{-1}\!\Bigl[\, s\bigl( \mathcal{F}[\tilde E\, e^{i\beta z_q}],\; \mathcal{F}[\tilde E\, e^{i\omega\tau}\, e^{i\beta z_q}] \bigr)\Bigr],\]

with the trace again \(T = \lvert\psi\rvert^2\). Here \(\mathcal{F}\) is the raw DFT; the model works internally in DFT-bin order and exposes centred arrays at its boundary (see FFT convention). Setting \(\beta=0\), \(L=0\) and a single node recovers the thin-medium formula exactly.

The propagation constant

\(\beta(\omega)\) comes from croak.materials.beta(). On the centred grid (relative to the carrier \(\omega_0\)):

\[\beta(\omega) \;=\; \frac{(\omega+\omega_0)\,n(\omega+\omega_0)}{c} \;-\; \beta_1(\omega_0)\,\omega,\]

where \(n\) is the material’s refractive index and the group-delay term \(\beta_1(\omega_0)\,\omega\) keeps the pulse centred in the time window. The built-in materials are listed in croak.materials.MATERIALS (SiO₂, BK7, CaF₂, BaF₂, MgF₂, sapphire and silicon, each optmat’s recommended model of the material); a model is extrapolated beyond its stated range, and the entries at non-positive frequency, which have no index, are set to zero. You are not limited to these — register a custom n(λ) and pass its name here just the same; see Materials and chirped mirrors.

Using it

Pass material, thickness, npoints and the carrier omega0:

omega0 = croak.maths.wlfreq(260e-9)        # 260 nm
trace = croak.maketrace(g.omega, delays, ew, "pg",
                       material="SiO2", thickness=10e-6, npoints=20, omega0=omega0)

Remember that dispersive propagation is supported for PG and SD only (their signal stays at the fundamental); SHG raises a ValueError. See Interactions.

Quadrature convergence

quadrature_nodes_weights() builds the depth integral. Gauss–Legendre quadrature is efficient for a smooth depth integrand; the "uniform" rule provides a useful independent discretisation check. Required node counts depend on bandwidth, thickness, dispersion and field structure. The solvers default to npoints=20, chosen for few-µm to tens-of-µm UV slabs. For a thicker or more dispersive medium, increase npoints until the trace and the retrieved field stop changing.

Note

The GUI ties npoints to a depth step Δz. The Retrieve stage’s Dispersive propagation group exposes a Δz step (µm) control (default 1 µm) alongside the thickness and Npoints, with npoints = round(thickness / Δz). Editing the thickness holds Δz fixed and grows Npoints so the node density stays constant; editing Npoints back-computes Δz to the actual spacing. Only Npoints reaches quadrature_nodes_weights() — Δz is a convenience for picking it, so a 40 µm slab automatically gets ~4× the nodes of a 10 µm one. The nodes are the (non-uniform) Gauss–Legendre abscissae, so Δz is a nominal step, not the literal spacing of every node.

Generation envelope (thick slabs / tight focusing)

The depth integral above assumes uniform generation density: every depth contributes with its quadrature weight, modified only by the physics the 1D model contains (the replicas disperse, so deeper slices are weaker). What a 1D model cannot see is the transverse geometry of the crossing beams: the focal spot growing away from the focus (Rayleigh range), the beamlets walking through each other, and the Gouy phase — which enters the coherent depth sum as a z-dependent phase, not just an amplitude. To good approximation all of this factorises into one smooth complex envelope \(g(z)\) multiplying the local generation.

depth_weight on make_param_trace_fn() (and the lbfgs-ad solver) accepts one complex weight per quadrature node and multiplies the depth-quadrature weights with it. The overall scale is irrelevant (absorbed by \(\mu\)); only the shape matters. Two rules of thumb:

• Negligible when \(L \lesssim z_R\) — which covers conventional 10–50 µm practice with ordinary focusing (for the example DUV geometry, \(z_R \approx 87\) µm: the envelope varies by <10 % across a 40 µm slab, and its smooth symmetric part is nearly pure renormalisation).

• Required once \(L\) approaches the beam-overlap length \(\sim w_0/\theta\), with \(\theta\) the beam crossing angle (the crossed beamlets walk through each other): beyond it, generation stops regardless of how much substrate remains, and a uniform-weight model would attribute signal to depths that produce none. A diagnostic signature: a uniform-weight retrieval with fit_thickness returns the effective generation length, not the physical slab.

With fit_thickness the envelope is pinned to the normalised nodes, i.e., it rescales with the fitted slab rather than staying fixed in absolute \(z\).

The adjoint

adjoint_single() is the exact reverse-mode adjoint of (1). Given the trace cotangent it returns \(\partial L/\partial\tilde E^*\) — the gradient the solvers descend on — at the cost of a single extra pass that reuses the recorded forward fields. This is what lets COPRA and L-BFGS run in pure NumPy, with no JAX on the hot path. The Gradients page derives it, sets out what it covers, and explains how it is cross-checked against automatic differentiation.

Geometric smearing

The equations above evaluate the polarisation at a single transverse point. In a non-collinear BOXCARS geometry that is not quite what a spectrometer sees: the three input beams arrive at the interaction plane with tilted pulse fronts, so their relative arrival times vary linearly across the focal spot, and the detector integrates over that spread. croak.smearing puts that instrument response inside the forward model.

A beam apertured at mask position \(\mathbf r_j\) and focused by a lens of focal length \(f\) crosses the focus at angle \(\mathbf r_j/f\), so at focal-plane position \(\mathbf r\) it arrives at \(t - \boldsymbol\alpha_j\cdot\mathbf r\) with \(\boldsymbol\alpha_j = -\mathbf r_j/(fc)\). Because \(|\mathcal F\{\cdot\}|^2\) is blind to a global time origin — and that origin may be chosen per \(\mathbf r\) — only two scalar combinations of the three arrival times survive:

(2)\[\begin{split}\mathrm{PG:}\quad S &= E(t)\,E(t-\tau+\delta+p/2)\,E^*(t-\tau+\delta-p/2) \\ \mathrm{SD:}\quad S &= E(t+p/2)\,E(t-p/2)\,E^*(t-\tau+\delta)\end{split}\]

Setting \(p=\delta=0\) recovers test * |gate|**2 and test**2 * conj(gate) exactly, so switching the kernel off changes nothing.

\(p\) separates the two replicas that coincide in the 1-D model, so it changes the shape of the signal and needs an explicit quadrature. \(\delta\) is a pure translation of the delay axis, so integrating over it is a convolution along \(\tau\), applied once to the finished trace. Both are linear in \(\mathbf r\), so with an isotropic Gaussian focal weight the pair \((p,\delta)\) is jointly Gaussian — which is all SmearingKernel stores.

Kernel width from the mask

The local signal amplitude is the product of the three focal-plane amplitudes, so position \(\mathbf r\) contributes with weight \(|A(\mathbf r)|^6\). For a uniformly illuminated circular hole of diameter \(D\) that is the Airy amplitude \(2J_1(u)/u\), \(u=\pi Dr/(\lambda f)\), whose rms per Cartesian component is \(\sigma_r = C\lambda f/D\) with \(C = 0.245673\) (AIRY6_RMS_COEFF). Since \(\boldsymbol\alpha_j \propto 1/f\) the focal length cancels — the crossing angle and the spot size scale oppositely. For the standard folded square BOXCARS with holes at \((\pm d,\pm d)\):

interaction	\(\sigma_p\)	\(\sigma_\delta\)	\(\rho\)
"pg" (delay on the gate pair)	\(0.49135\,(d/D)\,\lambda/c\)	\(0.54933\,(d/D)\,\lambda/c\)	\(1/\sqrt5\)
"sd" (delay on the conjugated arm)	\(0.69487\,(d/D)\,\lambda/c\)	\(0.34743\,(d/D)\,\lambda/c\)	\(0\)

Within this achromatic, full-collection approximation, the widths depend on the mask ratio \(d/D\) and wavelength; more complete models also depend on collection, focusing and depth evolution. The two rows differ because the arm that carries the scanned delay decides which pair of replicas \(p\) separates; only in the SD arrangement do \(p\) and \(\delta\) come out uncorrelated. Do not assume one geometry’s ratios carry over to another — kernel_from_arms() derives them from the actual hole coordinates.

For an example DUV geometry (1 mm holes, 0.5 mm edge to edge so \(d/D=0.75\), 260 nm) the PG kernel is \(\sigma_p = 0.320\) fs, \(\sigma_\delta = 0.357\) fs — a 0.84 fs FWHM blur along the delay axis, which matters for a 1 fs pulse. The simulated instrument used in validation has \(d/D=1.0\), with proportionally wider kernels.

Cost and accuracy

The kernel multiplies the depth quadrature rather than coupling to it: the pulse-front tilts are fixed before the medium, and the longitudinal walk-off inside a thin slab is far below a femtosecond, so the same \((p,\delta)\) applies at every depth node. (The depth-independence of the transverse phase is a separate approximation, and exactly the one the collected model’s depth_transverse flag relaxes — see the subsection below.) Because the quadrature nodes are exactly antisymmetric, each shifted field is built once and used twice, and the cost is about \((K+1)/2\) times the unsmeared model — roughly 3× at the default \(K=5\) nodes. The \(\delta\) integral is free: it is a multiplication of the delay-axis DFT by the analytic Gaussian transform, which is exact for any trace that is Nyquist-sampled in \(\tau\) (so the kernel need not span several delay steps). It does require a uniformly spaced delay axis.

Within the achromatic, fully collected reduction, the remaining approximation is replacing the exact \(|A|^6\) focal weight by a Gaussian with the same second moments. The leading trace correction is quadratic in \((p,\delta)\), so only those moments matter to first order. For the chirped test pulse and reference geometry in tests/test_smearing.py, the Gaussian and Airy-weighted focal-plane traces differ by about one part in \(10^5\) of the trace peak.

Note

Geometric smearing is supported through the JAX model by warm-lbfgs, lbfgs-ad, lm, lm-optx, lbfgs-optx and cma-es. The incoherent sum over quadrature nodes means the measured trace is no longer \(|\psi(\omega,\tau)|^2\) for a single signal, which breaks COPRA’s magnitude-replacement projection, and the hand-derived adjoint of (1) is not extended to it. Asking for smearing with any other algorithm raises rather than silently dropping the kernel.

See Geometric smearing for how to use it.

Beyond the reduction: the mixture and the collection aperture

Two assumptions in the paragraphs above are separable, and croak lets each be dropped.

croak.focal drops the achromatic one. The reduction freezes the \(|A|^6\) weight at the carrier, but the mask maps to transverse wavevector as \(x = k_\perp z_\text{mask} c/\omega_\text{abs}\), so \(A\) depends on the product \(r\omega_\text{abs}\) and different focal positions see different spectra. The mixture keeps the sum over \(\mathbf r\) explicit, each node applying its own \(A(r,\omega)\) and carrying its own definite \((p,\delta)\) instead of a distribution over them. Held at the carrier it reproduces the closed forms above exactly, which reduced_moments() checks.

croak.collection drops the full-beam collection one. “\(|\mathcal F\{\cdot\}|^2\) is blind to a global time origin, and that origin may be chosen per \(\mathbf r\)” is the step that produced \((p,\delta)\), and it is legitimate only because Parseval turns \(\int|\tilde S(\mathbf k)|^2\,\mathrm d^2k\) over all \(\mathbf k\) into \(\int|S(\mathbf r)|^2\,\mathrm d^2 r\). Through a finite hole the signal must be transformed before it is squared, and the discarded per-position phase comes back.

The reduction itself survives that intact. The exact focal-plane signal and croak’s per-node build differ by a rigid time shift, \(S(t,\mathbf r) = S_\text{croak}(t-\boldsymbol\alpha_1\cdot\mathbf r,\mathbf r)\), so the only thing lost is one phase — and after the signal’s own phase-matched carrier \(\boldsymbol\alpha_s = \boldsymbol\alpha_1+\boldsymbol\alpha_2-\boldsymbol\alpha_3\) is removed, what remains is \(e^{+i\Omega p}\), with \(\Omega=\omega_\text{abs}-\omega_0\) the centred frequency: exactly \(\Omega\) times one of the two reduced parameters (and \(-\Omega\delta\) for SD). The \((p,\delta)\) parameterisation is therefore sufficient for an apertured instrument as well as a fully collected one.

Whether it matters is a question about one ratio: the aperture’s k half-width \((\omega_\text{abs}/c)(D_\text{hole}/2)/z_\text{mask}\) against the signal’s own k content \(3(\omega_\text{abs}/c)(D/2)/f\), both at absolute optical frequency. For the example DUV geometry that is 6.04e4 against 3.63e5 m⁻¹ — the aperture is six times narrower, so the fully collected (incoherent) reduction does not describe that instrument and the collection must be modelled. See Modelling the collection aperture.

Depth-resolved transverse dephasing (depth_transverse)

One assumption survives even the collected model: every depth node of (1) is added with the same transverse phase — fully coherently in transverse wavevector. Real propagation dephases the angular spectrum. Signal generated at depth \(z\) accumulates, over the remaining slab, the extra phase

\[ \phi(\omega, k_\perp, z) = \bigl[k_z(\omega, k_\perp) - k_z(\omega, 0)\bigr](L - z) \;\approx\; -\frac{c\,k_\perp^2\,(L - z)}{2\,n(\omega)\,\omega}, \qquad k_z = \sqrt{(n\omega/c)^2 - k_\perp^2}, \]

relative to the on-axis propagation the model already applies — depth-accumulating, quadratic in \(k_\perp\), and chromatic (\(\sim 1/\omega\), red-weighted). The depth_transverse flag of make_param_trace_fn() applies it inside the collection transform: the depth sum moves inside the aperture integral, each depth node’s transformed field carrying \(e^{i\phi_j(\omega,z_q)}\) at the aperture node’s absolute wavevector before the coherent sum over depth. It is the k-dependent sibling of the generation envelope above: depth_weight carries the k-independent per-depth amplitude factor \(g(z)\), depth_transverse the k-dependent per-depth phase, and the two compose. Parameter-free, opt-in, and validated (sign included) against a per-depth dense-FFT reference with the exact \(k_z\) — see the collection how-to for usage, cost and caveats.

The depth-resolved mixture (fc-z, evolve_profiles)

depth_transverse repairs the exit side (how the generated signal reaches the aperture); the depth-resolved mixture repairs the entrance side, by allowing depth evolution: the chromatic beamlet profiles are not frozen at the entrance face but evolve linearly through the slab, and each of the three arms evolves about its own walked axis. With focal_mixture(..., evolve_profiles=True, material=..., thickness=..., npoints=...) every arm \(j\) at node \(\mathbf r_k\) and depth \(z_q\) applies its own complex profile

\[ A_j(r_k, \omega; z_q) = N(\omega)\!\int_0^{k_R(\omega)}\! S(k_\perp,\omega)\, J_0\!\bigl(k_\perp\,|\mathbf r_k - \bm\delta_j(z_q)|\bigr)\, e^{\,i[k_z(\omega,k_\perp)-k_z(\omega,0)]\,z_q}\, k_\perp\,\mathrm dk_\perp, \qquad \bm\delta_j(z) = -\frac{z\,\mathbf r_j}{f\,n_0}, \]

in place of the single shared entrance filter \(A(r_k,\omega)\): the aperture’s k-space (a disc of radius \(k_R = (\omega_\text{abs}/c)(D/2)/f\), or the Gaussian control’s \(e^{-k_\perp^2 w_0^2/4}\)) propagated by the transverse part of \(k_z\) only — the on-axis dispersion stays in the depth propagator, so \(z\to 0\) reduces exactly to the entrance model. The \((p_k,\delta_k)\) arrival offsets stay depth-independent (exact: transverse wavevector is conserved). Under fc-z the “Cost and accuracy” statement that the transverse weight commutes with the depth propagator no longer holds — deliberately: that non-commutation is what lets the model carry a gate-vs-probe arrival asymmetry that accumulates with depth. Cost is ~fc (the depth axis is already batched; only a seconds-level numpy filter build and a \((3,K,Q,N)\) table are new). PG only; fitted thickness is unsupported; the spec’s slab parameters must match the trace map’s (verified). Validated by an independent plane-wave propagator at \(10^{-6}\) of peak and by moment-pinning against a three-beam overlap/centroid/timing table measured on the 3-D instrument simulations — see tests/test_fcz.py.

The model hierarchy, and when each member applies

The options form a hierarchy, with different assumptions and cost:

Model	Included physics	Main assumptions to assess
Standard	Local instantaneous interaction	Propagation and spatial averaging negligible
Extended	Coherent slab generation and reduced timing kernel	Achromatic local profiles and full isolated-signal collection adequate
Full	Chromatic focal fields and coherent finite-aperture collection, with selected depth corrections	Specified beam profiles, paraxial geometry, material response and detector model adequate

The standard model is conventional FROG practice: right for thin media and multi-cycle pulses. Two dimensionless quantities decide whether more is needed: the fractional bandwidth (chromatic focal effects) and the ratio of collection aperture to generation aperture (coherent-collection effects). With the dispersive slab and smearing kernel, the extended model costs about the same as the standard one. For near-single-cycle pulses through a restrictive collection aperture, use the full model; opening the aperture reduces collection effects, but pump leakage limits how far it can be opened. In the simulated TG-FROG instrument the full model with the recorded 0.5 mm aperture recovers a 1.05 fs pulse to 1.047–1.079 fs at every thickness from 1 to 40 µm, while the standard model returns 1.35–2.41 fs; the validation chapter gives the details, and the gallery provides runnable illustrations.

API

croak.forward.ForwardModel, croak.forward.maketrace(), croak.forward.quadrature_nodes_weights(), croak.smearing, croak.focal, croak.collection and croak.materials in the API reference.