Supplementary code example for the research article: “Complete characterization of ultrafast optical fields by phase-preserving nonlinear autocorrelation”
Published in: “Light: Science & Applications” (2022)
Authors: Alexander Gliserin (1, 2, 4), Soo Hoon Chew (1, 2, 4), Seungchul Kim (1, 2), Dong Eon Kim (3, 4)
Software author: Alexander Gliserin
E-mails: alex@mpk.or.kr (primary), novgorod@novgorod.de (secondary)
Affiliations: 
1) Department of Optics and Mechatronics Engineering, College of Nanoscience and Nanotechnology, Pusan National University, Busan, South Korea
2) Department of Cogno Mechatronics Engineering, College of Nanoscience and Nanotechnology, Pusan National University, Busan, South Korea 
3) Department of Physics, Center for Attosecond Science and Technology, Pohang University of Science and Technology, Pohang, South Korea
4) Max Planck Center for Attosecond Science, Max Planck POSTECH/Korea Research Initiative, Pohang, South Korea

=============================================================
This work is published under the Creative Commons Attribution-NonCommercial 4.0 International license (CC BY-NC 4.0).
See full license information at:
https://creativecommons.org/licenses/by-nc/4.0/
The software is provided as-is. Use at your own risk.
=============================================================

Scope:

This software package provides a reference implementation of the Phase-Enabled Nonlinear Gating with Unbalanced Intensity (PENGUIN) method for optical field retrieval from an unbalanced nonlinear interferometric autocorrelation (IAC). It allows simulating an IAC from a reference optical pulse field (synthetic or imported from external files) and retrieving the original field from the IAC using the generalized multi-order retrieval algorithm (see section S2 in the Supplementary Information of the original publication for the mathematical details).

-------------------------------------------------------------

Prerequisites:

The software is written in National Instruments LabVIEW, version 2018. The source code can be viewed, edited and executed using the National Instruments LabVIEW Full Development System version 2018 or later. The software may run on the Linux and MacOS versions of the LabVIEW Development System, but it was only tested on the Windows 10 operating system. A fully featured 7-45 days trial version of the LabVIEW Development System can be obtained from National Instruments:
https://www.ni.com/en-us/support/downloads/software-products/download.labview.html
A compiled executable (.exe) can be requested from the author, which only requires the free LabVIEW Runtime Engine to be installed.

-------------------------------------------------------------

Quick start:

1. Open the main program file “PENGUIN_demo.vi” and run it (Ctrl+R). By default, the reference spectrum file (reference_spectral_intensity.txt) is automatically loaded from the program directory together with a synthetic 4th-order polynomial spectral phase, reproducing the exemplary electric field used in the original publication. The harmonic order is n = 2 and the balance factor is s = 0.5 by default.

2. Switch to the “Field retrieval” tab and start the retrieval using the default parameters, either one step at a time (“Step” button) or automatically (“Auto” button). Unpress the “Auto” button after the retrieval has converged.

The default settings are suitable for n = 2 and n = 3 with various balance factors. Set “Update strength (q)” to 1 for faster convergence in the n = 2 case (this value is unstable for n = 3, hence the default value is s = 0.5). Since the IAC is invariant with respect to the pulse's linear spectral phase (CEP and delay), the retrieved pulse will slightly differ from the reference pulse in this regard. If the “Fit linear phase” checkbox is enabled, the linear spectral phase is fitted after each iteration to match the CEP and delay to the simulated reference field, resulting in minimized field intensity-envelope and phase errors. Enabling the linear phase fit has a significant performance penalty.

-------------------------------------------------------------

User guide:

- Data entry: All numerical controls with units use the SI prefix notation for convenience. E.g., you can enter either “100T” or “1E14” for a value of 100 THz, or “100a” or “1E-16” for a value of 100 as. Decimal digits of numbers can be incremented/decremented with the up/down arrow keys. Note that a value of 0 defaults to a unit without an SI prefix (e.g., 0 Hz), thus the increment/decrement uses the base unit in this case unless an SI prefix is manually entered.

