Covariance Analysis of Impulsive Streaking

We first derive the impulsive model for $K[X,Y]$ in the absence of machine fluctuations, i.e. when the only variation in the measurement is $\bm{k}$, in which case $K[X,Y]=C[X,Y]$. As illustrated in Eqn. (3), in the impulsive regime $X(\bm{r})\simeq X^{0}(\bm{r}-\bm{k}_{X})$, and $Y(\bm{r})\simeq Y^{0}(\bm{r}-\bm{k}_{Y})$. Substituting the Taylor expansions of $X$ and $Y$ into the definition of covariance Eqn. (5), we find the leading order is

where $\nabla$ is the 2D momentum gradient operator, $R(-\phi)$ is the 2-by-2 matrix representing the counterclockwise rotation of $-\phi$ in the $xy$-plane. This leading order is the inner-product between $\nabla{X^{0}}$ and $R(-\phi)\nabla{Y^{0}}$, also referred to as the gradient inner-product (GIP), and the delay is encoded in the rotation.

In the impulsive regime, the streaking covariance defined in Eqn. (6) is given by a bilinear operation on the unstreaked distributions, $C[{\mathcal{D}}_{\bm{k}_{X}},{\mathcal{D}}_{\bm{k}_{Y}}]{X^{0}}{Y^{0}}$. As shown in Eqn. (7), the leading order in this bilinear operator $C[{\mathcal{D}}_{\bm{k}_{X}},{\mathcal{D}}_{\bm{k}_{Y}}]$ can be written as

where $\hat{G}(\phi)\equiv\lambda\sum_{i,j}R_{ij}(-\phi)\partial_{i}\otimes\partial_{j}$, with the direct product defined as $\partial_{i}\otimes\partial_{j}\equiv\frac{\partial}{\partial r_{q,i}}\frac{\partial}{\partial r_{p,j}}$ to perform partial differentiation with respect to momentum coordinates $\bm{r}_{q}$ and $\bm{r}_{p}$, respectively. The full expansion of $C[{\mathcal{D}}_{\bm{k}_{X}},{\mathcal{D}}_{\bm{k}_{Y}}]$ is

where $\hat{H}\equiv\left(\nabla^{2}\otimes\hat{I}+\hat{I}\otimes\lambda^{2}\nabla^{2}\right)/2$, $\hat{I}$ is the identity operator, and $\gamma_{n}\equiv\langle k^{2n}\rangle/(2^{2n}(n!)^{2})$. According to Eqn. (9), the streaking covariance expands into multiple direct products between the $n_{X}$-th order partial derivative of $X$ and the $n_{Y}$-th order partial derivative of $Y$. We number them by the differentiation orders $(n_{X}+n_{Y})$, e.g. the GIP term is the (1+1) order. The next orders are (1+3), (3+1), and (2+2), since all terms’ $(n_{X}+n_{Y})$ are even numbers, as shown by Eqn. (9).

In a measurement, we record electron yield with finite momentum resolution. We consider the covariance between the electron yield measured in two arbitrary regions $Q$ and $P$. Due to the bilinearity of covariance, this is identical to the regional integral of the covariance between the densities $X$ and $Y$:

