# Machine Learning Enhances Algorithms for Quantifying Non-Equilibrium
Dynamics in Correlation Spectroscopy Experiments to Reach Frame-Rate-Limited
Time Resolution

Tatiana Konstantinova  Lutz Wiegart  Maksim Rakitin  Anthony M DeGennaro
<EMAIL_ADDRESS>Andi M Barbour<EMAIL_ADDRESS>Brookhaven National
Laboratory
Upton, NY 11973, USA

###### Abstract

Analysis of X-ray Photon Correlation Spectroscopy (XPCS) data for non-
equilibrium dynamics often requires manual binning of age regions of an
intensity-intensity correlation function. This leads to a loss of temporal
resolution and accumulation of systematic error for the parameters quantifying
the dynamics, especially in cases with considerable noise. Moreover, the
experiments with high data collection rates create the need for automated
online analysis, where manual binning is not possible. Here, we integrate a
denoising autoencoder model into algorithms for analysis of non-equilibrium
two-time intensity-intensity correlation functions. The model can be applied
to an input of an arbitrary size. Noise reduction allows to extract the
parameters that characterize the sample’s dynamics with temporal resolution
limited only by frame rates. Not only does it improve the quantitative usage
of the data, but it also creates the potential for automating the analytical
workflow. Various approaches for uncertainty quantification and extension of
the model for anomalies detection are discussed.

††preprint: APS/123-QED

##  I Introduction

Technological advances, such as high-brightness X-ray photon sources [1, 2, 3]
and high-rate high-sensitivity detectors [4, 5], enable new discoveries by
means of X-ray Photon Correlation Spectroscopy (XPCS) [6]. Non-equilibrium
dynamics are signature of many systems studied with XPCS, including
superconductors [7], metallic glasses [8], proteins [9] and polymeric
nanocomposites [10]. Such dynamics are commonly represented via two-time
intensity-intensity correlation functions ($2TCF$s) [11], defined by the
expression:

$C2(\pmb{q};t_{1},t_{2})=\frac{\langle
I(\pmb{q};t_{1})I(\pmb{q};t_{2})\rangle}{\langle
I(\pmb{q};t_{1})\rangle\langle I(\pmb{q};t_{2})\rangle}$ (1)

where $I(\pmb{q};t)$ is the intensity of a detector pixel corresponding to the
wave vector $\pmb{q}$ at time $t$. The average is taken over pixels with
equivalent $\pmb{q}$ values. The quantitative analysis of a non-equilibrium
$2TCF$ (Fig. 1(A)) often starts with its binning along the sample’s age axis
_$t_{a}$_ into quasi-equilibrium regions [12, 13] prior to making cuts along
the time delay axis _$t_{d}$_ and averaging them, i.e. obtaining ‘aged’ one-
time intensity-intensity correlation functions ($1TCF$s)[14] of delay time:

$C1(\pmb{q},t_{a};t_{d})={\langle C2(\pmb{q};t_{a},t_{d})\rangle}_{t_{a}}$ (2)

The reason for binning and averaging is increasing the accuracy of further
analysis. The resulting $1TCF$s are fit to a functional form that is
characteristic for the dynamics under investigation.

Figure 1: Examples of intensity-intensity correlation functions. (A) A non-
equilibrium $2TCF$. $t_{1}$ and $t_{2}$ are experiment time measured in
frames. Yellow arrows indicate the directions of age axis _$t_{a}$_ and time
delay axis _$t_{d}$_ , defined according to the formulas in the inset. Pink
dashed lines and the grey solid line represent the edges and the bin center,
respectively, for the region used for generating the $1TCF$ in (B). (B)
$1TCF$s for a single cut along the bin center (grey solid) and average within
the bin of 50 frames(pink dashed).

The traditional analysis of experimental data with non-equilibrium dynamics
presents several challenges. Selection of the quasi-equilibrium regions is
often a non-trivial and time-consuming task, which requires interventions from
a human researcher for a visual inspection of the calculated $2TCF$ and
intermediate results. Moreover, the temporal resolution of the evolution of
dynamics’ parameters deteriorates due the binning procedure. Besides, with
increasing data collection rates, it often becomes unfeasible to properly do
on-the-fly data analysis during XPCS experiments when such analysis requires
inputs from a researcher. This can lead to inefficient use of beam facilities
and the efforts of scientific stuff as well as to unfortunate omission of
important parameter space regions during measurements. The sheer volume of
collected data can reach hundreds of thousands of $2TCF$s a day, making even
the offline data processing a challenging task. Automating the analytical
routine is an inevitable step towards high-throughput autonomous XPCS
experiments. Ability to perform the analysis with experimental frame-rate
resolution is the ultimate target.

Since the motivation behind binning and averaging the data is to curtail the
impact of noise, pre-processing a $2TCF$ with a noise-reducing algorithm would
increase the chance of achieving a frame-rate resolution along _$t_{a}$_
during quantitative analysis. Previously, we demonstrated an approach for
noise reduction in $2TCF$s with equilibrium dynamics, which is based on
Convolutional Neutral Network–Encoder-Decoder model (CNN–ED) [15]. The method
is shown to work well for equilibrium data and has a promising potential for
use with slowly ageing dynamics when applied as a sliding window along
_$t_{a}$_. However, there are two major concerns regarding the use of the
model for general cases of non-equilibrium dynamics. Firstly, if the rate of
the dynamics changes fast, an equilibrium approximation may not be valid, even
in narrow age regions. Thus, the application of the model would lead to loss
of experimental resolution and extraction of less accurate results. Secondly,
for a $2TCF$ containing a much larger number of frames than the input size of
the model, a considerable portion of it would still contain the original level
of noise after the model is applied along _$t_{a}$_.

Here, we build upon the previous model to address the reduction of noise in
arbitrary-sized $2TCF$s with non-equilibrium dynamics such as ageing. We
further demonstrate how the new model, integrated into a workflow for
quantitative analysis of XPCS data, eliminates the need of age binning and
allows wider parameter bounds for the fit of $1TCF$s. The results of such
analysis have temporal resolution, which is limited only by the experimental
acquisition rates. The methods for estimating the creditability of the results
– uncertainty scoring for denoised $2TCF$s and _trust regions_ for dynamics’
parameters – are discussed. Various options for analysis workflows are
demonstrated for several XPCS experiments.

##  II Methods

###  II.1 Data

For a denoising model to perform well for non-equilibrium $2TCF$s, data from
experiments with these types of dynamics are used for training. The data are
collected from 65 XPCS experiments at the CHX beamline of NSLS-II from three
different samples, recorded at different conditions and with different
detector collection rates. Some of the measurements are repeated multiple
times at the same conditions. The samples exhibit dynamics, common for many
materials, that at various rates monotonically accelerate with age,
monotonically decelerate with age or stay at quasi-equilibrium. For all
considered dynamics, individual $1TCF$s can be approximated with the
Kohlrausch-Williams-Watts (KWW) form:

$C1(t_{d})=C_{\infty}+\beta e^{-2(\Gamma t_{d})^{\alpha}}$ (3)

where $\Gamma$ is the rate of the dynamics, $\beta$ is the contrast factor,
$\alpha$ is the compression constant and $C_{\infty}$ is the baseline. All
parameters are functions of $t_{a}$ and $\pmb{q}$. Several regions of interest
in the reciprocal space are considered for each experiment when calculating
$2TCF$s. There is a total of 492 $2TCF$s, ranging from 134 to 2950 frames. The
original full-sized $2TCF$s were split between training (454) and validation
(38) sets prior to generating the inputs to the model and augmenting the data
as described below.

Inputs are obtained from raw experimental $2TCF$s by cropping the
100$\times$100 frames pieces along the $t_{a}$ in a similar way as was done
previously [15]. Prior to cropping, the diagonal values, containing high
errors, are replaced by the average of their nearest off-diagonal neighbors.
That is, the value at [j, j] is replaced by the mean of the values at [j-1, j]
and [j+1, j]. This approach preserves the information in the experiments with
slowly varying contrast. For data augmentation, additional inputs are
constructed by considering every 2nd, 3rd, 4th and 5th frame of the original
data when enough frames are available to form at least one 100$\times$100
example. The final size of the training and validation sets are 25222 and 1557
cropped examples respectively. Each input is scaled to have zero mean and unit
variance prior to passage through the model. A reverse transformation is
performed for the model output.

###  II.2 Denoising Autoecoder Model

The CNN-ED model architecture, used for equilibrium data, demonstrated several
advantages, such as simplicity, control of overfitting as well as fast
training and application. Trying the same model architecture with certain
adjustments for non-equilibrium data is a natural choice. The model presented
here consists of an encoder with two 10-channel convolutional layers,
8-dimensional latent vector and the decoder with two 10-channel transpose-
convolutional layers. Two modifications are implemented during the model
training to meet the peculiarities of non-equilibrium dynamics. Firstly, since
a ‘noise-free’ version of an input cannot be obtained by averaging multiple
cropped inputs collected from the same full-sized $2TCF$, the model is trained
in an autoencoder (AE) mode with an input and its target being the same.
Secondly, a non-equilibrium $2TCF$ cannot be described by a single $1TCF$ and
thus only the mean square loss between the output and the target are chosen as
the training loss function.

Despite the increased complexity of dynamics in the training data, the
architecture still appears appropriate for the denoising task based on the
model’s good performance for the validation set [see Appendix D, E]. The
effectiveness of models with convolutional kernel sizes from 1 to 17 is tested
and only minor improvements in the smoothness of the output and the level of
finite details preserved for larger kernels are discovered. The models with
larger kernels take longer to train and the computational time during their
application also becomes significant. Further in the text, unless stated
otherwise, a model with the kernel size 1 is considered. The fact that the
same model architecture works well for both equilibrium and non-equilibrium
ageing $2TCF$s confirms it robustness for different types of dynamics.

###  II.3 Adapting Model to Arbitrary Input Size

Since the model has a fully connected layer (the latent vector), it can only
be applied to a fixed-size input, i.e. 100$\times$100 frames. However, the
dimensions of an experimental $2TCF$ can vary considerably. The symmetry of a
$2TCF$ prevents the model to be applied as a raster scan over the entire
image. Application of the model in a sliding window fashion along the age axis
$t_{a}$ can help with the noise reduction in certain cases like in out
original quasi-equilibrium analysis [15]. Nevertheless, this approach is not
suitable for general cases of $2TCF$s, where dynamics may happen outside of
the first 100 delay frames along the $t_{d}$ axis.

Figure 2 shows a possible solution, which is inspired by the procedure of data
augmentation when the original $2TCF_{raw}$ is mapped down (_down-mapped_) to
a smaller size by considering every N-th frame horizontally and vertically,
where N = 2, 3, 4, 5. For model application, N is a whole part of division of
total number of frames in $2TCF_{raw}$ by 100. This way, the size of a newly
constructed $2TCF_{reduced}$ matches the required model’s input size. By its
definition (Eq. 1), a $2TCF$ is symmetric with respect to $t_{a}$. However,
when the starting point for _down-mapping_ is away from the $t_{a}$ axis
($t_{d}>$0) of $2TCF_{raw}$, the reduced version is not symmetrical. To ensure
the symmetry of the input, only the bottom part of the $2TCF_{reduced}$ is
considered and a reflection across the $t_{d}$=0 line is applied to it. For
all starting points, the $2TCF_{reduced}$ values at $t_{d}$=0 are extrapolated
by the neighboring points in the same way it is done for model’s training.
After application of the model, its output $2TCF_{out}$ is mapped back (_up-
mapped_) to the size of the $2TCF_{raw}$ by placing [i, j]-th point from the
$2TCF_{out}$ to [N*i, N*j]-th point in the resulting denoised
$2TCF_{denoised}$. Again, only the bottom half is filled and then reflected
across the $t_{d}$=0 line. The process is repeated for all possible starting
points, covering fully the $2TCF_{raw}$.