- Spectral intensity source: Use the tab to select either a synthetic (super-)Gaussian spectrum or load an external spectrum from a text file. The synthetic spectrum is Gaussian for a supergauss order of 1 and becomes progressively flat-top (approaching a rectangular shape) for large supergauss orders. The expected file format for an external spectral intensity source is a tab-separated ASCII table with the first column being the frequency (in THz) and the second column the spectral intensity (arb. units). The first row is considered a header and omitted. The frequency sampling does not need to be equidistant, since the spectral intensity is interpolated according to the selected timebase parameters in the program. If the file path is empty on startup, a reference spectral intensity file in the program directory will be used (reference_spectral_intensity.txt). If the selected file cannot be loaded, the directory control will be colored in red.

- Spectral phase source: Use the tab to select either a synthetic polynomial spectral phase or load an external phase from a text file. The polynomial coefficients of the Taylor expansion are in units of fs^n where n is the polynomial order (starting with 0). The 0-th order describes the CEP (in rad), the 1st order the group delay (in fs), the second order the constant group delay dispersion or chirp (in fs²), etc. Arbitrary many coefficients can be used by adding elements to the coefficient array. The expected file format for an external spectral phase source is a tab-separated ASCII table with the first column being the frequency (in THz) and the second column the spectral phase (in rad). The first row is considered a header and omitted. The frequency sampling does not need to be equidistant, since the spectral phase is interpolated according to the selected timebase parameters in the program. If the file path is empty on startup, a reference spectral phase file in the program directory will be used (reference_spectral_phase.txt). If the selected file cannot be loaded, the directory control will be colored in red.

- Simulation parameters: The time step and span define the timebase and the corresponding frequency sampling interval and range. The nonlinear order (n) and balance (s) define the calculated IAC. Gaussian white noise can be added to the simulated IAC trace in the time domain using the noise parameter. The preview graphs in the “IAC simulation” tab are updated when any of the simulation parameters are changed, allowing for live parameter tuning before field retrieval.

- Initial field guess: By default, the initial field guess is obtained from the fundamental peak of the simulated IAC spectrum. The “Gaussian field guess” checkbox replaces this initial guess with a Gaussian spectral magnitude of the given full-width-at-half-maximum (FWHM) bandwidth and flat spectral phase.

- Filter settings (“Field retrieval” tab): The field retrieval employs spectral filters to separate IAC harmonics, as well as a time-domain filter to optionally confine the E-field in time. The default width of the temporal filter far exceeds the default simulation range, effectively disabling it. The filters are Gaussian with an adjustable supergaussian order (making them rectangular for large orders) and FWHM bandwidth. In addition, the spectral filter allows shifting its central frequency with respect to the numerically calculated central frequency of the spectral intensity (“IAC simulation” tab, bottom left). The fundamental spectral filter is used to isolate the fundamental peak of the IAC spectrum as the initial field guess and to constrain the bandwidth of the retrieved field after each iteration. It is also used to derive the harmonic filters to isolate other harmonics of the IAC for multi-order retrieval. The central frequency of the harmonic filters is linked to the fundamental filter by a factor of m (the number of the harmonic order), and their bandwidth is linked to the fundamental bandwidth by a factor of sqrt(m) for m > 0 and sqrt(2) for m = 0, respectively. Harmonic filters are only displayed in the “Spectra” graph if the corresponding harmonic weights are non-zero.

- Update strength (q): The q parameter controls the rate of convergence by mixing the updated field, E_<i+1>, with the previous field, E_<i>, after each iteration via qE_<i+1> + (1-q)E_<i>. Fastest convergence is achieved for q = 1, but the retrieval can become numerically unstable for n > 2, large noise or high dispersion, which requires lowering q.

