Simulation

Model

The field is a complex envelope \(A(x, y, t)\) with carrier \(\exp[i(k_0 z - \omega_0 t)]\), propagated with the unidirectional pulse propagation equation (UPPE) in the frame moving at the group velocity at \(\lambda_0\):

\[\frac{\partial \tilde A(k_x, k_y, \Omega)}{\partial z} = i\left[k_z - k_0 - k_1\Omega\right]\tilde A + i\,\gamma(\Omega)\,\mathcal{F}_t\!\left[|A|^2 A\right], \qquad k_z = \sqrt{k(\omega)^2 - k_x^2 - k_y^2},\qquad \gamma = \frac{n_2\,\omega}{c}\,\frac{n(\omega_0)}{n(\omega)} .\]

Linear part: diffraction and the full dispersion of the gas (Sellmeier, Börzsönyi et al., Appl. Opt. 47, 4856, 2008).
Kerr part: instantaneous electronic response, with self-steepening (self_steepening=False gives \(\gamma = n_2\omega_0/c\)). It is tapered to zero at the edges of the spectral window.
Not modelled: ionisation and plasma, Raman response of molecular gases (N₂, air), the folded off-axis geometry (incidence astigmatism, real spot pattern), plates, windows and pressure gradients, vector effects.

Two solvers share the same pulses, diagnostics and files:

	lens frame (recommended)	standard
classes	`suggest_lens_grid()`, `LensFrameSimulation`	`suggest_grid()`, `MPCSimulation`
grid	follows the beam: a few waists wide, whatever \(d\)	fixed: a few mirror spots wide
diffraction	paraxial	non-paraxial
use	every cell, near-concentric ones above all	cross-checks

The lens frame is explained in The lens frame. The standard solver keeps the field in \((x, y, \Omega)\): the linear step needs only 2D spatial FFTs and the Kerr term only 1D temporal FFTs. They are combined by Strang splitting, with an RK4 Kerr step of size \(h = \varphi_\text{max}/(k_0 n_2 I_\text{peak})\).

The cell

MPC is a symmetric two-mirror cell, treated unfolded: the run starts right after a mirror and each pass goes to the cell centre, on to the next mirror, and bounces.

mirror ──(d/2)──▶ cell centre ──(d/2)──▶ next mirror │ bounce

The pulse is defined at the cell centre, and the diagnostics are taken there after every pass. That plane is the waist of the eigenmode with a flat wavefront, so a Hermite–Gauss decomposition there is exact. A bounce applies the mirror curvature (focal length \(R/2\)), its reflectivity, GDD, TOD or measured phase, and an optional soft clear aperture.

With \(g = 1 - d/R\), the Gouy phase per pass is \(\theta = \arccos g\) and

\[w_0^2 = \frac{\lambda}{2\pi}\sqrt{d(2R-d)},\qquad w_\text{mirror}^2 = \frac{\lambda R}{\pi}\sqrt{\frac{d}{2R-d}} .\]

MPC.herriott(R, N, k, medium) builds the cell that is re-entrant after \(N\) passes, \(d = R\,[1\pm\cos(2\pi k/N)]\) (+: concentric branch). A linear re-entrant cell images any input onto the output, top-hats included: run it once with \(n_2 = 0\) and check that it does.

Nonlinear mode matching. For a Gaussian at constant power \(P\), the Kerr lens changes the eigenmode as if \(\lambda \to \lambda\sqrt{1 - P/P_c}\), with \(P_c = \lambda^2/(2\pi n_0 n_2)\): the matched waist and spot are the linear ones times \((1-P/P_c)^{1/4}\). cell.match(pulse) returns them with the B-integral per pass. A pulse is matched at one power only: its wings see a mode closer to the linear one, a mismatch that the 4D simulation resolves.

The pulse

A pulse is a spectrum times a spatial beam, \(A(x, y, \Omega) = E(\Omega)\,U(x, y; \Omega)\).