Figure 2:  Scheme of the model application to an input of arbitrary size.

For regions of very fast dynamics and for large-sized $2TCF$s, this _down-and-
up mapping_ approach may impair the experimental temporal resolution. In such
situations, the values near the diagonal of $2TCF_{raw}$ can instead be
processed in a sliding-window fashion, when the model is successively applied
to parts consisting of cropped 100$\times$100 frames equally distanced along
ta, and the outputs are averaged. The results of the model’s application in
the _down-and-up mapping_ and sliding window fashion are combined to obtain
the final denoised version of the original $2TCF_{raw}$. This method reduces
noise in the entire $2TCF_{raw}$ and preserves the experimental temporal
resolution for both $t_{a}$ and $t_{d}$.

###  II.4 Uncertainty Quantification for Denoising

When applying a model to unseen data, there are three potential sources of
uncertainty: error (noise) of the input, model’s bias and model’s variance.
The bias of the model primarily originates from selection of the training
dataset and constricting the model’s architecture. Variance appears due to
random initializations of model’s weights and ordering of data butches during
model’s training. Approaches for quantifying uncertainty for a deep learning
model often involve random perturbations of the model during training and/or
application and obtaining the distribution of respective outputs [16].

We estimate the variance of the denoising AE by considering deviations of
outputs for models trained with different random initializations (see Appendix
E). It appears that the model’s variance is generally much smaller than the
corresponding values of $\beta$. Moreover, the fluctuations of values in the
neighboring points of a $2TCF_{denoised}$, caused by _down-and-up mapping_
during model application, scale linearly with the inherent variance of the
model. Thus, it is not necessary to separately estimate the model’s variance
for experiments with more than 100 frames, as the point-to-point variations in
a $2TCF_{denoised}$ already influence the quantitative analysis of its
$1TCF$s.

From the practical point of view, the main potential source of uncertainty of
the model is its bias. Obviously, a denoising AE must exhibit _some_ bias to
remove noise. Hence, to drive the decision about applicability of the model in
each case, the quantification of the bias ought to answer the question: how
certain can one be that the model’s output is a valid representation of the
underlying sample’s dynamics in the input? Naturally, a model returns the most
trustworthy results for the inputs that are very similar to the examples in
the training set. We suggest two quantitative measures – uncertainty scores –
for estimating the bias by comparing a new input to the training examples.

The first measure is based on the Euclidean distance between the latent vector
representation of an input and the center of the latent representations of all
training examples. The analysis of the pairwise distributions of the 8 latent
coordinates of the training set reveals that they form a single compact
cluster [see Appendix F]. The shape of the cluster is determined by the
distribution of dynamics and the noise level in the training set. However, the
values of each of the eight latent coordinates are mostly concentrated around
a single point. The distance from this central point is thus indicative of how
different a given input is from most of the training data and thus can be used
as an uncertainty score. To ensure an equal contribution of each latent
dimension to the uncertainty score, all coordinates are scaled to have a zero
mean and unit variance across the training dataset. To improve the
interpretability of the bias, the distance to the center is normalized by its
median value among the training set examples. The normalization of the
distance helps to set a general threshold for an acceptable level of the bias
without reference to a particular model or the training set.

The second measure of bias of the denoising is less abstract and does not
depend on the model architecture. It involves evaluation of trends in
residuals, i.e., the differences between the model’s targets and corresponding
outputs. When the output is a valid representation of the sample’s dynamics,
the residual should resemble random and sometimes correlated noise in form of
vertical and horizontal stripes, without trends along neither $t_{a}$ nor
$t_{d}$. Since the training examples mostly represent monotonically ageing
dynamics, the model tends to perform well for similar cases. A trend in a
residual can indicate a heterogeneity on top of the ageing dynamics or a
completely different type of dynamics. Such situations require additional
attention during analysis. To identify the trends, it is convenient to look at
the projection of a two-dimensional residual to the vertical (or horizontal)
direction. Without trends, the autocorrelation coefficients (ACCs) of this
projection are close to zero. In contrast, when a trend is present, the
absolute value of the ACCs grows. We calculate the first ACC for all the
examples of the training set and calculate probability density distribution
using Gaussian kernel density estimator [see Appendix G]. The probability
density for the first ACC of the residual for a new test example is the second
type of uncertainty scoring of the model prediction.

Since both bias measures are based on the distribution of the training
examples, they are suitable for detecting anomalies – inputs that are
different from the training set. The analysis workflow can flag a new $2TCF$
as an anomaly if its uncertainty score is above the user-defined threshold.

###  II.5 Extracting Dynamics Parameters

Quantitative analysis of XPCS experimental results with non-equilibrium
dynamics, such as ageing, involves taking time-slices along the delay axis
$t_{d}$ of a $2TCF$ at different sample’s ages $t_{a}$ and fitting the
resulting $1TCF$s with an appropriate model that describes the dynamics, which
in many cases is the KWW form Eq. 3. Typically, several adjacent slices within
age range $\Delta t_{a}$ need to be averaged to target a signal-to-noise ratio
that is sufficient for extracting dynamics parameters with reasonable
certainty. It is necessary that the dynamics do not change much within $\Delta
t_{a}$. An example of selecting a quasi-equilibrium bin within a $2TCF$ is
shown in Fig. 1. This approach leads to a loss of temporal resolution for the
parameters. Moreover, the procedure of binning $1TCF$s often requires an
repetitive evaluation of the fitting results and adjustment of bins’
boundaries by an expert researcher. For a large volume of collected data, it
becomes challenging to properly split $2TCF$s for each experimental region of
interest (ROI), leading to growing uncertainty of the results. Ability to
perform an adequate quantitative analysis for a single value of $\Delta t_{a}$
across the entire experiment would facilitate the automation of the analysis
process. The highest temporal resolution and maximum usage of experimental
data are achieved when considering cuts with $\Delta t_{a}$ = 1
$\frac{1}{frame\ rate}$ (or, $fr^{-1}$).

We compare the possibility of conducting a quantitative analysis of a
$2TCF_{raw}$ and a corresponding $2TCF_{denoised}$ while considering $1TCF$
cuts with $\Delta t_{a}$ = 1 $fr^{-1}$ without a prior knowledge of the
dynamic’s parameter ranges (Fig. 3). An example of a noisy $2TCF_{raw}$ (350
frames) from the validation set (Fig. 3A) is passed through the model,
resulting in the $2TCF_{denoised}$ with a significantly reduced level of noise
(Fig. 3B-C). The latent space representations of the 100$\times$100 frames
$2TCF_{r}educed$ are close to the center of the training set (Fig. 3D),
indicating that the model output is likely a valid representation of the
dynamics in the input. To extract the parameters, each of the $1TCF$ cuts is
fitted to Eq. 3. Cuts with less than 5 points are not considered. The results
of the fit are shown in Fig. 3(E-H). The fits for the $1TCF$ cuts from the
$2TCF_{raw}$ with $\Delta t_{a}$ = 35 $fr^{-1}$ are provided as ground truth
values. Apparently, high noise in the $2TCF_{raw}$ does not allow extracting
meaningful information about temporal evolution of dynamics parameters without
restricting the fit parameters within narrow regions or increasing the bins’
width. The $2TCF_{denoised}$, on the other hand, produces clear, slowly
evolving trends, matching the ground truth values.

Figure 3:  Comparison of quantitative analysis for $2TCF_{raw}$ and
$2TCF_{denoised}$. (A)$2TCF_{raw}$, (B) $2TCF_{denoised}$, (C) $1TCF$s taken
along the line in (A-B), (D) distribution of latent vector lengths for the
training set (green) and the $2TCF_{reduced}$s (orange), (E-J) plots for
dynamics parameters values according to the legend. The _trust regions_ are
calculated for fits of $2TCF_{denoised}$.

While binning and averaging $1TCF$s reduces the noise, it is not a universal
solution because material’s dynamics can vary considerably within an age bin,
and therefore, a single set of parameters cannot describe them. For the
example shown in Fig. 3, this becomes apparent when considering averaging all
available $1TCF$s and fitting the result to Eq. 3. Since $\beta$ and
$C_{\infty}$ are not changing in this experiment, the resulting values for
these parameters are close to the ground truth values. However, the values for
$\Gamma$ and $\alpha$ are not close to their corresponding average values.

For an automatic analysis, it is important to flag the results that cannot be
fully trusted, even when the goodness-of-fit to Eq. 3 is acceptable. If
dynamics are not fully captured by the experiment, parameters of the Eq. 3
become mutually dependent and the outputs of non-linear regression can be
misleading because several different sets of the dynamics’ parameters can
result in almost equally good fits. Therefore, we introduce the concept of
_trust regions_ for the parameters. A _trust region_ of a parameter is a
binary vector with the length equal to the number of bins along $t_{a}$, which
indicates whether the parameter is likely to be reliably identified within
each bin. The binary values are determined by several criteria including the
rate of the dynamics, correlations between parameters and parameters’ relative
errors.

When a material’s dynamics are slow with respect to the duration of the
experiment, the baseline $C_{\infty}$ cannot be reliably identified. Likewise,
the true value of $\beta$ cannot be extracted if the dynamics are not fully
captured by the experiment due to them being either too fast or too slow. The
quantitative measure for establishing respective thresholds for dynamics’ rate
is a half-time, i.e., the time it takes for the contrast to drop by half:

$T_{1/2}=\left(\frac{ln2}{2}\right)^{\frac{1}{\alpha}}\frac{1}{\Gamma}$ (4)

For slow dynamics, when the half-time is larger than a user-defined portion of
$1TCF$’s length, the _trust region_ values for $C_{\infty}$ and $\beta$ are
set to zero: not trustable. In case of fast dynamics, when the half-time is
less than a certain number of frames, the _trust region_ for $\beta$ is set to
zero. However, the $C_{\infty}$ can be reliably extracted for fast dynamics
and hence its _trust region_ values are set to one. Similarly, threshold-based
conditions are defined for the correlation coefficient between parameters,
relative errors of the parameters and $R^{2}$ measure of the fit.

##  III Results

###  III.1 Analysis Workflow

Depending on the application, there are various ways to design an analysis
workflow that includes $2TCF_{raw}$s and $2TCF_{denoised}$s as well as any
prior knowledge about the material’s dynamics. An example flowchart for
extracting dynamical parameters from a $2TCF_{raw}$ is shown in Fig. 4. In
this workflow, a $2TCF_{denoised}$ is used for flagging unusual observations
and narrowing parameter boundaries based on the initial fit. The new
boundaries are then used to fit the $2TCF_{raw}$ and the $2TCF_{denoised}$.
The last step of the analysis results in two sets of parameters with
corresponding errors and _trust regions_ : one for the $2TCF_{raw}$ and one
for the $2TCF_{denoised}$. Step-by-step dynamics extraction algorithms are
provided in Appendix H.

Figure 4:  Analysis workflow. Shape keys: rhombus – data supplied to the
algorithm, rounded shape - model, diamond shape – result of calculations.
Color keys: green – results of XPCS measurements, purple – supplied by a user,
orange – intermediate results, red – the final results.