As illustrated in Fig. 2, we define two sets of 2D momentum regions of interest (ROIs): $\{Q_{q}\}_{q=1}^{N_{Q}}$ for $X$ and $\{P_{p}\}_{p=1}^{N_{P}}$ for $Y$. We use subscripts $q,p$ on $X,Y$ as a shorthand for the regional integrals, e.g. $X_{q}\equiv\int_{Q_{q}}X(\bm{r})d^{2}\bm{r}$, as illustrated in Fig. 2. Fig. 3(a) shows an exemplary visualization of $K[X,Y]$, where the covariance has been calculated between $N_{Q}{=}N_{P}{=}180$ ROIs on the $p_{z}{=}0$ plane, each $2^{\circ}$-wide. The 180 angular bins are on a ring $p_{\mathrm{min}}<|\bm{r}|<p_{\mathrm{max}}$ and labelled by their central angular positions $\theta_{q}$ or $\theta_{p}$, thus they are also referred to as angular bins. In Fig. 3, both photoelectron features $X$ and $Y$ are simulated in the same ionization process as Fig. 1 but detected separately, and we set $\lambda=1$ to simulate the scenarios where the difference in photoemission delays is much shorter than the dressing field period $\omega|\tau_{X}-\tau_{Y}|\ll 1$. As shown in the $\phi=\pi/3$ example in Fig. 3(a), one feature of $K[X,Y]$ is that the most positive part is around $\theta_{p}-\theta_{q}\sim\phi$. At small streaking amplitudes, $K[X,Y]$ is well approximated by the GIP, as shown in Fig. 3(b-c). The the relative error $\|M_{\mathrm{GIP}}-K[X,Y]\|_{2}/\|K[X,Y]\|_{2}$ between panel (b) and (a) is evaluated as $0.03$, which we use to quantify the accuracy of the GIP. The accuracy of the GIP model is insensitive to $\phi$, but it degrades as streaking amplitude $k$ increases, due to the increased contribution from the higher order terms. As shown by the solid curve in Fig. 3(c), the root-mean-squared (rms) relative error over $\phi\in[-\pi,+\pi]$ becomes comparable to unity once $k$ exceeds the width of the photoelectron feature in momentum $\Delta p$ (quantified as the full-width at half-maximum, FWHM). For panels (a-b), the lower boundary $p_{\mathrm{min}}$ of the ROIs is set to the maximum point of the unstreaked MD gradient $p_{\mathrm{MG}}\equiv\arg\max_{r}\int d\theta|\nabla I^{0}|$  (red solid line in Fig. 3(d)), and the upper boundary is set to where the streaked MD has fallen to zero. This choice of limit for $p_{\min}$ optimizes the accuracy of the GIP model, because $p_{\mathrm{MG}}$ is close to the zero point of $\nabla^{2}I^{0}$. As shown in Fig. 3(c-d), when we alternate $p_{\min}$ by a fraction of $\Delta p$, the relative error is increased when $k/\Delta p<1$.

For the photoelectron features produced in direct ionization, the discrepancy between the GIP model and the streaking covariance $K[X,Y]$ generally depends on three characteristic dimensionless quantities: (1) the size of the momentum shift normalized to the momentum spread of the photoelectron feature $k/\Delta p$, (2) the normalized radius $p_{c}/\Delta p$, with $p_{c}$ denoting the central momentum of the photoelectron feature, and (3) the ratio of the ionizing pulse duration to the streaking laser period $\Delta t_{X}/T_{L}$. In the impulsive regime $\Delta t_{X}\ll T_{L}$, we find that the relative error is predominantly a function of $k/\Delta p$, with a much weaker dependence on $p_{c}/\Delta p$, as shown in Fig. 3(f). In the case shown in panels (a-e), both $X$ and $Y$ have a normalized radius of $p_{c}/\Delta p=6.6$, and when we vary $p_{c}/\Delta p$ by changing $\Delta t_{X}$ and/or the photon energy of the ionizing pulse, the relative error curve does not significantly change as long as we remain in the impulsive regime ($\Delta t_{X}\ll T_{L}$). The two features $X$ and $Y$ can have different $\Delta p_{X},\Delta p_{Y}$, but we note that by defining $\Delta p$ as the geometric mean $\Delta p=\sqrt{\Delta p_{X}\Delta p_{Y}}$, the accuracy of the GIP can be well described by $k/\Delta p$.

We differentiate the shot-to-shot fluctuations in angular streaking experiments into two categories, the fluctuation of the streaking vector $\bm{k}$, and everything else. The fluctuations of these other parameters (or “machine fluctuations”), such as pulse energy and central photon energy of the ionizing pulse, also affect the single-shot MD, but the relative timing between the ionizing and dressing pulses typically does not depend on these parameters. Thus in this work, machine fluctuations are assumed to be statistically independent from the random streaking vector. The streaking covariance $K[X,Y]$ defined in Eqn. (6) is therefore simplified by the impulsive condition as $K[X,Y]\simeq C[\mathcal{D}_{\bm{k}_{X}}\langle X^{0}\rangle,\mathcal{D}_{\bm{k}_{Y}}\langle Y^{0}\rangle]$. We note that the expected dressing-free MDs $\langle X^{0}\rangle,\langle Y^{0}\rangle$ are fixed, free from shot-to-shot variations, so the treatment of $K[X,Y]$ introduced in Sec. II.1 are also applicable when machine fluctuations are present, by simply replacing the stable $X^{0},Y^{0}$ with $\langle X^{0}\rangle,\langle Y^{0}\rangle$ respectively.