E = m.gaussian_spectrum(grid, fwhm_fs=350)                    # or fwhm_nm=...
E = m.measured_spectrum(grid, wl_nm, intensity, phase_rad=ph) # spectrometer (+ retrieved phase)
E = m.spectrum_from_temporal_field(grid, t_fs, field)         # retrieved pulse
E = m.add_spectral_phase(grid, E, gdd_fs2=500, tod_fs3=2e3)

beam = m.Beam(m.top_hat(4e-3),
              zernike={"astig_0": 0.05, "coma_x": 0.02},     # RMS waves at lambda0
              spatial_chirp_mm_per_nm=(0.01, 0), pft_fs_per_mm=(2.0, 0))
beam = m.Beam(m.knife_edge(m.gaussian(w), 0.5 * w, axis="x", side="+"))
beam = m.Beam(m.measured_profile(camera_image, pixel_size=5.5e-6))

Where the beam is defined matters: a top-hat imaged on the mirror and a top-hat focused at the centre are different inputs.

Builder	Beam defined
`Pulse.at_focus`	at the cell centre (matched Gaussian, tests)
`Pulse.from_near_field`	collimated, then focused into the cell by an achromatic lens (top-hat → Airy pattern)
`imaged_near_field_lens`	imaged on a plane of the cell, by default the first mirror, in physical coordinates there

pulse = m.Pulse.from_near_field(grid, E, m.Beam(m.top_hat(4e-3)), energy,
                                nf_size=12e-3, nf_points=192, target_waist=mm["w0"])
pulse = m.imaged_near_field_lens(grid, E, m.Beam(m.top_hat(1.12 * mm["w_mirror"])), energy, cell)

Top-hats, hard edges and Airy patterns need a larger box (box_factor 10–14).

Running

sim = m.LensFrameSimulation(pulse, cell, phi_max=0.02, gouy_step=0.05, mode_reference=mm["w0"],
                            profile_wavelengths=[1010e-9, 1030e-9, 1050e-9], profile_bandwidth=2e-9)
sim.run()                     # or sim.run(n_passes=1) to go pass by pass
sim.save("run.npz")

record_every, near_field=False (no collimated-beam diagnostics, a few seconds less per record), record_mirror: what is recorded.
profile_wavelengths: 1D profiles of single colours (bandpass profile_bandwidth) at the centre, in the collimated beam and on the mirror, at every pass.
Mirror records are in the frame of the beam: multiply positions by sim.mirror_scale and divide fluences by its square. The saved file keeps mirror_scale.

Long runs stop cleanly and resume:

if os.path.exists("ckpt.pkl"):
    sim.resume("ckpt.pkl")
sim.run(max_seconds=3600)     # stops after the pass that exceeds one hour
sim.checkpoint("ckpt.pkl")

Between passes, any element can be applied to sim.Aw (the field at the centre, physical in both solvers), for example an aperture:

aperture = (grid.r2 < (1.5e-3) ** 2)[:, :, None]
for p in range(cell.n_passes):
    sim.run(n_passes=1, verbose=False)
    sim.Aw *= aperture

New physics (ionisation, Raman, …) goes in the nonlinear term: subclass the solver and override nonlinear(), which works in \((x, y, \Omega) \leftrightarrow (x, y, t)\). LensFrameSolver inherits it; there \(|B|^2\) is the intensity times \(s^2\). Such a solver runs its RK4 in NumPy, not in the Numba kernels (Numba kernels (CPU)).

Results

sim.history has one entry per recorded pass: energy, B-integral, peak power and intensity, spectrum, fluence, spectrograms \(S(x,\lambda)\) and \(S(\theta,\lambda)\), Hermite–Gauss content, homogeneity, beam radii and \(M^2\) (Recorded history). The functions of thunder4d.diagnostics analyse any field:

import numpy as np
from thunder4d import diagnostics as dg

x, lam, S = dg.spectrogram_x_lambda(grid, sim.Aw)
c = dg.compress(grid, sim.Aw)                # best uniform GDD, compressed pulse
d = np.load("run.npz")                       # the saved history