Quantitative uncertainty determination of phase retrieval in RABBITT: supplement

Reconstruction of attosecond beating by interference of two-photon transitions (RABBITT) is one of the most widely used approaches to measure the time delays in photoionization. The time delay, which corresponds to a phase difference of two oscillating signals, is usually retrieved by cosine fitting or fast Fourier transform (FFT). We propose two estimators for the phase uncertainty of cosine fitting from the signal per se of an individual experiment: (i) σ(φ fit)≈B A2N, where B/A is the mean-value-to-amplitude ratio, and N is the total count number, and (ii) σ(φ fit)≈1-R 2 R 2 n bins, where nbins is the total number of bins in the time domain, and R2 is the coefficient of determination. The former estimator is applicable for the statistical fluctuation, while the latter includes the effects from various uncertainty sources, which is mathematically proven and numerically validated. This leads to an efficient and reliable approach to determining quantitative uncertainties in RABBITT experiments and evaluating the observed discrepancy among individual measurements, as demonstrated on the basis of experimental data.


Quantitative uncertainty determination of phase retrieval in RABBITT: supplemental document 1. GENERAL FORMULA FOR LINER FITTING AND THE DISTRIBUTION OF THE CO-EFFICIENT OF DETERMINATION
A vector of F = [F(t 1 ) . . .F(t n bins )] T is fitted to F fit (t) = α cos 2π T t + β sin 2π T t + B by: . . . ) F fit (t 1 ) . . .
The fitted signal solely depends on the column space of M, denoted as C(M), but the basis can be arbitrary.By Schmidt orthogonalization, we can find a 2-dimensional linear space V = span{v 1 , v 2 } that satisfies span{1, v 1 , v 2 } = C(M) and 1 ⊥ V, where 1 is the first column of M. Therefore, F = F • 1 + v + ε ε ε = F fit + ε ε ε, where v ∈ V and ε ε ε is the fitting residual perpendicular to C(M).Strictly speaking, F = B only when there are integer number of periods; nonetheless, ε ε ε is always perpendicular to v. Hence, A perfect cosine signal gives ε ε ε = 0 and R 2 = 1.On the other hand, the expected value of R 2 for a white noise is 2/n bins , instead of 0, because the white noise is uniformly distributed in the frequency domain.In particular, R 2 n bins ∼ χ 2 (k = 2), since the positive-and negativefrequency components are uncorrelated.Filtering the signal by FFT does not reduce the effect from background noise, since only the contamination at the frequency of interest is relevant to the phase fitting, which is unaffected by filtering.
If both the 2ω-and the 4ω-components are included, R 2 n bins ∼ χ 2 (k = 4), and for white noise we have:

NUMERICAL SIMULATIONS
The simulated RABBITT signals are generated as following.Let the probability density function be: Two random variables t and u are drawn from uniform distributions in [0, n T T] and [0, 1], respectively.If u ≤ f (t)/(A 2ω + A 4ω + B), then this t is accepted as an observed event at time t.The generation process is repeated until N events are accepted.These events are binned as a histogram, as shown in Fig. S1.The distribution of fitted α, β, and ϕ are compared in Fig. S2 (no external background noise or shifting is added).The difference between the sampling B/A (of PDF) and the fitted B/A is caused by the binning effect, and the numerical results highly agree with our prediction.Figure S3 demonstrates that the expected R 2 of the simulated RABBITT signal agrees with the theoretical value, and the distribution is narrower for smaller B/A ratio and greater N, indicating that the uncertainty estimated from the signal with low uncertainty level is more precise on single-experiment basis.The effect of correlated jitter for the ideal cosine signal and the simulated RABBITT signal are illustrated in Fig. S4 and Fig. S5, respectively.Both the expected values and the variances of fitted α, β, and A are well-described by the lowest-order expansion for σ s /T ≤ 0.05, which corresponds to 67 as.The effect of the 4ω-component is shown in Fig. S6, where the R 2 -estimator that includes both the 2ω-and the 4ω-component is more precise at the low-uncertainty regime.However, this leads to a lower plateau at the high-uncertainty regime, as predicted by theory.Nonetheless, the estimation is "safe" when the individual uncertainty does not exceed ∼ 0.3 rad, which corresponds to ∼ 60 as.
The period-fitting uncertainty for the simulated RABBITT signals is plotted in Fig. S7.The numerical results indicate that σ(π∆T/T), i.e. the phase change caused by the deviation of the fitted period, is proportional to N −1/2 B/A.The effect on R 2 is shown in Fig. S8.For cosine fitting, this effect is negligible, since the fitting uncertainty is dominated by the target channel with several orders of magnitude fewer counts than the reference channel.For FFT, however, the spectral leakage will result in significant reduction of R 2 , which is determined by the non-integer