Both machine fluctuations and fluctuations in $\bm{k}$ contribute to the covariance $C[X,Y]$ defined in Eqn. (5). According to the law of total covariance [24], $C[X,Y]=K[X,Y]+L[X,Y]$ consists of the streaking covariance $K[X,Y]$ defined in Eqn. (6) and the contribution from machine fluctuations:

where in the approximation, we used the impulsive streaking condition and the statistical independence between machine fluctuations and streaking vector.

To retrieve the delay, the GIP model of $K[X,Y]$ relies on the gradients $\nabla X_{q},\nabla Y_{p}$:

One way to reconstruct the gradients $\nabla X_{q},\nabla Y_{p}$ is through numerical differentiation (ND). Dressing-free shots provide the sample mean of the measured distributions $\overline{X^{0}}$ and $\overline{Y^{0}}$ to estimate the expectations $\langle X^{0}\rangle$ and $\langle Y^{0}\rangle$, respectively. With adequate momentum resolution, we can use a finite difference scheme to numerically differentiate $\overline{X^{0}}$, and then integrate over $Q_{q}$ to measure the gradient $\nabla X_{q}$. An example is shown in Fig. 5(a) for the case in Fig. 1 &3. Obtaining the gradients $\nabla X_{q},\nabla Y_{p}$ from numerical differentiation, we fit the GIP model to the sample streaking covariance $\mathrm{K}_{XY}$, with the factor $\langle k^{2}\rangle\lambda$ treated as a free parameter in Eqn. (12). Because the amplitude ratio $\lambda$ has been absorbed into the free scaling parameter, not knowing $\lambda$ does not affect this method. We benchmark the accuracy of the ND method in the noiseless limit: i.e. with sufficient number of shots $N_{s}\to\infty$ to suppress statistical noise, complete removal of machine fluctuations, and ignoring readout noise in the electron yield. Thus we equate $\mathrm{K}_{XY}$ to the underlying $K[X,Y]$ and obtain the delay retrieval error $\Delta\phi=\phi_{\mathrm{fit}}-\phi$ at various “ground-truth” $\phi$, Fig. 5(b). This systematic error is zero when $\phi$ is a multiple of $\pi/2$ and reaches a maximum at $\phi\,{\sim}\pm\pi/4$. The rms error over $\phi$ increases with $k$ due to the increase in higher-order terms, as shown in panel (c), but within $k<\Delta p$ this rms error is $<0.03\,\mathrm{rad}{=}1.7^{\circ}$, which converts to a time delay error $<30\,\mathrm{as}$ for a $2\,\mathrm{\mu m}$-wavelength dressing field.

Another way to reconstruct the gradient is “rank-reduction” (RR), described in Algorithm 1, which employs a low-rank approximation of $K[X,X]$ and $K[Y,Y]$ to estimate the gradient. This method can be used when each set of ROIs complete a loop on the high-momentum side of the corresponding photoelectron feature, e.g. the rings delineated in Fig. 4. Within the GIP model, $K[X_{q},X_{q^{\prime}}]=\langle k^{2}/2\rangle(\nabla X_{q})^{T}\nabla X_{q^{\prime}}$, whose rank is at most 2 since $\nabla X_{q}$ is 2D. Thus it appears we can solve the following optimization problem given $\mathrm{K}_{XX}$:

where $\xi=(\xi_{1},\cdots,\xi_{N_{Q}})\in\mathbb{R}_{2{\times}N_{Q}}$ represents the gradient field, $a_{q}$ is the arc length of the angular bin at $\theta_{q}$ and $W$ is a $N_{Q}\times N_{Q}$ matrix weights, which can be configured to combat the effect of readout noise as described below. The reconstructed gradient field is given by the optimal point $\xi^{*}$ of Eqn. (13) up to a global factor, $\xi_{q}^{*}=\sqrt{\langle k^{2}/2\rangle}\nabla X_{q}$. Applying the same procedure to $\mathrm{K}_{YY}$, we obtain the optimal point $\eta^{*}_{p}=\sqrt{\langle k^{2}/2\rangle}\lambda\nabla Y_{p}$. Thus the GIP model described in Eqn. (7) can be rewritten as $M_{\mathrm{GIP}}=(\xi^{*})^{T}R(-\phi)\eta^{*}$, which we then substitute into Eqn. (12) to retrieve $\phi_{\mathrm{fit}}$.