Naturally, when reporting experimental outcomes, the fitting results of
unprocessed $2TCF_{raw}$ have a higher priority with respect to the fitting
results of $2TCF_{denoised}$ because signal processing, such as noise removal,
introduces an uncertainty. According to the presented analysis workflow,
$2TCF_{denoised}$ can be used for evaluating the parameters’ boundaries for
fitting the raw data during the second pass with the same or greater bin
width. However, it is also possible to supplement the fitting results for
$2TCF_{raw}$ with fitting results for $2TCF_{denoised}$ for regions in
parameter space where the low signal-to-noise ratio in raw data prevents from
obtaining meaningful parameter values. Note that removing the noise is just an
additional step in extracting parameters from the raw experimental data.
Therefore, the fitting results for $2TCF_{denoised}$ can be tested against
$2TCF_{raw}$ with goodness-of-fit measures such as $R^{2}$. Alternatively, one
can use the fitting results for $2TCF_{denoised}$ to select quasi-equilibrium
regions for the analysis of $2TCF_{raw}$.

Upon calculating the uncertainty score for the $2TCF_{denoised}$, one may wish
to extend the workflow to investigate possible dynamic heterogeneities of the
sample by relying on the denoised output and the results of the $1TCF$s’ fit.
The bias of the denoising model provides the opportunity to separate the
average ‘envelope’ dynamics and stochastic heterogeneities by subtracting
either the $2TCF_{denoised}$ or the $2TCF$, computed based on the fit
parameters for the $2TCF_{denoised}$, from the original $2TCF_{raw}$.

###  III.2 Application to New Data: Standard Analysis

We test a workflow for XPCS analysis with the use of denoising AE on
previously reported data for 3D printing with solution of Lithium Titanate
particles [17]. The experiment is not a part of the training or the validation
datasets used in the current work. However, the far-from-equilibrium dynamics
of the ink, exhibited during its deposition and recovery, are similar to the
types of dynamics used during the models’ training. Thus, the experiment
presents an intended use case for the denoising AE.

Figure 5:  Comparison of use $2TCF_{raw}$ and $2TCF_{denoised}$ from 3D
printing experiments [17] for quantitative analysis. (A) $2TCF_{denoised}$ and
bins of width 1 frame (dashed lines). (B) $1TCF$s corresponding to the bins in
(A) (open circles) and their fits to Eq. 3 (solid lines). (C) $2TCF_{raw}$ and
bins of width 1 frame (dashed lines). (D) $1TCF$s corresponding to the bins in
(C) (open circles) and fits to denoised data from (B) (solid lines). (E)
$2TCF_{raw}$ and bins of width 100 frames (dashed lines). (F) $1TCF$s
corresponding to the bins in (E) (open circles) and fits to denoised data from
(B) (solid lines). (G)drift velocity obtained from a set of experimental
scattering series, calculated using $2TCF_{denoised}$ (solid teal circles) and
from Ref. [17]. Results from $2TCF_{denoised}$ have 7854 points and the
previous results have 19 points.

There are several factors that limit the signal-to-noise ratio in this
experiment: operando character of the measurements, beam-sensitivity of the
ink and anisotropy of the dynamics, which requires selection of small ROIs in
the reciprocal space for analysis. In the original work, unevenly spaced age
bins of various width have been selected to obtain good quality ‘aged’ $1TCF$
slices [17]. Here, we test the advantages of applying the denoising AE for
conducting the same analysis with bin width $\Delta t_{a}$ = 1 $fr^{-1}$.

Figure 5 provides an example of the denoising AE being applied to the
$2TCF_{raw}$ for one of the ROIs. Since the dynamics is relatively fast, the
points near the $t_{a}$ axis are denoised when the model is applied in a
sliding window fashion and the points away from $t_{a}$ are denoised via
_down-and-up mapping_ approach. The difference in level of noise for
$2TCF_{raw}$ and $2TCF_{denoised}$ can be seen from corresponding $1TCF$s
taken at different ages. Bin width of 1 frame has been considered for
$2TCF_{denoised}$ and resulting $1TCF_{denoised}(t_{a})$ are fitted to Eq. 3
to obtain $1TCF_{fit}(t_{a})$. $1TCF_{raw}(t_{a})$ for bin widths 1 and 100
are compared against the same $1TCF_{fit}(t_{a})$. As expected, the fit to
denoised version of correlation function describes the raw data well. As the
width of the bins increases, the values of $1TCF_{raw}(t_{a})$ lie closer to
the corresponding $1TCF_{fit}(t_{a})$ lines. The model captures well not only
the characteristic time of the signal decorrelation, but also the change in
the contrast factor $\beta$, caused by initial fast non-ergodic dynamics with
timescales outside of the experimental time window.

The sample’s dynamics for multiple ROIs is quantified according to the scheme
in Fig. 4. We calculate the drift velocity from relaxation times at different
wave vectors along deposition direction $\Phi=-90^{\circ}$ using only denoised
$2TCF$s and compare them to the values obtained in Ref. [17], as shown in Fig.
5G. The values obtained from $2TCF_{denoised}$ are very close to results of
the original analysis, confirming the denoising AE does not distort the data.
Instead, it eliminates the need to carefully select the age bins and ensures
the ultimate temporal resolution of a single acquisition period.

###  III.3 Application to New Data: Outlier Detection

By using latent space representation of a denoising AE, it is possible to
quantitatively compare $2TCF$s from the same group of measurements even if the
model was not exposed to this type of dynamics during training. This is
suitable for anomaly detection in series of consecutive XPCS measurements for
the same material. We test the application of the model for instabilities
detection during data collection for static dynamics of a LSAT sample (MTI
Corporation) at the CSX beamline of NSLS-II. Six data series, 7200 frames
(exposure 1.0875 seconds) each are purposefully collected during synchrotron
accelerator testing to introduce a disturbance of experimental conditions.
Three of the series are collected in a stable beam regime, one series is
recorded with small disturbance of the beam energy and two series are recorded
during strong beam energy disturbance. From each series, 72 model’s inputs
(size 100$\times$100 frames) are generated by considering every 72nd frame
with different starting points. Each of the _down-mapped_ inputs from a single
series contain the same information, but different noise. To make the
comparison sensitive to the mean value and variance of the contrast, the
inputs are not scaled, but are clipped within values of 1 and 2. The values at
$t_{d}$ = 0 are extrapolated with the neighboring values as it is done for a
routine model application.

Figure 6: Detection of anomalies in static experiments by encoding data with
the denoising AE model. (A) points in the latent space corresponding to inputs
generated from the experimental series. The points are mapped to 2D space
using linear principal components. B,D,F – $2TCF$s for stable(B), slightly
unstable(D) and very unstable(F) experimental conditions. C,E,F – average
intensities on the detector for the respective experiments on the left.

Figure 6 shows the representations of stable and unstable series in the latent
space. Inputs formed from the same full-sized $2TCF$ are located in close
proximity of each other, forming a tight cluster. This confirms the
contractive property of the denoising AE - the fact that similar inputs are
close in the latent space. Having multiple $2TCF$s representing the same
dynamics through _down-mapping_ procedure allows estimating the characteristic
size of the cluster, to which all other points in the latent space should be
compared to. The means and variances of contrasts for the inputs collected
from experiment with a small disturbance are similar to the undisturbed
series, resulting in a small distance in the latent space, comparable to the
size of the clusters. However, when the disturbance is so strong, that the
beam is lost for part of the series, the latent space representations of a
$2TCF$ are much further away from the undisturbed series. In fact, the
strongly disturbed series form a separate cluster, where points are much more
spread out than for stable series. Selecting a threshold distance between
points helps to filter out the measurements taken at very unstable
experimental conditions. Thus, it is possible to identify unusual observations
in a set of XPCS series by looking into the denoising AE’s latent
representation of $2TCF$s and comparing the distance between series.

##  IV Discussion

Analysis of XCPS data is a multistep, often iterative, process that requires
continuous evaluation of intermediate results by a domain expert. Large
volumes of observations and subjectivity in alleviating low signal-to-noise
ratio complicate proper extraction of valuable information from experiments.
Computationally reducing noise in $2TCF$s as a data processing step helps to
achieve several goals: automating the analysis workflow, improving temporal
resolution of parameters that quantify the system’s dynamics and increasing
quantitative usage of data with high cost of collection.

Here, we demonstrate how a denoising AE model can be included into analysis
workflow of XPCS data for experiments with ageing dynamics. Quantification of
the model’s bias helps driving the decision about the use of its outputs as
well as flagging unusual observations, such as heterogeneities or dynamics
that are very different from ageing. The denoised correlation function can be
used for complementing the fit of the raw data. The concept of _trust regions_
combines the assessment of the fit quality as well as domain expertise, which
helps to not only report the most reliable results, but to automate sequential
narrowing of the fit parameters boundaries.

The model’s performance for unseen 3D printing data (not included in the
training/validation datasets) demonstrates the advantages of its application
to actual XPCS experiments. For measurements where the material’s dynamics are
similar enough to the model’s training dataset, the analysis does not require
a human-in-the-loop after all the thresholds are selected prior to analysis.
Moreover, the analysis is suitable for autonomous data acquisition. It
provides values of dynamics parameters, which are important for making
decisions about adapting data acquisition parameters such as acquisition rate,
exposure time, duration of data acquisition or about changing the sample
and/or the processing parameters.

We further demonstrate how encoded representations of a $2TCF$ can be used for
quantitative comparison of two or more scattering time series, which can be
useful for identifying anomalies such experimental instabilities or phase
transitions. The comparison can even be done for the types of dynamics, not
present in the training set.

In conclusion, in this work we demonstrate how a CNN-based denoising AE can be
used for $2TCF$s with non-equilibrium dynamics for experimental scattering
series of arbitrary size. Addition of the denoising AE to XPCS analysis along
with quantifying uncertainty helps automating the analysis and improves
temporal resolution of extracted parameters. Besides, analysis of residuals
and latent space representations of the inputs helps detect anomalous dynamics
that go beyond monotonic ageing. This property can be employed for recognizing
heterogeneities or phase transitions. Several examples of incorporating the
denoising AE in analysis workflows for actual experimental data, collected at
NSLS-II, demonstrate its effectiveness for unseen data and diversity of its
applications. Denoising and encoding properties of the model are auspicious
for various online and offline data analysis tasks at XPCS facilities.

###### Acknowledgements.

This research used the CHX and CSX beamlines and resources of the National
Synchrotron Light Source II, a U.S. Department of Energy (DOE) Office of
Science User Facility operated for the DOE Office of Science by Brookhaven
National Laboratory(BNL) under Contract No. DE-SC0012704 and under a BNL
Laboratory Directed Research and Development (LDRD) project 20-038 ”Machine
Learning for Real-Time Data Fidelity, Healing, and Analysis for Coherent X-ray
Synchrotron Data”

##  Appendix A Denoising Model Training

The denoising model is written in Pytorch (version 1.9.0+cuda10.2). The model
is trained using CUDA accelerator Nvidia GeForce RTX 2070 Super. Early
stopping is used for model overfitting. Aside a few changes, the same model
architecture and training procedure is used as described in [15].

The model is trained in the autoencoder regime, i.e. each input serves as its
target. The cost function is the mean square loss between the model’s output
and its target(input). Adam [18] optimizer is used for training. Its starting
learning rate of 1e-4 (batch size 8) is decreased at each epoch by a factor
0.9995.