EXAMPLE WITH REAL EXPERIMENTAL DATA A. Attosecond Beamline
The attosecond beamline is based on High Harmonic Generation (HHG) producing an XUV attosecond pulse train (APT), which is co-linearly overlapped with a 35 fs short IR pulse in a Mach-Zehnder like geometry.The laser system is a Coherent Legend Elite Duo HE+, producing a 10 W IR beam with 800 nm central wavelength, 5 kHz repetition rate and 35 fs short pulses.This IR beam is split into two pathways with a 70/30 beamsplitter.The more intense part of the beam is focused (f=45 cm) into a 6 mm long gas cell filled with 10 mbar of xenon.The generated HHG pass through a 100 nm thick aluminum filter to filter out the remaining IR and are then focused into a gas jet in the coincident spectrometer by a nickel-coated toroidal mirror (f=50 cm).The less intense part of the beam is co-linearly overlapped with the APT by a holey mirror and focused via a lens (f=1.1 m), so that the focal positions of both XUV APT and IR are centered in the gas jet.After the coincident spectrometer the XUV APT enters an XUV spectrometer, based on a reflective grating (1200 lines/mm), an MCP (micro-channel-plate), a phosphor screen and a CCD camera to characterize the high harmonics.The pathlength of the IR arm can be controlled with attosecond precision using a piezo-electric motor and is actively stabilized in an interferometric approach using a frequency-stabilized helium-neon laser coupled into the XUV-IR Mach-Zehnder interferometer [1], resulting in a time jitter ranging from 40 to 100 as.
To detect the charged fragments we use a COLTRIMS (COLd Target Recoil Ion Momentum Spectroscopy) spectrometer [2,3], where the electrons and ions are guided by weak homogeneous electric fields of 2.6 V/cm and 2.8 V/cm, respectively, onto two opposite time-of-flight and position sensitive detectors.The detectors each consist of two 75-mm diameter micro-channel plates (Photonis) in Chevron configuration and a three-layer delay-line anode (HEX) with an active radius of 40 mm and a crossing angle of 60 degrees between the adjacent layers (RoentDek).A pair of Helmholtz coils is used to apply and additional homogeneous magnetic field (6.4 G) and cancel the influence of the earths magnetic field.On the electron side the region with the weak electric field present (extraction region) is 70-mm long, followed by a 140-mm electric-field-free region.On the ion side there is only a 400 mm long extraction region.Using COLTRIMS we gain access to the full three-dimensional momentum vectors of electrons and ions in 4π solid angle.The count rate was set to a maximum ion count of 0.3 per laser shot, to ensure coincidence detection.The argon gas is supersonically expanded into the vacuum chamber through a 50-µm diameter nozzle with a backing pressure of 0.5 bar.Two conical skimmers with orifices of 200 µm and 1 mm are located 10 mm and 30 mm respectively.The gas jet crosses the IR and XUV beams before it is captured by a differentially pumped beam dump.The description of the setup combining COLTRIMS with APT can be found in ref [1].

B. Angular-resolved time delays
The time delays at SB-12, 14, and 16 with weighting by the B/Aand the R 2 -estimator are compared in Tab.S1 ∼ S5.The comparison of individual measurements are shown in Fig. S9.The uncertainty of the weighted mean value calculated from individual signals is consistent with the statistical approach.

Fig. S2 . 2 N
Fig. S2.Simulated RABBITT signal with N = 10 5 , B/A = 4, n T = 5, n p = 10, ϕ = 1.0.(a) 2-D distribution of fitted α and β.The correction takes the binning effect into account .(b) The distribution and standard deviation (std) of fitted ϕ; the theoretical std is B A 2 N = 0.0179.(c) (e) show the distributions and variances of fitted α, β, and A, respectively; the theoretical variances are all equal to 2B/n bins = 80.