- Linear spectral phase: The (unbalanced) IAC is invariant with respect to the original field's linear spectral phase (i.e., its CEP and global delay), thus the linear spectral phase of the retrieved field can differ from the reference field, resulting in an apparent delay or CEP mismatch. To assess the true performance of the field retrieval, the linear phase mismatch can be removed after each iteration by enabling the “Fit linear phase” checkbox. This sets the CEP of the retrieved field to the CEP of the reference field and also performs a numerical search for the global delay value which minimizes the field intensity-envelope error in the time domain. This fit is computationally expensive and substantially slows down the retrieval. For best performance, the linear phase fit can be done once on a single iteration after the retrieval has converged (minimized IAC error). The linear phase fit may fail if the delay mismatch is substantial (many cycles).

- Data export: The graph data on the “Field retrieval” tab can be exported using the corresponding “Save” buttons. “Save fields” exports the reference and retrieved fields as real parts, magnitudes and temporal phases with the linear phase (2pi*f_0*t, f_0 being the central frequency) subtracted as well as the temporal filter. “Save IAC” exports the real parts of the reference and retrieved IAC traces. “Save spectra” exports all spectral magnitudes and (unwrapped) phases shown in the “Spectra” graph as well as the currently enabled harmonic spectral filters. “Save error” exports the convergence graph with the time-domain IAC, intensity-envelope, and phase errors for each iteration step. All graphs are exported as tab-separated ASCII files with the first column being time/delay (in fs), frequency (in THz), or iteration step, respectively, and all other columns representing the graph data. The first row contains the column names.

- Multi-order retrieval: Numerically (near-)exact field retrieval is only achievable for low-dispersion pulses with all “Harmonic weights” set to zero except for the fundamental-order weight (see section S1 in the Supplementary Information of the original publication for a discussion of the single-order retrieval limit). Single-order retrieval yields a wrong solution for pulses above a dispersion limit, which depends on the particular spectral shape and dispersion type, and therefore requires using additional IAC harmonic peaks as a constraint to force the correct solution (see section S2 in the Supplementary Information of the original publication). This comes at the expense of amplifying numerical errors, since all harmonics of the IAC except for the fundamental retrieve some higher power of the electric field, which requires taking roots to obtain the fundamental field. The reference implementation of the multi-order retrieval method uses only the first non-trivial cross-correlation component for each harmonic number m (i.e., k = 1 for m = 0 and k = 0 for m > 0, following the formalism in the original publication) because higher-order components cannot improve the retrieval quality in practice and introduce excessive systematic error. A good retrieval quality with the multi-order method in the high-dispersion regime requires carefully adjusting the filter parameters (including the temporal filter in some cases), the harmonic weights and the update strength parameter (q). Depending on the actual pulse shape, the IAC parameters (n and s), and the amount and type of dispersion as well as noise, the retrieval may exhibit oscillatory or chaotic behavior and convergence cannot be achieved. In this case, the retrieval result is obtained at the iteration step which yields the lowest IAC error. If the “Hold best” checkbox is enabled, the displayed retrieval result is only updated if the IAC error could be reduced after the current iteration, otherwise the previous result with the lowest error is displayed (the E-field guess which is used internally for retrieval is still updated at every iteration). Since the IAC error behavior is non-monotonous in this case, the stopping criteria must be judged by the number of iterations since the last time the IAC error could be reduced.

Optimum retrieval parameters for the high-dispersion pulse used in section S2 in the Supplementary Information of the original publication, consisting of the reference spectral intensity and a purely quadratic spectral phase with a chirp of 40 fs² (no time filter was used):
For n = 2, s = 0.5: Filter width = 320 THz, Filter shift = 15 THz, Filter order = 40, DC weight = 200, Fund. weight = 1, 2nd order weight = 100, q = 1
For n = 3, s = 0.5: Filter width = 286 THz, Filter shift = 7 THz, Filter order = 10, DC weight = 500, Fund. weight = 1, q = 0.75