One training epoch takes from 28 seconds (kernel size 1) to 320 seconds
(kernel size 17). 20-30 epochs are required for training, depending on random
initialization of the model’s weights and the order of training examples
supplied to the model.

##  Appendix B Fitting a $2TCF$

Each individual $1TCF$ is fit to Eq. 3 using nonlinear regression implemented
in Scipy package [19] according to the Trust Region Reflective algorithm [20].
Generally, when a $1TCF$ is calculated for a single-frame-wide cut of a
$2TCF$, each of its points is given equal weight during the fit and when
several-frames-wide cuts are considered, the weights are assigned inversely
proportionally to the standard errors of correlation function values at each
$t_{d}$ point.

The optimization algorithm for nonlinear regression is sensitive to the
initial guess of the parameters. The initial values for $\beta$ and
$C_{\infty}$ can be estimated from characteristics of the experimental setup
and test measurements. The algorithm seems to be robust with respect to the
initial value of $\alpha$, which is set to 1. However, a proper initial guess
of $\Gamma$ is important and $\Gamma$ values can vary significantly between
experiments and ROIs in the reciprocal space. To automate selection of the
initial guess for $\Gamma$, several values are attempted and the fit with the
smallest $R^{2}$ is returned. The values range between 0.01 and 1 _fr_. If
neither of the values results in better fit than a line $1TCF(t_{d})=const$,
then the error (_‘not a number’_) is returned for all parameters.

It is common to have a situation when a signal decorrelates within a small
delay and most of the points along the $t_{d}$ axis in a $1TCF$ are at the
baseline level. In such case, points at the tail of the $1TCF$ have high
leverage, but carry little information regarding the dynamics’ rate, the
contrast factor and the compression constant. The equidistant delay points in
$1TCF$s obtained through cuts of a $2TCF$ contrast with logarithmic-spaced
delay points of $1TCF$s, calculated directly from the speckled images using
the traditional multi-tau algorithm [21]. To reduce the influence of the
points at the tail of the $1TCF$, lower weights can be assigned to them during
the fit with nonlinear regression. For implementing this, at the first pass,
the data are fit to Eq. 3 with points weighted either equally or according to
standard errors, as described above. Based on the results of the first fit,
the points for which the fitted values are larger than $1+\frac{\beta}{e}$ are
assigned the weight of 1 during the second fit and the rest of the points are
assigned the weight $\frac{1}{log(t_{d})}$. The same data are fit for the
second time to Eq. 3 with the fixed value $C_{\infty}$ from the previous fit
and the new points’ weights. Both weighted and unweighted results are returned
by the fitting procedure.

The nonlinear regression algorithm returns the parameter values as well as
their covariance matrix. The values of $\Gamma$ and $\alpha$ for the best fit
are used for calculating the half-time (Eq. 4) and establishing initial _trust
regions_. Square roots of variances (the diagonal elements of the covariance
matrix) are attributed to the parameters’ errors. The off-diagonal elements
are converted to correlations. Optionally, the parameters’ errors and
correlations and the $R^{2}$ values are used for narrowing the _trust regions_
based on user-defined thresholds.

##  Appendix C Output Files Formats

$2TCF_{denoised}$ and uncertainty scores are saved as datasets in a _hdf5_
[22] file, which allows accessing data (loading into memory) in parts.
Parameters, parameter’s errors, parameter’s correlations and the $R^{2}$ of
the fits at each age are recorded as nested dictionaries into a _json_ [23]
file. Besides, the parameter boundaries and the initial guesses are also saved
in a _json_ file.

##  Appendix D Different Kernel Sizes

Comparison of the denoising AEs with different convolutional kernels shows
that models with larger kernel sizes have less bias and closer resemble the
dynamics even if they slightly deviate from the KWW form. Examples of model
outputs for the validation data are shown in Fig. 7 and Fig. 8. The figures
show the model output for the $2TCF$s with 300 frames, rather than for single
100-frame inputs (see the main text for details of applying the model to an
arbitrary-sized input). Generally, with increasing the kernel size of the
model, the trends in the residuals become less pronounced.

Figure 7:  Comparison of denoising with models with different kernel sizes
(1,3,11,17) for relatively slow dynamics. (A) $2TCF_{raw}$, (B-E)
$2TCF_{denoised}$, (F-J) residuals, (I-M) $1TCF$s for cut at
$t_{a}=150frames$, (N-Q) $1TCF$ for $t_{a}=50frames$. $\Delta t_{a}=1frame$.
Color scale is the same within each row. Figure 8:  Comparison of denoising
with models with different kernel sizes (1,3,11,17) for relatively fast
dynamics. (A) $2TCF_{raw}$, (B-E) $2TCF_{denoised}$, (F-J) residuals, (I-M)
$1TCF$s for cut at $t_{a}=150frames$, (N-Q) $1TCF$ for $t_{a}=50frames$.
$\Delta t_{a}=1frame$. Color scale is the same within each row.

However, the difference between the models’ outputs is generally not
significant for quantitative analysis, which involves further fitting of the
$1TCF$s. We compare the error of parameters extracted from the $2TCF$s passed
through the models with different kernels for all full-sized $2TCF$s from the
validation set. For this, we fit single-frame-wide cuts of a $2TCF_{raw}$ and
its denoised version $2TCF_{denoised}$ in an iterative process according to
Fig. 4. The initial fitting of $2TCF_{denoised}$ is used for establishing
parameter boundaries for the contrast factor $\beta$ and the baseline
$C_{\infty}$, applied during the second fit. The ground truth values are
obtained by separately fitting 20-frames-wide cuts of $2TCF_{raw}$ and
interpolating the results of the fit across the entire range of ages. The
ground truth fits involve human intervention for establishing the parameters’
boundaries.

One $2TCF$ is removed from consideration because its dynamics are too fast to
be captured with the _down-and-up mapping_ approach. Two other $2TCF$s are
removed from consideration because the _trust regions_ for $\beta$ from the
first fit of $2TCF_{denoised}$ contained less than two points and thus the
parameter boundaries could not be established. Those $2TCF$s contained very
slow dynamics. Overall, 8100 single-frame $1TCF$s have been fit for each
kernel value.

The distribution of errors for the models with kernels 1, 3 and 17 is shown in
Fig. 9. For $\Gamma$ the relative error (RE$\Gamma$) is considered since
dynamics’ rate can be in a broad range, including values close to zero. The
absolute error is considered for $\alpha$ (AE$\alpha$), $\beta$ (AE$\beta$)
and $C_{\infty}$ (AE$C_{\infty}$). Table 1 shows the comparison of median
values of the errors. The analysis procedures involving all models produce
comparable errors.

Figure 9:  Distribution of errors for parameters extracted from $2TCF$s in the validation set, denoised with models having kernel size 1 (circle), 3(triangle) and 17(square). Relative error for dynamics rate $\Gamma$ (A); absolute error for stretching constant $\alpha$ (B), $\beta$ (C) and $C_{\infty}$ (D). kernel |  RE$\Gamma$(model) |  AE$\beta$(model) |  AE$C_{\infty}$(model) |  AE$\alpha$(model)
---|---|---|---|---
1 |  0.133 |  0.0076 |  0.0023 |  0.21
3 |  0.135 |  0.0081 |  0.0026 |  0.23
7 |  0.127 |  0.0075 |  0.0023 |  0.21
11 |  0.121 |  0.0069 |  0.0023 |  0.21
17 |  0.132 |  0.0075 |  0.0025 |  0.22
kernel |  RE$\Gamma$(raw) |  AE$\beta$(raw) |  AE$C_{\infty}$(raw) |  AE$\alpha$(raw)
NA |  0.195 |  0.0088 |  0.0025 |  0.35
Table 1:  Errors for parameters extracted from a $2TCF$ denoised with
different kernels and from the raw data for examples from the validation set.

The errors produced by fitting $2TCF_{raw}$ have larger median errors for both
RE$\Gamma$ and AE$\alpha$, even when the parameter bounds for $\beta$ and
$C_{\infty}$ are narrowed. The difference in errors for $\beta$ and
$C_{\infty}$ between fits for $2TCF_{raw}$ and $2TCF_{denoised}$ is not
significant because in both cases the parameter boundaries were established
using $2TCF_{denoised}$ at the previous step.

When considering denoising AEs with different kernel sizes, one needs to
account for the computational time. The comparison of computational times for
denoising models with different kernel sizes is shown in Fig. 10. These times
are measured for a single 100-frame $2TCF$ input (average of 50 repeated
measurements). Depending on the acceleration hardware and the need to transfer
data to/from a CPU/GPU Pytorch Tensor, the difference between computation
times for models with kernel=1 and kernel=17 can be 60 times (CPU) or can be
non-existent (GPU, without transfer).

Figure 10:  Computation times for application of the denoising AE model
(single 100x100 frames input) with/without GPU acceleration and with/without
data transfer from CPU-based numpy arrays to CUDA-based Pytorch tensors. Error
bars are standard deviations for measurements repeated 50 times.

##  Appendix E Model Variance

We estimate the variance of a denoising AE model by applying an ensemble of 9
models trained with different random initializations to examples from the
validation set. For each full-sized $2TCF$ from the validation set, a reduced
version $2TCF_{reduced}$ is generated by considering each 3rd frame
horizontally and vertically, starting at point [0,0], and then cropping the
first 100x100 frames part. The resulting $2TCF_{reduced}$ is passed to each of
the models in the ensemble and the frame-wise variance ($\sigma^{2}$) is
calculated. That is, each of the 100$\times$100 frames input produces $10^{4}$
variances. We refer to these variances as generated via drawing from _model
space_.

When a single denoising AE model is applied to a $2TCF$ with 200 or more
frames, the _down-and-up mapping_ process also introduces variance between
neighboring values. The difference in underlying dynamics represented in the
_down-mapped_ $2TCF$ originating from starting points [i, i] (i, j = 0,1,2) is
often negligible. Therefore, the _down-mapped_ $2TCF$s can be viewed as an
approximation of realization of the same dynamics. Thus, it is possible to
estimate the variance of a model via drawing from _input space_ by applying
the model to each of such realizations and estimating the variances of the
outputs. Here, we estimate the variance by considering inputs, obtained from
the $2TCF$s from the validation set by considering each 3rd frame horizontally
and vertically, starting at points: [0,0], [1,1], [2,2] and then cropping the
first 100$\times$100 frames part from each of the _down-mapped_ $2TCF$. A
single model is applied to each of the 3 inputs and the frame-wise variance is
calculated.

As a result, for each example in the validation set the frame-wise variance is
calculated by drawing from both _model space_ and _input space_. The
variances, calculated by both methods, are comparable. (Fig. 11) The variance
values are considerably smaller than both the noise in the original raw inputs
and the typical values of contrast factor $\beta$, making the model’s variance
negligible for the quantitative analysis.

Figure 11:  Variance of the model. (A,C) A $1TCF$ with $\Delta t_{a}=1fr^{-1}$
obtained from a $2TCF_{raw}$ (green) of two different experiments and from
respective average of $2TCF_{denoised}$s for ensembles drawn from the _model
space_ (orange) and the _input space_ (blue). (B) Standard deviations of
$2TCF_{denoised}$ values for ensembles drawn from _model space_ versus
respective values for ensembles drawn from the _input space_. Each point
corresponds to a value of a $2TCF$ at a single point $[t_{1},t_{2}]$. Red line
is a linear fit between the variables. (D) Distribution of the standard
deviation values from (B).