Fig. S3 .
Fig. S3.Numerical test for the R 2 method.Simulated RABBITT signals in one oscillation period are generated with default: B/A = 4.0, n p = 100, and N = 10 4 .Each simulation scans one parameter.The mean values and error bars of R 2 correspond to 100 independent trials for each green circle.(a) Scan over B/A ratio.(b) Scan over n bins = n p ; (c) Scan over N. The expected values are shown as the black dashed lines.

Fig. S4 .
Fig. S4.Correlated jitter of an ideal signal with n bins = 500.The relative standard deviations of α fit , β fit , and A fit correspond to 5000 independent trials at each point.(a), (b), and (c) show the relatively small shifts, while (d), (e), and (f) show the relatively large shifts.The black dashed lines indicate the theoretical values under the lowest order expansion.

Fig. S5 .
Fig. S5.Fitted parameters of a simulated RABBITT signal with N = 10 4 , B/A = 4.0, n T = 5, n p = 10, ϕ = 0 under different correlated jitter levels.The mean values and error bars correspond to 5000 independent trials at each point.The black dashed lines indicate the theoretical expected values with correction of bin migration effect.

Fig. S6 .
Fig. S6.Correlated jitter of a simulated RABBITT signal with probability density function proportional to [4 + cos(2ωt) + 0.3 cos(4ωt)] and n bins = 50.Each red circle corresponds to 5000 independent trials.The line shadow indicate the mean value and standard deviation of the estimated σ(ϕ fit ) for each trial, respectively.(a) and (b) use the R 2 that includes both the 2ωand the 4ω-component, while (c) and (d) use the R 2 that uses only the 2ω-component.

Fig. S7 .
Fig. S7.Simulated RABBITT signals are generated with default: N = 10 4 , B/A = 4.0, n T = 5, n p = 10.Both the phase and oscillation period are fitted with non-linear least-squares fitting.Each simulation scans one parameter.The standard deviation is calculated based on 100 independent trials for each magenta circle.(a) Scan over N (in logarithm scale).(b) Scan over B/A ratio.(c) Scan over n T (with interval of 5/3).(d) Scan over n p (with interval of 17/19).

Fig. S8 .
Fig. S8.The discrepancy between the discrete frequencies of FFT and non-integer oscillation of the real signal results in ∆ϕ of the fitted phase and reduces the R 2 value, as shown by the red circles.The clue curve corresponds to the lowest-order expansion: R 2 = 1 − 2−3/π 2 6 (∆ϕ) 2 .

Table S1 . Angular-resolved Time Delays for Ar at SB-12 from 18 Individual Experiments with (w) and Without (r) Weighting. The Weighting Is Based on the B/A-method θ
range (°) ∆t r (as) ∆t w (as) s r ∆t (as) σ w ∆t (as) s w ∆t (as) M eff S

Table S2 . Angular-resolved Time Delays for Ar at SB-14 from 18 Individual Experiments with (w) and Without (r) Weighting. The Weighting Is Based on the R 2 -method θ
range (°) ∆t r (as) ∆t w (as) s r ∆t (as) σ w ∆t (as) s w ∆t (as) M eff S

Table S3 . Angular-resolved Time Delays for Ar at SB-14 from 18 Individual Experiments with
(w) and Without (r) Weighting.The Weighting Is Based on the B/A-method θ range (°) ∆t r (as) ∆t w (as) s r ∆t (as) σ w ∆t (as) s w ∆t (as) M eff S

Table S4 . Angular-resolved Time Delays for Ar at SB-16 from 18 Individual Experiments with (w) and Without (r) Weighting. The Weighting Is Based on the R 2 -method θ
range (°) ∆t r (as) ∆t w (as) s r ∆t (as) σ w ∆t (as) s w ∆t (as) M eff S

Table S5 . Angular-resolved Time Delays for Ar at SB-16 from 18 Individual Experiments with (w) and Without (r) Weighting. The Weighting Is Based on the B/A-method θ
range (°) ∆t r (as) ∆t w (as) s r ∆t (as) σ w ∆t (as) s w ∆t (as) M eff S