The minimization of the objective function $f_{\mathrm{RR}}$ is closely related to principal component analysis (PCA). The weight matrix can be configured as uniform $W_{qq^{\prime}}=1$ when the signal-to-noise ratio is high in the electron yield readout process. When the readout signal-to-noise ratio is low, we recommend ignoring the diagonal band of $\mathrm{K}_{XX}$ by setting the diagonal band in $W$ to zero (e.g. $W_{qq^{\prime}}=1-\delta_{qq^{\prime}}$). Guidance for defining is detailed in Sec. IV.2. With a uniform $W$, a minimal point of $f_{\mathrm{RR}}$ is given by the PCA result $\xi^{\mathrm{P}}=\left(\sqrt{s^{(1)}}v^{(1)},\sqrt{s^{(2)}}v^{(2)}\right)^{T}$ [25], where $s^{(\alpha)}$ is the $\alpha$-th largest eigenvalue of $\mathrm{K}_{XX}$ and $v^{(\alpha)}$ is the corresponding eigenvector (in column vector form). When $W$ is non-uniform, e.g. $W_{qq^{\prime}}=1-\delta_{qq^{\prime}}$, we can minimize $f_{\mathrm{RR}}$ by optimizing $\xi^{\mathrm{P}}$ using gradient descent.

The objective function $f_{\mathrm{RR}}$ is invariant under any global orthogonal transformation $O{:}~{}\xi_{q}\mapsto O\xi_{q}$, but only certain minimal points can satisfy the zero-loop constraint Eqn. (13b), which is a general property of a gradient field $\oint\nabla X^{0}\cdot d\bm{l}=0$. Thus obtaining one minimal point $\xi^{\mathrm{P}}$, we solve for an orthogonal matrix $O$ such that $l_{\mathrm{RR}}(O\xi^{\mathrm{P}})=0$ This is straightforwardly achieved by parameterizing $O$ with the rotation angle and parity. An alternative to the zero-loop constraint is maximizing the gradient-flux $j_{\mathrm{RR}}(O\xi^{\mathrm{P}})\equiv-\sum_{q}(\cos\theta_{q},~{}\sin\theta_{q})O\xi^{\mathrm{P}}_{q}a_{q}$ with $O$, which necessarily satisfies the zero-loop constraint and additionally breaks the remaining discrete degeneracy, as proven in Supplemental Sec. 2.

Similar to the results in Fig. 5(b-c), we benchmark the RR method under various $\phi$ and $k$ in the noiseless limit as mentioned above, as shown in panels (d-f). The gradient fields $\nabla X_{q}$ and $\nabla Y_{p}$ are reconstructed from $\mathrm{K}_{XX}$ and $\mathrm{K}_{YY}$, respectively, and we show the resultant $\nabla X_{q}$ in panel (d). Similar to the ND method, the delay retrieval error of the RR method is minimized when $\phi$ is a multiple of $\pi/2$. The rms error increases with $k$, but the magnitude of the error is notably smaller than with the ND method. As shown in Fig. 5(f), the rms error is ${<}1$ mrad within $k{<}\Delta p$, corresponding to ${<}1$ as for a $2\,\mathrm{\mu m}$ dressing field.

The main reason for the higher accuracy of the RR method lies in the reconstructed “gradient”. The ND gradients give the first-order derivatives of the dressing-free MD $\nabla X_{q}$ and $\nabla Y_{p}$, independent from the streaking amplitude $k$. The gradient reconstructed by the RR procedure, in contrast, deviates from the first-order derivative as $k$ increases. To the next lowest order, the RR gradients are