Moreover, there is a linear dependence between the variances calculated by the
two methods (the correlation coefficient is 0.73), indicating that the model
variance typically is already reflected in the _down-and-up mapping_ process
of the model application. Therefore, it is not necessary to separately
estimate the variance of the model using an ensemble method.

##  Appendix F Latent Coordinates Distribution

Latent space coordinates can be used for applying distance-based similarity
measures for the model inputs. The model is likely to perform well for new
inputs that are similar to the training set examples. The pair-wise
distributions of the latent coordinates for the training set (Fig. 12) reveals
that the majority of the data is approximately centered around a single point.
As a result, the measure of similarity between a new input and the training
set can be expressed as the Euclidean distance from the input to the central
point of the training set coordinates distribution. Prior to calculation of
the distance, the coordinates are standardized using the mean values and
variances for each coordinate calculated for the training set.

Figure 12: Pairwise distribution of latent coordinates for the training set.
The coordinates are scaled to have zero mean and unit variance.

##  Appendix G Quantifying Bias With Residuals

For general applications, investigating trends in residuals is a common test
for a model’s bias. If a residual resembles a ‘random’ noise, the bias of the
model is low. In an opposite case of trends in the residual, the model does
not fully represent the process in the input, i.e., it exhibits a bias.

For a denoising AE the situation is not straightforward and depends on the
perception of noise for each application. For a $2TCF$, in cases of low
detector count or frame-by-frame instabilities, the residual for the denoised
output should not have any trends (Fig. 13). However, in cases when an
average, ‘envelope’, dynamics is under interest, the dynamics heterogeneities
are treated as noise, leaving a pattern in the residual (Fig. 14), that can be
investigated separately. It is convenient to measure noise using
autocorrelation coefficient (ACC) at lag=1 of residuals, projected at one of
the time axes .

Figure 13: Example from training set where the residual does not have trends.
Figure 14: Example from training set where the residual has a well pronounced
trend.

Just like the latent space coordinates, the ACC of residuals for each new
example can be compared to the analogous values for examples in the training
set. This helps to understand, what values of ACC are unusual. For quantifying
the comparison, the probability density function (Fig. 15) is estimated for
ACC values of the training set and is calculated for each new test example.

Figure 15: Distribution of the model uncertainty (ACC(lag=1) of residuals) for
the training set. Solid blue line is the probability density function for the
distribution calculated with Gaussian Kernel density estimator.

##  Appendix H Examples of workflows for Online and Offline Data Analysis

Algorithm 1 (online analysis):

  * •

apply denoising AE to the $2TCF_{raw}$ to obtain $2TCF_{denoised}$

  * •

calculate an uncertainty score through latent space and/or residuals

  * •

unless observation is unusual:

    * –

do the first round of fit for $2TCF_{denoised}$ (optionally with wider time
cuts) using wide parameter boundaries

    * –

calculate new parameter boundaries based on values within _trust regions_

    * –

if needed, repeat the fits for $2TCF_{denoised}$ for single-frames cuts with
new parameter boundaries

  * •

if the observation is unusual:

    * –

use $2TCF_{raw}$ and wide-bin $1TCF$s for the fits

  * •

report results of the fit and _trust regions_

Algorithm 1 (offline analysis):

  * •

apply denoising AE to the $2TCF_{raw}$ to obtain $2TCF_{denoised}$

  * •

calculate an uncertainty score through latent space and/or residuals

  * •

based on the uncertainty score, alarm about an unusual observation if
encountered

  * •

unless observation is unusual:

    * –

do the first round of fit for $2TCF_{denoised}$ (optionally with wider time
cuts) using wide parameter boundaries

    * –

calculate new parameter boundaries based on values within _trust regions_

    * –

fit both $2TCF_{raw}$ and $2TCF_{denoised}$ with new parameter bounds

  * •

if the observation is unusual:

    * –

use $2TCF_{raw}$ and wide-bin $1TCF$s for the fits

    * –

report results of the fit and _trust regions_

## References

  * Bergmann  _et al._ [2009] U. Bergmann, J. Corlett, S. Dierker, R. Falcone, J. Galayda, M. Gibson, J. Hastings, B. Hettel, J. Hill, Z. Hussain, _et al._ ,  _Science and technology of future light sources: A white paper_ , Tech. Rep. (Stanford Linear Accelerator Center (SLAC), 2009).
  * Couprie [2014] M. E. Couprie, New generation of light sources: Present and future, Journal of Electron Spectroscopy and Related Phenomena 196, 3 (2014).
  * Vartanyants and Singer [2016] I. Vartanyants and A. Singer, Synchrotron light sources and free-electron lasers, edited by ej jaeschke, s. khan, jr schneider & jb hastings (2016).
  * Grybos  _et al._ [2015] P. Grybos, P. Kmon, P. Maj, and R. Szczygiel, 32k Channels readout IC for single photon counting detectors with 75 µm pitch, ENC of 123 e- rms, 9 e- rms offset spread and 2 % rms gain spread, in  _2015 IEEE Biomedical Circuits and Systems Conference (BioCAS)_ (2015) pp. 1–4.
  * Llopart  _et al._ [2002] X. Llopart, M. Campbell, R. Dinapoli, D. San Segundo, and E. Pernigotti, Medipix2: A 64-k pixel readout chip with 55-/spl mu/m square elements working in single photon counting mode, IEEE transactions on nuclear science 49, 2279 (2002).
  * Zhang  _et al._ [2018] Q. Zhang, E. M. Dufresne, and A. R. Sandy, Dynamics in hard condensed matter probed by x-ray photon correlation spectroscopy: Present and beyond, Current Opinion in Solid State and Materials Science 22, 202 (2018).
  * Ricci  _et al._ [2020] A. Ricci, G. Campi, B. Joseph, N. Poccia, D. Innocenti, C. Gutt, M. Tanaka, H. Takeya, Y. Takano, T. Mizokawa, _et al._ , Intermittent dynamics of antiferromagnetic phase in inhomogeneous iron-based chalcogenide superconductor, Physical Review B 101, 020508 (2020).
  * Evenson  _et al._ [2015] Z. Evenson, B. Ruta, S. Hechler, M. Stolpe, E. Pineda, I. Gallino, and R. Busch, X-ray photon correlation spectroscopy reveals intermittent aging dynamics in a metallic glass, Physical review letters 115, 175701 (2015).
  * Begam  _et al._ [2021] N. Begam, A. Ragulskaya, A. Girelli, H. Rahmann, S. Chandran, F. Westermeier, M. Reiser, M. Sprung, F. Zhang, C. Gutt, _et al._ , Kinetics of network formation and heterogeneous dynamics of an egg white gel revealed by coherent x-ray scattering, Physical Review Letters 126, 098001 (2021).
  * Yavitt  _et al._ [2020] B. M. Yavitt, D. Salatto, Z. Huang, Y. T. Koga, M. K. Endoh, L. Wiegart, S. Poeller, S. Petrash, and T. Koga, Revealing nanoscale dynamics during an epoxy curing reaction with x-ray photon correlation spectroscopy, Journal of Applied Physics 127, 114701 (2020).
  * Brown  _et al._ [1997] G. Brown, P. A. Rikvold, M. Sutton, and M. Grant, Speckle from phase-ordering systems, Physical Review E 56, 6601 (1997).
  * Madsen  _et al._ [2020] A. Madsen, A. Fluerasu, and B. Ruta, Structural dynamics of materials probed by x-ray photon correlation spectroscopy, Synchrotron Light Sources and Free-Electron Lasers: Accelerator Physics, Instrumentation and Science Applications , 1989 (2020).
  * Bikondoa [2017] O. Bikondoa, On the use of two-time correlation functions for x-ray photon correlation spectroscopy data analysis, Journal of applied crystallography 50, 357 (2017).
  * Madsen  _et al._ [2010] A. Madsen, R. L. Leheny, H. Guo, M. Sprung, and O. Czakkel, Beyond simple exponential correlation functions and equilibrium dynamics in x-ray photon correlation spectroscopy, New Journal of Physics 12, 055001 (2010).
  * Konstantinova  _et al._ [2021] T. Konstantinova, L. Wiegart, M. Rakitin, A. M. DeGennaro, and A. M. Barbour, Noise reduction in x-ray photon correlation spectroscopy with convolutional neural networks encoder-decoder models, Scientific Reports 11, 14756 (2021).
  * Abdar  _et al._ [2021] M. Abdar, F. Pourpanah, S. Hussain, D. Rezazadegan, L. Liu, M. Ghavamzadeh, P. Fieguth, X. Cao, A. Khosravi, U. R. Acharya, _et al._ , A review of uncertainty quantification in deep learning: Techniques, applications and challenges, Information Fusion (2021).
  * Lin  _et al._ [2021] C.-H. Lin, K. Dyro, O. Chen, D. Yen, B. Zheng, M. T. Arango, S. Bhatia, K. Sun, Q. Meng, L. Wiegart, _et al._ , Revealing meso-structure dynamics in additive manufacturing of energy storage via operando coherent x-ray scattering, Applied Materials Today 24, 101075 (2021).
  * Kingma and Ba [2014] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Branch  _et al._ [1999] M. A. Branch, T. F. Coleman, and Y. Li, A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems, SIAM Journal on Scientific Computing 21, 1 (1999).
  * Wohland  _et al._ [2001] T. Wohland, R. Rigler, and H. Vogel, The standard deviation in fluorescence correlation spectroscopy, Biophysical journal 80, 2987 (2001).
  * The HDF Group [2010] The HDF Group, Hierarchical data format version 5 (2000-2010).
  * Pezoa  _et al._ [2016] F. Pezoa, J. L. Reutter, F. Suarez, M. Ugarte, and D. Vrgoč, Foundations of json schema (2016).


[1,2,3]M. Faverzani These authors contributed equally to this work. These
authors contributed equally to this work. These authors contributed equally to
this work. These authors contributed equally to this work.

1]Physics Department, University of Milano - Bicocca, Milano, Italy 2]INFN
Sezione di Milano - Bicocca, Milano, Italy 3]Bicocca Quantum Technologies
(BiQuTe) Centre, Milano, Italy 4]Fondazione Bruno Kessler, Trento, Italy
5]INFN - TIFPA, Trento, Italy 6]Physics Department University of Salerno,
Salerno, Italy 7]INFN Sezione di Napoli, Napoli, Italy 8]Physics Department
University of Salento, Lecce, Italy 9]INFN Sezione di Lecce, Lecce, Italy
10]INRiM, Torino, Italy 11]INO-CNR BEC Center, Trento, Italy 12]Physics
Department University of Trento, Trento, Italy 13]INFN Laboratori Nazionali di
Frascati, Frascati, Italy 14]IFN-CNR, Trento, Italy 15]Department of Science
and Technology University of Sannio, Benevento, Italy

# Broadband parametric amplification in DARTWARS

<EMAIL_ADDRESS>P. Campana  R. Carobene  M. Gobbo  F. Ahrens  G.
Avallone  C. Barone  M. Borghesi  S. Capelli  G. Carapella  A. P. Caricato  L.
Callegaro  I. Carusotto  A. Celotto  A. Cian  A. D’Elia  D. Di Gioacchino  E.
Enrico  P. Falferi  L. Fasolo  E. Ferri  G. Filatrella  C. Gatti  D.
Giubertoni  V. Granata  C. Guarcello  A. Irace  D. Labranca  A. Leo  C. Ligi
G. Maccarrone  F. Mantegazzini  B. Margesin  G. Maruccio  R. Mezzena  A. G.
Monteduro  R Moretti  A. Nucciotti  L. Oberto  L. Origo  S. Pagano  A. S.
Piedjou Komnang  L. Piersanti  A. Rettaroli  S. Rizzato  S. Tocci  A. Vinante
M. Zannoni  A. Giachero  [  [  [  [  [  [  [  [  [  [  [  [  [  [  [

###### Abstract

Superconducting parametric amplifiers offer the capability to amplify feeble
signals with extremely low levels of added noise, potentially reaching
quantum-limited amplification. This characteristic makes them essential
components in the realm of high-fidelity quantum computing and serves to
propel advancements in the field of quantum sensing. In particular, Traveling-
Wave Parametric Amplifiers (TWPAs) may be especially suitable for practical
applications due to their multi-Gigahertz amplification bandwidth, a feature
lacking in Josephson Parametric Amplifiers (JPAs), despite the latter being a
more established technology. This paper presents recent developments of the
DARTWARS (Detector Array Readout with Traveling Wave AmplifieRS) project,
focusing on the latest prototypes of Kinetic Inductance TWPAs (KITWPAs). The
project aims to develop a KITWPA capable of achieving $20\,$ dB of
amplification. To enhance the production yield, the first prototypes were
fabricated with half the length and expected gain of the final device. In this
paper, we present the results of the characterization of one of the half-
length prototypes. The measurements revealed an average amplification of
approximately $9\,$dB across a $2\,$GHz bandwidth for a KITWPA spanning
$17\,$mm in length.

###### keywords:

Quantum noise, parametric amplifier, traveling wave, detector array readout,
qubits readout

##  1 Introduction

The rapid evolution of quantum computing demands for amplifiers featuring wide
bandwidth, large dynamic range, and quantum-limited added noise. Devices with
these characteristics would also offer significant advantages to quantum
sensing for the readout of low-temperature detectors, ultimately enhancing the
performance of next-generation experiments.

Superconducting parametric amplifiers exploit a non-linear circuit to amplify
signals. Among these devices, Josephson Parametric Amplifiers (JPAs) [1] stand
out as the best-established technology. While JPAs have demonstrated their
ability to achieve quantum-limited noise performances, they are constrained by
other limitations, primarily associated with their relatively narrow bandwidth
and limited dynamic range of amplification.

Another class of devices comprises Traveling Wave Parametric Amplifiers
(TWPAs), which have the potential to expand the bandwidth and dynamic range to
values close to the ones of existing standard commercial amplifiers. A TWPA
features a non-linear transmission line that is either embedded with Josephson
junctions (TWJPA) [2] or made with high-kinetic inductance materials (KITWPA)
[3]. Both of these approaches are being pursued [4, 5, 6, 7, 8] by the
DARTWARS (Detector Array Readout with Traveling Wave AmplifieRS) project,
which aims to reach gain, noise and saturation power of around $20\,$dB,
$600\,$mK, and $-50\,$dBm respectively.

This paper focuses on the latest developments of DARTWARS, providing a more
comprehensive characterization of the KITWPA prototypes introduced in [9],
which demonstrated significant amplification in preliminary tests. A KITWPA
consists of a weakly dispersive phase-matched transmission line, controlled by
dispersion engineering, which is required to create exponential amplification
[10]. The prototype characterized in this paper is an artificial transmission
line organized in multiple cells, each one composed of a stub-loaded coplanar
waveguide (CPW) with an inductance ($L$) dependent on the kinetic inductance
of the material and two interdigitated capacitors (IDC) to establish
capacitance ($C$) towards the ground. Dispersion engineering was achieved by
organizing the transmission line into $N_{sc}=523$ super-cells, each
comprising $60$ unloaded cells with a characteristic impedance
$Z_{0}=\sqrt{{L}/{C}}=50\,\Omega$, and $6$ loaded cells with
$Z_{0}=80\,\Omega$, for a total of $L_{sc}=66$ cells in each supercell (fig.
2).

The $17\,$mm transmission line length corresponds to an expected amplification
of approximately $10\,$dB. The design of the KITWPAs was optimized to produce
these devices with a bandwidth centered at $4-5\,$GHz. The next fabrication
will produce devices with a higher ratio of loaded to unloaded cells, shifting
the bandwidth center up to $6\,$GHz. The KITWPAs were produced by depositing a
thin film of NbTiN on a high-resistivity silicon wafer (fig. 2). Further
details on the fabrication process are reported in [9].

Figure 1: The diagram (not to scale) of the lumped element artificial
transmission line. The dark blue areas indicate the high-resistivity
superconductor, whereas the light blue areas stand for the substrate.

Figure 2: Scanning electron microscope images of the patterned artificial
transmission line. Dark areas represent the NbTiN film, whereas light areas
are the Si substrate.

Figure 3: A schematic diagram of the experimental setup. At room temperature,
the VNA and pump rf signals are combined using a power combiner. A dc current
produced by an isolated voltage source is then added through a bias tee. The
resulting rf/dc combined signal is routed inside the cryostat to the KITWPA
located on the mixing chamber at $60\,$mK. The output signal is terminated
with a $50\,\Omega$ load using a second bias tee at room temperature before
the input port of the VNA.

##  2 Characterization of the KITWPA prototype

In the present work, one of the devices with $I_{c}>1\,$mA and a clearly
identifiable stop-band was selected for low-temperature characterization at
$60\,$mK. The experimental setup shown in fig. 3 was used to scan the $S_{21}$
response of the TWPA while controlling its dc bias and pump signal. The
instruments used for performing the characterization were controlled using the
Qtics Python package [11]. The KITWPA chip was mounted in a general-purpose
copper box, equipped with CPW launchers to feed the device.

###  2.1 Measurement of non-linearity parameters

At the lowest order, the inductance of the KITWPA increases quadratically as a
function of the dc bias $I$. This non-linear response can be measured by
considering the phase of the $S_{21}$ parameter at a fixed frequency $f$:

$\frac{\theta(I)-\theta_{0}}{\theta_{r}}=-\frac{1}{2}\left(\frac{I}{I^{*}}\right)^{2},$
(1)

where $\theta_{0}$ is the phase at zero bias, while $I^{*}$ is an intrinsic
parameter of the material known as the scaling current. Finally,
$\theta_{r}=2\pi f\Delta t_{\text{TWPA}}$ represents the phase delay over the
entire transmission line.

The time required by a signal to travel across the device $\Delta
t_{\text{TWPA}}=33\,$ns was measured with a Time-Domain-Reflectometer (TDR).
This allowed to estimate the phase delay and the phase velocity
$v_{p}=(N_{\text{sc}}L_{\text{sc}})/\Delta t_{\text{TWPA}}$. The latter was
also cross-checked with an alternative estimate based on the relation
$v_{p}=\lambda_{sb}\,f_{sb}$, where $\lambda_{sb}=2L_{sc}$ and
$f_{sb}=7.7\,$GHz refer respectively to the wavelength corresponding to the
spatial separation between subsequent periodic loadings and the central
frequency of the measured stop-band. These approaches yielded a phase-velocity
of $v_{p}=1016\,\text{cells/ns}$ ($1.67\cdot 10^{-3}\,c$) and
$v_{p}=1030\,\text{cells/ns}$ ($1.69\cdot 10^{-3}\,c$) respectively,
demonstrating consistent results between the two techniques.

To measure the scaling current $I^{*}$, the data were acquired by recording
the $S_{21}$ response with the VNA for various dc bias values while keeping
the rf pump turned off. As shown in fig. 5, fitting the measured and rescaled
phase at a fixed frequency with eq. 1 allows to determine $I^{*}$. Repeating
the fit for different frequencies over the whole amplification bandwidth
resulted in an average value of $I^{*}=(6.72\pm 0.03)\,$mA, where the
uncertainty was estimated as the standard deviation of the fit results.

The critical current $I_{c}$ of the device was instead measured by considering
the amplitude of $S_{21}$ as a function of the dc bias. For a fixed frequency,
$I_{c}$ can be easily identified as a significant discontinuity in the
transmitted VNA signal amplitude (fig. 5). Consistent values of $I_{c}$ were
found when considering different frequencies across the entire amplification
bandwidth, resulting in $I_{c}=(1.41\pm 0.01)\,$mA.

In a preliminary characterization [9], another batch of prototypes was tested
in liquid helium, which led to a measured scaling and critical currents of
$I^{*}=(5.3\pm 0.1)\,$mA and $I_{c}=(1.5\pm 0.2)\,$mA respectively. Comparing
these measurements with our device, the critical currents are in good
agreement while the scale currents have slightly different values, but still
maintain the same order of magnitude.

Figure 4: Measured relative phase of a signal transmitted through the TWPA
with the pump off as a function of bias current. The quadratic fit allows to
estimate the scaling current $I^{*}$.

Figure 5: Measured amplitude of $S_{21}$ with the pump off as a function of
bias current. The critical current $I_{c}$ is identified by a sharp
discontinuity in the values.

###  2.2 Gain measurements

The best gain and corresponding bandwidth were characterized through a VNA
frequency sweep in the $1-8\,$GHz range, while varying the pump rf signal
within the $7.9-8.3\,$GHz range, which is just above the stop-band. Additional
parameters, such as the dc biasing current or the pump power, were fine-tuned
separately to $1.2\,$ mA and $-39.4\,$ dBm. In fig. 7, the measured gain is
plotted as a function of the VNA and the pump rf signals frequencies. As
predicted by the phase-matching relations [12], a splitting of the bandwidth
becomes apparent as the pump rf signal frequency increases. A mean gain of
$9.11\,$dB was achieved over a $3-5\,$GHz range at a pump rf signal frequency
of $7.99\,$GHz (fig. 7). However, the gain does not appear to be uniform due
to large ripples (inset in fig. 7). By performing a Fast Fourier Transform
(FFT) of the transmitted signal, it can be seen that the dominant frequency
corresponds to a wavelength approximately twice the length of the device. This
explains the ripples as resulting from impedance mismatches introduced by the
line adaptation and connections between the KITWPA and the packaging, which
lead to the creation of standing waves across the device.

Figure 6: Measured gain as a function of VNA rf signal frequency and pump
frequency. Near a pump frequency of $8.1\,$GHz, a bandwidth splitting can be
observed

Figure 7: Maximal gain profile in a $1-7\,$GHz frequency range. The apparently
noisy profile is instead produced by large ripples, as shown in the inset.

The final measurements focused on assessing the power handling capabilities of
the device. While keeping the pump and dc bias parameters fixed, the mean gain
was evaluated as a function of the VNA rf signal power injected into the TWPA.
While the gain remains constant for low VNA powers, an increase in power
causes the TWPA to saturate its output, leading to a gradual degradation of
the gain. It is common to extract the $1\,$dB compression point, which is
identified as the input power necessary to deviate by $1$ dB from the maximum
stable gain. The measured mean gain as a function of the TWPA input power is
shown in fig. 8. A value of $P_{-1\text{dB}}=-47.4\,$dB was obtained
considering the mean gain in the $3-5\,$GHz frequency range, which
demonstrates the significant power handling capabilities of these devices with
respect to other parametric amplifiers.

Figure 8: Measured gain as a function of the input signal power. From this
plot, it is possible to extrapolate the $1\,$dB compression point
$P_{-1\,\text{dB}}$, defined as the power where the gain has decreased of
$1\,$dB in respect to its nominal value, due to saturation.

##  3 Conclusions

This paper presented the latest developments of the DARTWARS project,
regarding the advancements of KITWPAs for applications in quantum computing
and sensing. We characterized one of the latest prototypes at cryogenic
temperatures ($T\sim 60\,$mK) and measured a maximum mean gain of $9.11\,$dB
in the frequency range from $3$ to $5\,$GHz. The examined device exhibited a
critical current of $I_{c}=(1.41\pm 0.01)\,$mA and a scaling current of
$I^{*}=(6.72\pm 0.03)\,$mA, with its $1\,$dB compression point estimated at
$P_{-1\text{dB}}=-47.4\,$dB. The measured values are in reasonable agreement
with the design goals of the prototype. Future phases of the project will
involve noise performance characterization, producing additional short
prototypes with an amplification bandwidth centered at higher frequencies, and
eventually fabricating full-length devices with $20\,$dB of amplification.

## Acknowledgement

This work is supported by DARTWARS, a project funded by the Italian Institute
of Nuclear Physics (INFN) within the Technological and Interdisciplinary
Research Commission (CSN5), by European Union’s H2020-MSCA Grant Agreement No.
101027746, by the Italian National Centre for HPC Big Data and Quantum
Computing (PNRR MUR project CN0000013-ICSC) and by the Italian National
Quantum Science and Technology Institute (PNRR MUR project PE0000023-NQSTI).
We acknowledge the support of the FBK cleanroom team for the fabrication. We
also acknowledge useful discussions with Jiansong Gao, Michael Vissers, Jordan
Wheeler, and Maxime Malnou. CB and SP acknowledge support from University of
Salerno - Italy under the projects FRB19PAGAN, FRB20BARON and FRB22PAGAN.

## References

  *   * Castellanos-Beltran et al. [2008] Castellanos-Beltran, M.A., Irwin, K.D., Hilton, G.C., Vale, L.R., Lehnert, K.W.: Amplification and squeezing of quantum noise with a tunable Josephson metamaterial. Nature Physics 4, 929 (2008) https://doi.org/10.1038/nphys1090
  * Macklin et al. [2015] Macklin, C., et al.: A near quantum-limited Josephson traveling-wave parametric amplifier. Science 350(6258), 307–310 (2015) https://doi.org/10.1126/science.aaa8525
  * Ho Eom et al. [2012] Ho Eom, B., et al.: A wideband, low-noise superconducting amplifier with high dynamic range. Nature Physics 8(8), 623–627 (2012) https://doi.org/10.1038/nphys2356
  * Borghesi et al. [2023] Borghesi, M., et al.: Progress in the development of a KITWPA for the DARTWARS project. Nucl. Instrum. Meth. A 1047, 167745 (2023) https://doi.org/10.1016/j.nima.2022.167745 arXiv:2208.10101 [quant-ph]
  * Rettaroli et al. [2023] Rettaroli, A., et al.: Ultra low noise readout with traveling wave parametric amplifiers: The DARTWARS project. Nucl. Instrum. Meth. A 1046, 167679 (2023) https://doi.org/10.1016/j.nima.2022.167679 arXiv:2207.12775 [quant-ph]
  * Pagano et al. [2022] Pagano, S., et al.: Development of Quantum Limited Superconducting Amplifiers for Advanced Detection. IEEE Transactions on Applied Superconductivity 32(4), 1–5 (2022) https://doi.org/10.1109/TASC.2022.3145782
  * Granata et al. [2023] Granata, V., et al.: Characterization of Traveling-Wave Josephson Parametric Amplifiers at T = 0.3 K. IEEE Transactions on Applied Superconductivity 33(1), 1–7 (2023) https://doi.org/10.1109/TASC.2022.3214656
  * Guarcello et al. [2023] Guarcello, C., et al.: Modeling of josephson traveling wave parametric amplifiers. IEEE Transactions on Applied Superconductivity 33(1), 1–7 (2023) https://doi.org/10.1109/TASC.2022.3214751
  * Mantegazzini et al. [2023] Mantegazzini, F., et al.: High kinetic inductance NbTiN films for quantum limited travelling wave parametric amplifiers. Physica Scripta 98(12), 125921 (2023) https://doi.org/10.1088/1402-4896/ad070d
  * Chaudhuri et al. [2017] Chaudhuri, S., et al.: Broadband parametric amplifiers based on nonlinear kinetic inductance artificial transmission lines. Applied Physics Letters 110(15), 152601 (2017) https://doi.org/10.1063/1.4980102
  * Carobene et al. [2024] Carobene, R., Campana, P., Gobbo, M., Giachero, A.: Qtics - Quantum Technologies Instrumentation ControlS (Version 0.0.2) (2024) https://doi.org/10.5281/zenodo.10450507
  * Malnou et al. [2021] Malnou, M., Vissers, M.R., Wheeler, J.D., Aumentado, J., Hubmayr, J., Ullom, J.N., Gao, J.: Three-Wave Mixing Kinetic Inductance Traveling-Wave Amplifier with Near-Quantum-Limited Noise Performance. PRX Quantum 2, 010302 (2021) https://doi.org/10.1103/PRXQuantum.2.010302


$\displaystyle\equiv\frac{e^{-\bm{x}^{2}/2\sigma_{I}^{2}}}{(2\pi\sigma_{I}^{2})^{d/2}}I(t)$
Suppose also that we have linearized response function $\chi^{\lambda J}$ of
$\lambda$ to $I(t)$ of the form: $\langle\lambda(\bm{k},\omega)\rangle\equiv
e^{-k^{2}\sigma_{I}^{2}/2}\chi^{\lambda J}(\bm{k},\omega)I(\omega)$ Then the
response function of $O(t)$ with respect to $I(t)$ is given by:
$\chi^{OI}(\omega)\approx\frac{\sigma_{O}^{d}}{(2\pi
r^{2})^{\frac{d}{2}}}\int\mathrm{d}^{d}\bm{u}\,e^{i\bm{u}\cdot\bm{\hat{r}}}\chi^{\lambda
J}\left(\frac{\bm{u}}{r},\omega\right)$ (98) where the approximation holds
when $\sigma_{I}$, $\sigma_{O}\ll r$. When $\chi^{\lambda J}$ has rotational
symmetry, the result for $d=2$ becomes:
$\displaystyle\chi^{OI}(\omega)\approx\frac{\sigma_{O}^{2}}{r^{2}}\int_{0}^{\infty}\mathrm{d}u\,uJ_{0}(u)\chi^{\lambda
J}\left(\frac{u}{r},\omega\right)$ (99) where $J_{0}$ is a Bessel function.
The result for $d=3$ becomes:
$\displaystyle\chi^{OI}(\omega)\approx\sqrt{\frac{2}{\pi}}\frac{\sigma_{O}^{3}}{r^{3}}\int_{0}^{\infty}\mathrm{d}u\,u\sin
u\,\chi^{\lambda J}\left(\frac{u}{r},\omega\right)$ (100)

Proof: We start by writing down $O(\omega)$ and expanding
$\lambda(\bm{x},\omega)$ into $k$ modes:

$\displaystyle\langle O(\omega)\rangle$
$\displaystyle=\int\mathrm{d}^{d}\bm{x}\,e^{-(\bm{x}-\bm{r})^{2}/2\sigma_{O}^{2}}\langle\lambda(\bm{x},\omega)\rangle$
$\displaystyle=\int\frac{\mathrm{d}^{d}\bm{x}\mathrm{d}^{d}\bm{k}}{(2\pi)^{d}}\frac{e^{i\bm{k}\cdot\bm{x}}}{e^{(\bm{x}-\bm{r})^{2}/2\sigma_{O}^{2}}}\langle\lambda(\bm{k},\omega)\rangle$

Solving the $\bm{x}$ integral yields

$\displaystyle\langle O(\omega)\rangle$
$\displaystyle=\frac{(2\pi\sigma_{O}^{2})^{d/2}}{(2\pi)^{d}}\int\mathrm{d}^{d}\bm{k}\frac{e^{i\bm{k}\cdot\bm{r}}}{e^{k^{2}\sigma_{O}^{2}/2}}\langle\lambda(\bm{k},\omega)\rangle$

Then we insert the linear response formula for
$\langle\lambda(\bm{k},\omega)\rangle$:

$\displaystyle\langle O(\omega)\rangle$
$\displaystyle=\frac{(2\pi\sigma_{O}^{2})^{d/2}}{(2\pi)^{d}}\int\mathrm{d}^{d}\bm{k}\frac{e^{i\bm{k}\cdot\bm{r}}}{e^{k^{2}(\sigma_{I}^{2}+\sigma_{O}^{2})/2}}\chi^{\lambda
J}(\bm{k},\omega)I(\omega)$

We then change to the dimensionless integration variable
$\bm{u}=|\bm{r}|\bm{k}$ and find:

$\displaystyle\langle O(\omega)\rangle$
$\displaystyle=\frac{\sigma_{O}^{d}}{(2\pi
r^{2})^{d/2}}\int\mathrm{d}^{d}\bm{u}\frac{e^{i\bm{u}\cdot\bm{\hat{r}}}}{e^{\frac{u^{2}}{2}\left(\frac{\sigma_{I}^{2}}{r^{2}}+\frac{\sigma_{O}^{2}}{r^{2}}\right)}}\chi^{\lambda
J}\left(\frac{\bm{u}}{r},\omega\right)I(\omega)$

where $r=|\bm{r}|$. The linear response function $\chi^{OI}$ of $O(\omega)$ to
$I(\omega)$, defined by $\langle O(\omega)\rangle=\chi^{OI}(\omega)I(\omega)$,
is everything in the above expression except $I(\omega)$.

Whenever $\sigma_{I}$, $\sigma_{O}$ $\ll r$, we can drop the exponential term
$e^{u^{2}\sigma^{2}/2r^{2}}$. This is a very weak assumption: since
$\sigma_{I}$ and $\sigma_{O}$ are the sender and receiver sizes, the task of
sending information is not well defined when $\sigma_{I}$, $\sigma_{O}$ $\geq
r$. From this approximation we obtain our desired result:

$\displaystyle\chi^{OI}(\omega)$ $\displaystyle=\frac{\sigma_{O}^{d}}{(2\pi
r^{2})^{d/2}}\int\mathrm{d}^{d}\bm{u}\,e^{i\bm{u}\cdot\bm{\hat{r}}}\chi^{\lambda
J}\left(\frac{\bm{u}}{r},\omega\right)$

In $d=2$ dimensions, when $\chi^{\lambda J}(\bm{k},\omega)$ is rotationally
invariant, we can convert to polar coordinates and perform the $\theta$
integral:

$\displaystyle\chi^{OI}(\omega)$ $\displaystyle=\frac{\sigma_{O}^{2}}{2\pi
r^{2}}\int\mathrm{d}u\,u\chi^{\lambda
J}\left(\frac{u}{r},\omega\right)\int_{0}^{2\pi}\mathrm{d}\theta\,e^{iu\cos\theta}$
$\displaystyle=\frac{\sigma_{O}^{2}}{r^{2}}\int_{0}^{\infty}\mathrm{d}u\,uJ_{0}(u)\chi^{\lambda
J}\left(\frac{u}{r},\omega\right)$

where $J_{0}(u)$ is a Bessel function. Likewise in $d=3$ dimensions (with
rotational invariance), the result is obtained by moving to spherical
coordinates:

$\displaystyle\chi^{OI}(\omega)$
$\displaystyle=\frac{2\pi\sigma_{O}^{3}}{(2\pi
r^{2})^{3/2}}\int_{0}^{\infty}\mathrm{d}u\,u^{2}\chi^{\lambda
J}\left(\frac{u}{r},\omega\right)\int_{0}^{\pi}\mathrm{d}\theta\,\sin\theta
e^{iu\cos\theta}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\frac{\sigma_{O}^{3}}{r^{3}}\int_{0}^{\infty}\mathrm{d}u\,u\sin
u\,\chi^{\lambda J}\left(\frac{u}{r},\omega\right)$

###  E.2 Output fluctuations from density fluctuations

Lemma: $S^{O}(\omega)$ formula Consider a system in $d$ dimensions with a
density field $\lambda$ with mean zero and an output signal $O(t)$ defined as
the density integrated over a Gaussian shell around a point $\bm{r}$:
$\displaystyle
O(t)\equiv\int\mathrm{d}^{d}\bm{x}\,e^{-(\bm{x}-\bm{r})/2\sigma_{O}^{2}}\lambda(\bm{x},t)$
Then the power spectrum of $O(t)$ is related to the power spectrum of
$\lambda(\bm{x},t)$ by:
$S^{O}(\omega)=\sigma_{O}^{d}\int\mathrm{d}^{d}\bm{u}\,e^{-u^{2}}S^{\lambda}\left(\frac{\bm{u}}{\sigma_{O}},\omega\right)$
(101) When $S^{\lambda}$ has rotational symmetry, the result for $d=2$
trivially becomes:
$S^{O}(\omega)=2\pi\sigma_{O}^{2}\int_{0}^{\infty}\mathrm{d}u\,ue^{-u^{2}}S^{\lambda}\left(\frac{u}{\sigma_{O}},\omega\right)$
(102) For $d=3$ it becomes
$S^{O}(\omega)=4\pi\sigma_{O}^{3}\int_{0}^{\infty}\mathrm{d}u\,u^{2}e^{-u^{2}}S^{\lambda}\left(\frac{u}{\sigma_{O}},\omega\right)$
(103)

Proof: We start with step 2 from the proof of lemma (98):

$\displaystyle O(\omega)$
$\displaystyle=\frac{\sigma_{O}^{d}}{(2\pi)^{d/2}}\int\mathrm{d}^{d}\bm{k}\frac{e^{i\bm{k}\cdot\bm{r}}}{e^{k^{2}\sigma_{O}^{2}/2}}\lambda(\bm{k},\omega)$

Inserting this into the expectation value $\langle
O(\omega)O(\omega^{\prime})\rangle$, we find:

$\displaystyle\langle O(\omega)O(\omega^{\prime})\rangle$
$\displaystyle=\frac{\sigma_{O}^{2d}}{(2\pi)^{d}}\int\mathrm{d}^{d}\bm{k}\mathrm{d}^{d}\bm{k}^{\prime}\frac{e^{i(\bm{k}+\bm{k}^{\prime})\cdot\bm{r}}}{e^{\frac{\sigma_{O}^{2}}{2}(k^{2}+k^{\prime
2})}}\langle\lambda(\bm{k},\omega)\lambda(\bm{k}^{\prime},\omega^{\prime})\rangle$

Next we insert the identity $\langle
f(\omega)f(\omega^{\prime})\rangle=2\pi\delta(\omega+\omega^{\prime})S^{f}(\omega)$
(see eq. 2.54 in [54]) to both sides of the above expression (where we apply
it to both $\bm{k}$ and $\omega$ on the right hand side):

$\displaystyle S^{O}(\omega)$
$\displaystyle=\sigma_{O}^{2d}\int\mathrm{d}^{d}\bm{k}\mathrm{d}^{d}\bm{k}^{\prime}\frac{e^{i(\bm{k}+\bm{k}^{\prime})\cdot\bm{r}}}{e^{\frac{\sigma_{O}^{2}}{2}(k^{2}+k^{\prime
2})}}\delta^{d}(\bm{k}+\bm{k}^{\prime})S^{\lambda}(\bm{k},\omega)$
$\displaystyle=\sigma_{O}^{2d}\int\mathrm{d}^{d}\bm{k}\,e^{-k^{2}\sigma_{O}^{2}}S^{\lambda}(\bm{k},\omega)$

We then change to the dimensionless integration variable
$\bm{u}=\sigma_{O}\bm{k}$, and find the desired result:

$\displaystyle S^{O}(\omega)$
$\displaystyle=\sigma_{O}^{d}\int\mathrm{d}^{d}\bm{u}\,e^{-u^{2}}S^{\lambda}\left(\frac{\bm{u}}{\sigma_{O}},\omega\right)$

###  E.3 Dissipation kernel from response function

Lemma: $\mathcal{D}(\omega)$ formula This gives a formula for the dissipation
kernel $\mathcal{D}(\omega)$, which is essentially a generalization of the
real part of electrical impedance for arbitrary systems. Consider a system in
$d$ dimensions with a density field $\lambda$ with mean zero and an input
signal $J(\bm{x},t)$ defined via a Gaussian smeared density current:
$J(\bm{x},t)=\frac{e^{-\bm{x}^{2}/2\sigma_{I}^{2}}}{(2\pi\sigma_{I}^{2})^{d/2}}I(t)$
Suppose that we have linearized response function $\chi^{\lambda J}$ of
$\lambda$ to $J$, defined by
$\langle\lambda(\bm{k},\omega)\rangle\equiv\chi^{\lambda
J}(\bm{k},\omega)J(\bm{k},\omega)$, which is even in $\bm{k}$ and Hermitian in
$\omega$: $\displaystyle\chi^{\lambda J}(-u,\omega)=\chi^{\lambda
J}(u,\omega)$ $\displaystyle\chi^{\lambda
J}(u,-\omega)=\overline{\chi^{\lambda J}}(u,\omega)$ And finally, suppose that
we can write the dissipative work required to produce this input signal in
terms of some factor $\mathcal{C}$, the input current, and the density field:
$\langle\textnormal{Diss}(t)\rangle=\mathcal{C}\int\mathrm{d}^{d}\bm{x}\,\langle
J(\bm{x},t)\lambda(\bm{x},t)\rangle$ Then the dissipation kernel
$\mathcal{D}(\omega)$ of $S^{I}(\omega)$ is given by
$\begin{split}\langle\textnormal{Diss}\rangle&=\int\mathrm{d}\omega\,\mathcal{D}(\omega)S^{I}(\omega)\\\
\mathcal{D}(\omega)&=\frac{\mathcal{C}}{\sigma_{I}^{d}}\int\frac{\mathrm{d}^{d}\bm{u}\,e^{-u^{2}}}{(2\pi)^{d+1}}\,\textnormal{Re}\left[\chi^{\lambda
J}\left(\frac{\bm{u}}{\sigma_{I}},\omega\right)\right]\end{split}$ (104)

Proof: We start by expanding in $\bm{k}$ modes and using Parseval’s theorem:

$\displaystyle\langle\textnormal{Diss}\rangle$
$\displaystyle=\mathcal{C}\int\mathrm{d}^{d}\bm{x}\,\langle
J(\bm{x},t)\lambda(\bm{x},t)\rangle=\mathcal{C}\int\frac{\mathrm{d}^{d}\bm{k}}{(2\pi)^{d}}\langle
J(\bm{k},t)\lambda(-\bm{k},t)\rangle$
$\displaystyle=\mathcal{C}\int\frac{\mathrm{d}^{d}\bm{k}\mathrm{d}\omega\mathrm{d}\omega^{\prime}}{(2\pi)^{d}(2\pi)^{2}}e^{i(\omega+\omega^{\prime})t}\langle
J(\bm{k},\omega)\lambda(-\bm{k},\omega^{\prime})\rangle$

We now focus on the Fourier space term, which we evaluate by inserting the
linear response function:

$\displaystyle\langle J(\bm{k},\omega)\lambda(-\bm{k},\omega^{\prime})\rangle$
$\displaystyle=\langle
J(\bm{k},\omega)J(-\bm{k},\omega^{\prime})\rangle\chi^{\lambda
J}(-\bm{k},\omega^{\prime})$ $\displaystyle=\langle
I(\omega)I(\omega^{\prime})\rangle e^{-k^{2}\sigma_{I}^{2}}\chi^{\lambda
J}(\bm{k},\omega^{\prime})$

where we have used that $\chi$ is even in $\bm{k}$ and the Fourier transform
of the Gaussian $J(\bm{k},\omega)=e^{-\bm{k}^{2}\sigma_{I}^{2}/2}I(\omega)$.
Next we use the identity $\langle
I(\omega)I(\omega^{\prime})\rangle=2\pi\delta(\omega+\omega^{\prime})S^{I}(\omega)$
(see eq. 2.54 in [54]):

$\displaystyle\langle J(\bm{k},\omega)\lambda(-\bm{k},\omega^{\prime})\rangle$
$\displaystyle=2\pi\delta(\omega+\omega^{\prime})S^{I}(\omega)e^{-k^{2}\sigma_{I}^{2}}\chi^{\lambda
J}(\bm{k},-\omega)$

Plugging this into our original expression, we find:

$\displaystyle\langle\textnormal{Diss}\rangle$
$\displaystyle=\mathcal{C}\int\mathrm{d}\omega\,S^{I}(\omega)\int\frac{\mathrm{d}^{d}\bm{k}}{(2\pi)^{d+1}}e^{-k^{2}\sigma_{I}^{2}}\chi^{\lambda
J}(\bm{k},-\omega)$

To make the integral dimensionless we change variables to
$\bm{u}=\sigma_{I}\bm{k}$,

$\displaystyle\langle\textnormal{Diss}\rangle$
$\displaystyle=\frac{\mathcal{C}}{\sigma_{I}^{d}}\int\mathrm{d}\omega\,S^{I}(\omega)\int\frac{\mathrm{d}^{d}\bm{u}}{(2\pi)^{d+1}}e^{-u^{2}}\chi^{\lambda
J}\left(\frac{\bm{u}}{\sigma_{I}},-\omega\right)$

Finally, we make the observation that if $I(t)$ is real, one can check that
$S^{I}(\omega)$ is even in $\omega$. From this and the fact that
$\chi^{\lambda J}$ is Hermitian in $\omega$, we can see that only the real
portion of $\chi^{\lambda J}$ will contribute to the dissipation:

$\displaystyle\int_{-\infty}^{\infty}\mathrm{d}\omega\,S^{I}(\omega)\chi^{\lambda
J}(-\omega)$
$\displaystyle=\int_{-\infty}^{\infty}\mathrm{d}\omega\,S^{I}(\omega)\textnormal{Re}\left[\chi^{\lambda
J}(\omega)\right]$

This allows us to write down the dissipation kernel $\mathcal{D}(\omega)$,
which is analogous to the real part of electrical impedance:

$\displaystyle\mathcal{D}(\omega)$
$\displaystyle=\frac{\mathcal{C}}{\sigma_{I}^{d}}\int\frac{\mathrm{d}^{d}\bm{u}}{(2\pi)^{d+1}}e^{-u^{2}}\textnormal{Re}\left[\chi^{\lambda
J}\left(\frac{\bm{u}}{\sigma_{I}},\omega\right)\right]$