which include the third-order derivatives in addition to the first-order (see Supplementary Sec. 1.2). As a result, when using the RR gradients, the inner product $\xi_{q}^{T}R(\phi)\eta_{p}$ not only captures the (1+1) order of $K[X_{q},Y_{p}]$, but also the (1+3) and (3+1) orders. In contrast when using the ND gradient, only the (1+1) order is modelled. Since the errors in the noiseless limit arise from the finite accuracy of the GIP model and are insensitive to the normalized radius $p_{c}/\Delta p$, the result in Fig. 5 is general across $p_{c}/\Delta p$.

In some experiments, gradient reconstruction may be challenging. For example, the numerical differentiation becomes infeasible when the detector lacks angular resolution, and RR becomes infeasible when the available ROIs cannot complete a loop. In this case, it is often possible to identify a part of the streaking correlation matrix $R[X_{q},Y_{p}]\equiv{K[X_{q},Y_{p}]}/{\sqrt{K[X_{q},X_{q}]K[Y_{p},Y_{p}]}}$ that can be used to extract the photoemission delay without requiring knowledge of the gradient.

In the angular regions where the gradient is predominantly along the radial direction $|r^{-1}\partial_{\theta}X_{q}|\ll|\partial_{r}X_{q}|$, $|r^{-1}\partial_{\theta}Y_{p}|\ll|\partial_{r}Y_{p}|$ (also referred to as the “radial regions”), the streaking covariance is approximated as

The $\cos(\phi-\theta_{p}+\theta_{q})$ factor explains the positive ridge around $\theta_{p}-\theta_{q}{\sim}\phi$ in the covariance matrix. The radial gradients in Eqn. (15), including the factor $\lambda$, cancel out, i.e. $R\left[X_{q},Y_{p}\right]\approx\cos(\phi-\theta_{p}+\theta_{q})$. Therefore we define the gradient-free model as $M_{\mathrm{GF}}\equiv a\cos(\phi-\theta_{p}+\theta_{q})$, with free parameters $a$ and $\phi$. From the measured K-estimators, we calculate the sample correlation matrix $\mathrm{Corr}_{XY}=\mathrm{K}_{XY}/\sqrt{\mathrm{K}_{XX}\mathrm{K}_{YY}}$, and minimize the mean-squared error in the radial region between $\mathrm{Corr}_{XY}$ and the model $M_{\mathrm{GF}}$.

The radial region is determined from the dressing-free MD. For the case shown in Figs. 3 and 5, the anisotropy parameter is $\beta_{2}=2$, the radial region defined by $|r^{-1}\partial_{\theta}X_{q}|<0.2|\partial_{r}X_{q}|$ is $138^{\circ}$-wide around each antinode of the dipole feature. For lower anisotropy parameters, the radial region is larger, since an isotropic feature has no angular gradient component. As shown in Fig. 6(a), within the radial region, the correlation matrix is well described by the gradient-free model $M_{\mathrm{GF}}$. By fitting $M_{\mathrm{GF}}$ to the measured $\mathrm{Corr}_{XY}$ in the radial region, we obtain $\phi_{\mathrm{fit}}$, and the error $\phi_{\mathrm{fit}}-\phi$ is shown in Fig. 6(b & c). Similar to the gradient-based methods, the error is zero at $\phi=0,\pm\pi/2$, but the rms error across $\phi\in[0,2\pi]$ does not vanish when $k\to 0$. This residual error at small streaking amplitude arises from ignoring the non-radial regions, and in the case shown in Fig. 6, it amounts to $20$ mrad rms. The systematic error of the gradient-free method also increases with $k$, but the increase is slower than for the ND method. In this case, comparing Fig. 6(c) to Fig. 5(c), we find the gradient-free method is more accurate than the ND method when $k\gtrsim\Delta p$, but remains less accurate than the RR method. Another difference between the gradient-free and gradient-based methods is the dependence on $p_{c}/\Delta p$: the gradient-free method is generally more accurate with higher $p_{c}/\Delta p$, since the magnitude of the radial gradient relative to the angular gradient increases.

When implementing analysis procedures that make use of fluctuations in the measured data, it is important to understand how instabilities in the measurement affect the measured correlation, in this case $\mathrm{C}_{XY}$. When using an attosecond x-ray free electron laser there may be fluctuations in the x-ray parameters, for example the pulse energy and/or the central photon energy. In this section we describe how to manage the effect of these additional fluctuations in the delay extraction procedure.