# Surrogate-data-enriched Physics-Aware Neural Networks

Raphael Leiteritz  Institute for Parallel and Distributed Systems
University of Stuttgart, Germany  Patrick Buchfink  Institute of Applied
Analysis and Numerical Simulation
University of Stuttgart, Germany  Bernard Haasdonk  Institute of Applied
Analysis and Numerical Simulation
University of Stuttgart, Germany  Dirk Pflüger  Institute for Parallel and
Distributed Systems
University of Stuttgart, Germany

###### Abstract

Neural networks can be used as surrogates for PDE models. They can be made
physics-aware by penalizing underlying equations or the conservation of
physical properties in the loss function during training. Current approaches
allow to additionally respect data from numerical simulations or experiments
in the training process. However, this data is frequently expensive to obtain
and thus only scarcely available for complex models. In this work, we
investigate how physics-aware models can be enriched with computationally
cheaper, but inexact, data from other surrogate models like Reduced-Order
Models (ROMs). In order to avoid trusting too-low-fidelity surrogate
solutions, we develop an approach that is sensitive to the error in inexact
data. As a proof of concept, we consider the one-dimensional wave equation and
show that the training accuracy is increased by two orders of magnitude when
inexact data from ROMs is incorporated.

##  1 Introduction

Design, optimization or control of complex phenomena are tasks that are
critical for applications such as CO2 storage, e.g. [6], or biomechanical
simulations, e.g. [15]. For computationally expensive high-fidelity
simulations, such tasks are prohibitive. With computationally efficient
surrogate models these tasks can be carried out in an approximative fashion.

Machine Learning techniques can be used to derive such surrogates. Neural
networks are one class of such data-driven methods that have successfully been
applied to the solution of partial differential equations (PDEs) in various
settings [7, 11, 12]. In a recent work [8], it has even been demonstrated that
a data-driven method can outperform a numerical method in both, accuracy and
runtime, to solve an inverse uncertainty quantification problem.

Already a few decades ago, methods have been proposed that use neural networks
to solve PDEs by constraining the loss function with the underlying equations
[7]. In recent years, this original idea has seen a renaissance in the form of
so-called Physics-Informed Neural Networks (PINNs) [12]. These have meanwhile
been successfully applied to a variety of problems such as reconstructing
pressure and velocities from visual flow data, simulating blood-flow in
cardiovascular structures [13] or subsurface flow [14]. In contrast to methods
learning directly on simulation data [4] using e.g. CNNs [10] or LSTMs [9],
PINNs add a term to the loss function which penalizes predictions that do not
satisfy the underlying PDE. In the scope of this paper, we differentiate these
loss terms by their nature. The _physical (or BVP) loss_ aims to minimize (a)
the PDE residual on interior data points and (b) the (initial) boundary data
on boundary data points. PINNs featuring an additional _data loss_ on interior
data points are referred to as _data-enriched_ PINNs.

So far data-enriched PINNs in literature are based on expensive measurements
which are either obtained from real or numerical experiments. An example is
[12] where PINNs are trained on experimental data and afterwards are used to
estimate model parameters of the PDE to solve this inverse problem.

In this work, we aim to integrate comparably cheap data from surrogate models
in the data loss of the PINN. In a slightly more general view, we call this
data _inexact data_ as solutions from surrogates are approximative solutions.
We additionally assume that this inexact data is provided with an error bound
that quantifies the error with respect to the exact solution. This setting is
quite natural for many surrogate models like e.g. Reduced-Order Models (ROMs).
As our main contribution, we propose the notion of _error-sensitive_ PINNs
(see Fig. 1) as a generalization of data-enriched PINNs. The core idea is that
the error bound from the inexact data is taken into account during the
training by relaxing the optimization goal if the error with respect to the
inexact data is smaller than the error bound. This approach comes with two key
advantages: Firstly, the additional knowledge on the solution within the
solution domain may provide a boost to training times as well as prediction
accuracy as it now offers the optimizer more data to find a correct solution.
This is crucial since, in a simulation setting, data is usually scarce due to
their high computational costs. Secondly, since surrogates can be of low
fidelity, the error-sensitive part does not force the PINN to fit the inexact
data but instead gives it leeway to improve over the inexact data. Thereby,
PINNs are encouraged to refine the given inexact data.

Figure 1: Illustration of different PINN approaches within this paper.

Previous studies for PINNs have shown [17, 16] that the prediction accuracy is
sensitive to the weighting of the different loss terms. In order to apply [16]
to data-enriched PINNs, we generalize this idea to the case of more than two
loss terms. Moreover, we include a comparison of the different weightings in
our numerical experiments.

The performance of PINNs on inexact data and the new error-sensitive approach
are studied in a numerical experiment based on the one-dimensional linear wave
equation. We show that the error-sensitive PINN outperforms a standard non-
data-enriched PINN. Moreover, the experiments show that weighting the losses
correctly is essential.

The rest of the paper is structured as follows: In Section 2, we introduce the
essentials of scientific machine learning for PDEs using PINNs. Subsequently,
we present the loss weighting strategies, the error-sensitive PINNs and ROM-
data-enriched PINNs in Section 3. The new methods are compared to classical
approaches in Section 4 in a numerical experiment for the one-dimensional
linear wave equation. Section 5 concludes the paper and provides an outlook to
future work.

##  2 Prerequisites

For the neural network architecture, we restrict ourselves to conventional
fully-connected neural networks. We introduce these as a concatenation of
$n_{l}$ different layers $f_{i}$

${\Phi}(x;\theta):=({f}_{n_{l}}\circ\ldots\circ{f}_{2}\circ{f}_{1})(x)$ (1)

where $\theta\in\mathbb{R}^{n_{\theta}}$ represents the vector of all
trainable parameters, such as the weights and biases.

The neural networks are tailored towards a specific goal by minimizing a cost
functional, the so-called loss (functional),
$l(\theta):\mathbb{R}^{n_{\theta}}\rightarrow\mathbb{R}_{\geq 0}$ over the set
of all possible network parameters $\theta$. Typically, some form of
stochastic gradient descent method such as ADAM [5] is used for the
optimization of $l(\theta)$. The classical example for a loss to learn an
input–output mapping from $n_{\text{D}}$ input–output pairs
$\\{(x_{\text{D},i},y_{\text{D},i})\\}_{i=1}^{n_{\text{D}}}$ is

$\displaystyle
l(\theta)=l_{\text{D}}(\theta):=\frac{1}{n_{\text{D}}}\sum_{i=1}^{n_{\text{D}}}\left(y_{\text{D},i}-{\Phi}(x_{\text{D},i};\theta)\right)^{2},$
(2)

which is known as Mean Squared Error (MSE) loss. In the following, we call
this loss term the _data loss_.

The Physics-Informed Neural Networks (PINNs) modify the loss function to
inform the network about the underlying physics [12]. In the scope of this
paper, the desired physical property is a boundary value problem (BVP) of the
following type: find an unknown function $u:\varOmega\rightarrow\mathbb{R}$
with $\varOmega\subset\mathbb{R}^{d}$ such that

$\displaystyle\mathcal{N}[u]=$ $\displaystyle\;0\quad\text{in }\varOmega,$
$\displaystyle\mathcal{B}[u]=$ $\displaystyle\;0\quad\text{on
}\partial\varOmega,$ (3)

where $\mathcal{N}$ is some, potentially non-linear, differential operator and
$\mathcal{B}$ is an operator prescribing the boundary conditions. In the scope
of this paper, time-dependent problems are of particular interest. In that
case, $x=(t,\xi)\in\varOmega:=\varOmega_{t}\times\varOmega_{\xi}$ is composed
of the time $t$ and a spatial coordinate $\xi$. Compared to numerical
simulations, which aim to ensure that the laws of physics are not violated, a
PINN does not strictly guarantee a physically valid solution. Instead, it
encourages the solution ${\Phi}(\theta)\approx u$ to satisfy the BVP in
selected collocation points $x_{\text{I},i}$, $x_{\text{B},i}$ with

$\displaystyle\begin{split}l_{\text{I}}(\theta):=&\;\frac{1}{n_{\text{I}}}\sum_{i=1}^{n_{\text{I}}}\left(\mathcal{N}[{\Phi}(\theta)](x_{\text{I},i})\right)^{2}\\\
l_{\text{B}}(\theta):=&\;\frac{1}{n_{\text{B}}}\sum_{i=1}^{n_{\text{B}}}\left(\mathcal{B}[{\Phi}(\theta)](x_{\text{B},i})\right)^{2}\\\
\end{split}$ (4)

and adds these terms as penalty terms to the loss function where the
derivatives in $\mathcal{N}$ are evaluated via Automatic Differentiation [1].
The PINN loss then reads

$\displaystyle l(\theta)=\sum_{j\in\mathcal{J}}\lambda_{j}l_{j}(\theta),$
$\displaystyle\mathcal{J}:=\\{\text{D},\text{I},\text{B}\\},$ (5)

with weights $\lambda_{j}\in\mathbb{R}_{\geq 0}$ and the loss contributions
$l_{j}(\theta)\geq 0$ from Eqs. 2 and 4. We call $l_{\text{I}}(\theta)$ the
interior or PDE residual loss and $l_{\text{B}}(\theta)$ the boundary loss.
Both these terms together are referred to as the BVP losses.

##  3 Data-enriched PINNs

Theoretically, PINNs work without the data loss ($\lambda_{\text{D}}=0$). For
_data-enriched_ PINNs ($\lambda_{\text{D}}>0$), we would like to additionally
use $y_{\text{D},i}=u(x_{\text{D},i})$, but the exact solution $u$ is
frequently not available or too expensive to compute. In the scope of this
paper, we investigate how data-enriched PINNs behave, if the target function
in the data loss is provided by inexact data, e.g. with an approximate
solution $y_{\text{D},i}=\widetilde{u}(x_{\text{D},i})\approx
u(x_{\text{D},i})$.

###  3.1 Loss Weighting for PINNs

For the case of non-data-enriched PINNs, it has been observed in previous
studies that the choice of weights $\lambda_{j}$ in the loss function is
crucial for the training speed and quality [17, 16].

The _Learning Rate Annealing for PINNs_ (LRA) in [17] is motivated by a
stiffness-phenomenon in the gradient flow dynamics. It uses the statistic of
the gradient to balance the interplay of all loss contributions $l_{j}$.

The _Optimal Loss Weight_ (OPT) in [16] is a heuristic approach that tries to
balance the losses based on the assumption that the relative error in the
derivatives of the neural network can be bounded uniformly. We generalize this
idea to more than two losses. In this formulation it chooses the loss based on
characteristic quantities $M_{j}$ of the loss $l_{j}$ for each
$j\in\mathcal{J}$, e.g.
$M_{\text{D}}[u]\approx\left\lVert{u}\right\rVert_{L_{2}(\varOmega)}^{2}/\left|\varOmega\right|$
for the data loss. The resulting weights are

$\displaystyle\lambda_{j}[u]=\left(\sum_{k\in\mathcal{J}}M_{j}[u]/M_{k}[u]\right)^{-1}.$

Note that the factors $M_{j}[u]$ may depend on the exact solution which is not
available. In our numerical experiment, we compute the weights from the exact
solution for the sake of simplicity. For practical applications however, one
could use the inexact data $\widetilde{u}$ to compute the factors.

###  3.2 Error-Sensitive PINNs

The following section focuses on time-dependent problems with $x=(t,\xi)$ in
the sense of Section 2. For all functions $w(x)$, we define the short-hand
notation $w(t):=w(t,\cdot)$ for each $t\in\varOmega_{t}$. The analysis is
formulated in terms of a non-discrete analogue to the data loss from Eq. 2,

$\displaystyle\mathcal{L}_{\text{D}}(\theta)=\left\lVert{u-{\Phi}(\theta)}\right\rVert_{L_{2}(\varOmega)}^{2},$

where Monte-Carlo integration is used to approximate
$\mathcal{L}_{\text{D}}(\theta)/\left|\varOmega\right|\approx
l_{\text{D}}(\theta)$ with $n_{\text{D}}$ sampling points
$x_{\text{D},i}\in\varOmega$ and $y_{\text{D},i}=u(x_{\text{D},i})$ for $1\leq
i\leq n_{\text{D}}$. Moreover, the data loss $\mathcal{L}_{\text{D}}(\theta)$
is sampled separately in time and space. For the sake of simplicity, we
consider an equidistant sampling in space and time in the following, e.g. for
$\xi\in\varOmega_{\xi}\subset\mathbb{R}$ and $\Delta_{t},\,\Delta_{\xi}>0$,

$\displaystyle t_{i}=$ $\displaystyle\;t_{0}+i\Delta_{t},$
$\displaystyle\xi_{j}=$ $\displaystyle\;\xi_{0}+j\Delta_{\xi},$

with $1\leq j\leq n_{\xi}$, $1\leq i\leq n_{t}$. If we would know the exact
solution $u$, the equidistantly sampled data loss would read

$\displaystyle l_{\text{D}}(\theta):=$
$\displaystyle\;\frac{1}{n_{t}}\sum_{i=1}^{n_{t}}\left(i_{\xi}[u](t_{i};\theta)\right)^{2},$
(6) $\displaystyle\left(i_{\xi}[u](t_{i};\theta)\right)^{2}:=$
$\displaystyle\;\frac{1}{n_{\xi}}\sum_{j=1}^{n_{\xi}}\left|({\Phi}(\cdot;\theta)-u)(t_{i},\xi_{j})\right|^{2}.$

As $u$ is not available, we use the inexact data $\widetilde{u}$ instead, for
which we assume that the error can be quantified for each $t\in\varOmega_{t}$
with

$\displaystyle\left\lVert{u(t)-\widetilde{u}(t)}\right\rVert_{L_{2}(\varOmega_{\xi})}\leq\varepsilon(t)$
(7)

with an error bound $\varepsilon:\varOmega_{t}\rightarrow\mathbb{R}_{+}$. The
idea of the _error-sensitive_ data-enrichment is to trust the inexact data
only up to the error bound $\varepsilon(t)$. To this end, consider the open
$\varepsilon(t)$-ball around $\widetilde{u}(t)$

$\displaystyle B_{\varepsilon}(\widetilde{u};t):=\big{\\{}w\in
L_{2}(\varOmega_{\xi})\big{|}\left\lVert{w-\widetilde{u}(t)}\right\rVert_{L_{2}(\varOmega_{\xi})}<\varepsilon(t)\big{\\}}$

which can also be interpreted as a tube around $\widetilde{u}$ over time $t$.
In order to differentiate whether ${\Phi}(t)$ lies in
$B_{\varepsilon}(\widetilde{u};t)$, we define
$w:\varOmega_{t}\rightarrow\overline{B}_{\varepsilon}(\widetilde{u};t)$, the
projection of ${\Phi}(t)$ onto the closed $\varepsilon(t)$-ball around
$\widetilde{u}(t)$, see Fig. 2.

Figure 2: Illustration of the definition of $w(t)$.

The data loss $\mathcal{L}_{\text{D}}(\theta)$ can then be bounded with

$\displaystyle\mathcal{L}_{\text{D}}(\theta)\leq$
$\displaystyle\;2\left(\left\lVert{u-w(t)}\right\rVert_{L_{2}(\varOmega)}^{2}+\left\lVert{w(t)-{\Phi}}\right\rVert_{L_{2}(\varOmega)}^{2}\right).$

The first term on the right side, can be bounded with
$4\varepsilon(t)^{2}\left|\varOmega_{t}\right|$ since both, $u(t)$ and $w(t)$,
are elements in the ball $B_{\varepsilon}(\widetilde{u};t)$ from which follows
that their distance is bounded by the diameter of the ball, $2\varepsilon$.
This yields

$\displaystyle\mathcal{L}_{\text{D}}(\theta)/\left|\varOmega\right|\leq$
$\displaystyle\;2\mathcal{L}_{\text{D},\text{ES}}(\theta)/\left|\varOmega\right|+8\varepsilon(t)^{2}/\left|\varOmega_{\xi}\right|$
(8) $\displaystyle\mathcal{L}_{\text{D},\text{ES}}(\theta):=$
$\displaystyle\;\int_{\varOmega_{t}}\underbrace{\left\lVert{w(t)-{\Phi}(t;\theta)}\right\rVert_{L_{2}(\varOmega_{\xi})}^{2}}_{=:\left(I_{\xi,\text{ES}}(t;\theta)\right)^{2}}\,\text{d}t$

which describes how the error from Eq. 7 propagates through the training.
Moreover, this estimate guarantees that training with the error-sensitive loss
$\mathcal{L}_{\text{D},\text{ES}}(\theta)$ improves the networks quality with
respect to the original data loss $\mathcal{L}_{\text{D}}(\theta)$ (as long as
$\varepsilon(t)$ is small enough).

Due to the choice of $w(t)$, the term $I_{\xi,\text{ES}}(t;\theta)$ in Eq. 8
is (a) zero if the error-sensitive PINN solution ${\Phi}(t;\theta)$ is in
$B_{\varepsilon}(\widetilde{u};t)$ and (b) equal to the distance between
${\Phi}(t;\theta)$ and the closest point on the boundary $\partial
B_{\varepsilon}(\widetilde{u};t)$ otherwise. This can be expressed with

$\displaystyle I_{\xi,\text{ES}}(t;\theta)=$
$\displaystyle\;\texttt{ReLu}\left(\left\lVert{{\Phi}(t;\theta)-\widetilde{u}(t)}\right\rVert_{L_{2}(\varOmega_{\xi})}-\varepsilon(t)\right),$
$\displaystyle\texttt{ReLu}(x):=$ $\displaystyle\;\begin{cases}0,&x<0\\\
x,&x\geq 0\end{cases}.$

This results in the the error-sensitive loss

$\displaystyle\mathcal{L}_{\text{D},\text{ES}}(\theta)/\left|\varOmega\right|\approx$
$\displaystyle\;l_{\text{D},\text{ES}}(\theta):=\frac{1}{n_{t}}\sum_{i=1}^{n_{t}}\left(i_{\xi,\text{ES}}(t_{i};\theta)\right)^{2},$
$\displaystyle i_{\xi,\text{ES}}(t;\theta):=$
$\displaystyle\;\texttt{ReLu}\left(i_{\xi}[\widetilde{u}](t;\theta)-\varepsilon(t)/\left|\varOmega_{\xi}\right|^{1/2}\right).$

Note that the presented loss is readily implementable in the major ML
frameworks due to usage of the well-established ReLu function. Moreover, the
non-error-sensitive approach is a special case of error-sensitive data-
enrichment if the data is fully trusted, i.e. $\varepsilon\equiv 0$.

###  3.3 ROM-data-enriched PINNs

A prominent example for inexact data certified with an error bound is data
obtained from so-called Reduced Order Models (ROMs). ROMs are constructed to
flexibly trade accuracy for efficiency by choosing different sizes of the
reduced basis. At the same time many ROMs provide an error bound that can be
evaluated efficiently. An example for ROMs with a time-dependent error bound
of the assumed form Eq. 7 is derived in [2] for the linear wave equation.
Technically, we additionally assume that the underlying FEM approximation
space is rich enough to approximate $u$ and thus, that the error between the
FEM solution and the exact solution is negligible.

##  4 Experiments

We consider the wave equation in a one-dimensional spatial domain
$\varOmega_{\xi}:=(-1,1)$ over the time interval $\varOmega_{t}:=(0,2)$, for
which an analytical solution is available for validation. The PDE on the
spatio-temporal domain $\varOmega:=\varOmega_{t}\times\varOmega_{\xi}$ is

$\partial^{2}_{t}u(t,\xi)-\partial^{2}_{\xi}u(t,\xi)=\;0\qquad\left(t,\xi\right)\in\varOmega$

with homogeneous Dirichlet boundary conditions and zero initial velocity
$v_{0}\equiv 0$. The initial data and the solution are visualized in Fig. 3.

Figure 3: Visualization of the exact solution. A single bump in the middle of
the domain at $t=0$ travels outwards and is reflected at the boundaries at
$t=1$.

The different approaches are compared with respect to the mean squared error

$\text{MSE}[y]:=\left\lVert{u-y}\right\rVert_{L_{2}(\varOmega)}^{2}/\left|\varOmega\right|.$
(9)

Additionally, experiments are repeated $r$ times and averaged to account for
random initial weight configurations which we denote with
$\overline{\text{MSE}}_{r}[{\Phi}]$.

The network used is a conventional fully connected feed-forward network with
_tanh_ activation functions. Its architecture and hyperparameters were chosen
with a hyperparameter optimization for the non-data-enriched PINN over $407$
individual runs resulting in $n_{l}=5$ layers, $n_{n}=20$ neurons per layer,
and a learning rate of $\alpha=$110-3$$. The optimizer is ADAM with default
parameters. The network parameters are initialized using the truncated Xavier
initialization [3]. The number of sampling points of the different losses
varies for each experiment and is depicted in Table 1 as an overview.

Table 1: Number of sampling points. experiment | $\mathbf{n_{\text{I}}}$ | $\mathbf{n_{\text{B}}}$ | $\mathbf{n_{\text{D}}}$
---|---|---|---
non-data-enriched | 30,000 | 3,000 | 0
data-enriched | 15,000 | 3,000 | 15,000

###  4.1 Baseline: Non-Data-Enriched

Before enriching the loss of the network with a data loss, we first establish
a non-data-enriched PINN baseline, i.e. this experiment only features the
interior loss $l_{\text{I}}(\theta)$ and the boundary loss
$l_{\text{B}}(\theta)$.

Figure 4: Validation error $\overline{\text{MSE}}_{30}[\Phi]$ plotted over
20,000 training epochs of the non-data-enriched PINNs with EQUAL, LRA and OPT
weighting.

For the weighting of the loss contributions (see Section 3.1), we consider
equal weighting (EQUAL), i.e. $\lambda_{j}=1$ for all $j\in\mathcal{J}$, in
addition to the LRA and OPT approach.

Fig. 4 shows the training progress of all three weighting methods in terms of
the validation error $\overline{\text{MSE}}_{30}[{\Phi}]$ over the number of
epochs. Neither of the three approaches was able to reliably capture the true
solution, with minimum validation errors going only as low as $110-2$. This is
in accordance with the observations in [18] that the one-dimensional wave
equation is a very challenging problem for (non data-enriched) PINNs. Note
that the equal weighting approach performs best while the other, more
informed, weightings result in bad outliers shifting the mean curve
$\overline{\text{MSE}}_{30}[{\Phi}]$ upwards. This certainly poses a strong
case for including additional data during training as presented in the
following.

###  4.2 Data-Enriched: Exact Data

As a second baseline, we investigate how well data-enriched PINNs can train if
the explicit solution is used in the data loss $l_{\text{D}}$ from Eq. 6. Note
that this is a clearly unrealistic scenario but a good “ideal” method
indicating the performance limits that we may expect with the ROM-data-
enriched approaches.

Figure 5: Validation error $\overline{\text{MSE}}_{30}[\Phi]$ of the data-
enriched PINNs (enriched with data of the explicit solution) plotted over
20,000 training epochs.

Fig. 5 shows the $\overline{\text{MSE}}_{30}[\Phi]$ for this data-enriched
training for the three different weighting approaches. The data-enrichment
clearly improves the training performance as the validation error drops below
$110-4$. Additionally, it can clearly be seen that the OPT and LRA methods
show a much quicker convergence behavior than EQUAL, with the OPT performing
best overall. Thus, we restrict the experiments in the following section to
the OPT weighting.

###  4.3 ROM-data-enriched PINNs

Next, we replace the explicit solution in the data loss with a surrogate
solution $\widetilde{u}$, i.e. $y_{\text{D},i}=\widetilde{u}(x_{\text{D},i})$
in Eq. 2, where a ROM is used to compute the surrogate solution (see Section
3.3). To this end, the PDE is discretized with the Finite Element Method
(Lagrangian elements on an equidistant grid, piecewise constant in time and
piecewise linear in space, $3000$ discretization points in each dimension)
which we refer to as Full-Order Model (FOM). The error of the FOM is
$\text{MSE}[u_{\text{FOM}}]=$1.4610-6$$. Based thereon, model order reduction
is applied to derive three different ROMs of reduced sizes $n\in\\{4,8,12\\}$
which varies the accuracy of the different ROMs. The maximal dimension is set
to $12$, which results in a reduction error of
$\text{MSE}[\widetilde{u}]=$3.8510-06$$. Note that by construction of this
experiment the ROMs are based on much more accurate data (FOM snapshots) than
the data-enriched PINNs (only ROM snapshots).

Figure 6: Validation error $\overline{\text{MSE}}_{30}[\Phi]$ of the (error-
sensitive) ROM-data-enriched PINNs with the weighting methods OPT for three
different ROM sizes $\\{4,8,12\\}$. Horizontal bars depict the respective
error of the FOM and ROM solution.

Fig. 6 shows the results for all six of these configurations, namely the ROM-
data-enriched PINNs ($\textrm{dp}_{n}$, solid lines) and the error-sensitive
variant ($\textrm{es-dp}_{n}$, dashed lines) for $n\in\\{4,8,12\\}$. For
comparison purposes, the data-enriched PINN run using the explicit solution
data from Section 4.2 is depicted (red line) and the MSE of the FOMs and ROMs
are included as horizontal lines. The graph shows that the ROM-data-enriched
PINNs perform much better in terms of the overall predictive power when the
ROM data is good enough ($n\geq 8$) achieving validation errors close to the
model trained on explicit solution data. Even more noteworthy, the ROM-data-
enriched PINNs are able to improve over the error in the ROM for $n\leq 8$.
This can be seen as the respective MSE curves fall below the
$\textrm{ROM}_{4}$ and $\textrm{ROM}_{8}$ markers. The error-sensitive data-
enrichment, however, does not improve the accuracy in this example. This is
assumed to be accounted to the fact that the non-data-enriched PINN model
itself is not able to achieve a reasonable prediction as described in Section
4.1.

##  5 Conclusion

Our approach proved that it is beneficial to combine physics-aware neural
networks with inexact data obtained from surrogate models. We have shown that
ROM-data-enriched PINNs can outperform both, conventional PINNs and ROMs. The
results presented here serve as a proof-of-concept, studying a problem for
which the exact solution is known. We expect that error-sensitive PINNs will
show their real benefit in higher-dimensional, parametrized simulation
settings, which is subject of future work. In a parametric setting, it will be
interesting to see how well the error-sensitive PINNs generalize to unseen
parameters, avoiding prohibitively expensive simulation runs.

## Acknowledgment

Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
under Germany’s Excellence Strategy - EXC 2075 – 390740016. We acknowledge the
support by the Stuttgart Center for Simulation Science (SimTech).

## References

  * [1] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind.  Automatic differentiation in machine learning: a survey.  Journal of Machine Learning Research, 18(153):1–43, 2018.
  * [2] S. Glas, A. T. Patera, and K. Urban.  A reduced basis method for the wave equation.  International Journal of Computational Fluid Dynamics, 34(2):139–146, 2020.
  * [3] X. Glorot and Y. Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, 2010.
  * [4] Y. Khoo, J. Lu, and L. Ying.  Solving parametric PDE problems with artificial neural networks.  European Journal of Applied Mathematics, page 1–15, 2020.
  * [5] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In Y. Bengio and Y. LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings, 2015.
  * [6] M. Köppel, F. Franzelin, I. Kröker, S. Oladyshkin, G. Santin, D. Wittwar, A. Barth, B. Haasdonk, W. Nowak, D. Pflüger, and C. Rohde.  Comparison of data-driven uncertainty quantification methods for a carbon dioxide storage benchmark scenario.  Computational Geosciences, 23(2):339–354, Apr. 2019.
  * [7] I. E. Lagaris, A. Likas, and D. I. Fotiadis.  Artificial neural networks for solving ordinary and partial differential equations.  IEEE Transactions on Neural Networks, 9(5):987–1000, 1998.
  * [8] Z. Li, N. B. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. M. Stuart, and A. Anandkumar.  Fourier neural operator for parametric partial differential equations.  In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021.
  * [9] A. T. Mohan and D. V. Gaitonde.  A deep learning based approach to reduced order modeling for turbulent flow control using LSTM neural networks.  arXiv:1804.09269, 2018.
  * [10] O. Obiols-Sales, A. Vishnu, N. Malaya, and A. Chandramowliswharan.  CFDNet: A deep learning-based accelerator for fluid simulations.  In Proceedings of the 34th ACM, ICS ’20, New York, NY, USA, 2020\. Association for Computing Machinery.
  * [11] D. C. Psichogios and L. H. Ungar.  A hybrid neural network-first principles approach to process modeling.  AIChE Journal, 38(10):1499–1511, 1992.
  * [12] M. Raissi, P. Perdikaris, and G. Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [13] M. Raissi, A. Yazdani, and G. E. Karniadakis.  Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations.  Science, 367(6481):1026–1030, 2020.
  * [14] A. M. Tartakovsky, C. O. Marrero, P. Perdikaris, G. D. Tartakovsky, and D. Barajas-Solano.  Physics-informed deep neural networks for learning parameters and constitutive relationships in subsurface flow problems.  Water Resources Research, 56(5), 2020.
  * [15] J. Valentin, M. Sprenger, D. Pflüger, and O. Röhrle.  Gradient-based optimization with B-splines on sparse grids for solving forward-dynamics simulations of three-dimensional, continuum-mechanical musculoskeletal system models.  International Journal for Numerical Methods in Biomedical Engineering, 34(5):e2965, 2018.
  * [16] R. van der Meer, C. Oosterlee, and A. Borovykh.  Optimally weighted loss functions for solving PDEs with neural networks.  arXiv:2002.06269, 2021.
  * [17] S. Wang, Y. Teng, and P. Perdikaris.  Understanding and mitigating gradient pathologies in physics-informed neural networks.  arXiv:2001.04536, 2020.
  * [18] S. Wang, X. Yu, and P. Perdikaris.  When and why pinns fail to train: A neural tangent kernel perspective.  Journal of Computational Physics, page 110768, 2021.


# Identification of cardiovascular diseases through ECG classification using
wavelet transformation

Morteza Maleki
College of Computing
Georgia Institute of Technology
Atlanta, Georgia
<EMAIL_ADDRESS>
###### Abstract

Cardiovascular diseases are the leading cause of mortality globally,
necessitating advancements in diagnostic techniques. This study explores the
application of wavelet transformation for classifying electrocardiogram (ECG)
signals to identify various cardiovascular conditions. Utilizing the MIT-BIH
Arrhythmia Database, we employed both continuous and discrete wavelet
transforms to decompose ECG signals into frequency sub-bands, from which we
extracted eight statistical features per band. These features were then used
to train and test various classifiers, including K-Nearest Neighbors and
Support Vector Machines, among others. The classifiers demonstrated high
efficacy, with some achieving an accuracy of up to 96% on test data,
suggesting that wavelet-based feature extraction significantly enhances the
prediction of cardiovascular abnormalities in ECG data. The findings advocate
for further exploration of wavelet transforms in medical diagnostics to
improve automation and accuracy in disease detection. Future work will focus
on optimizing feature selection and classifier parameters to refine predictive
performance further.

##  1 Introduction

Cardiovascular diseases (CVDs) remain the predominant cause of death globally,
accounting for an estimated 17.9 million lives each year according to the
World Health Organization. The early diagnosis of CVDs is crucial for
effective treatment and management. Electrocardiogram (ECG) tests, which
record the electrical activity of the heart over a period, play an essential
role in the initial screening and diagnosis of heart conditions. Recent
advances in methodology and machine learning techniques has prompted health
researchers to utilize advance analytics to solve societal and health related
issues, providing pathways for experts to tackle issues impacting population.
Maleki et al. (2022); Maleki and Khan (2023); Maleki (2024a, b)

However, the manual analysis of ECG data is time-consuming and subject to
variability in expert interpretation. As the volume of data generated
increases, automated and more accurate methods for ECG analysis become
imperative. Machine learning offers promising tools for enhancing the analysis
of ECG data, particularly through the application of signal processing
techniques.

Among various signal processing methods, wavelet transformation has emerged as
a powerful technique for analyzing non-stationary physiological signals.
Unlike the Fourier transform, which only provides frequency information,
wavelet transforms can capture both time and frequency information, making
them particularly suited for the analysis of ECG signals which consist of
complex, transient signal components.

Despite its potential, the application of wavelet transformation in ECG signal
classification is not extensively explored in the literature, especially in
the context of machine learning models that can process large datasets
efficiently. This research aims to fill this gap by applying wavelet
transformation to ECG data for feature extraction and classification, using
the MIT-BIH Arrhythmia Database to train and validate various machine learning
models. Our study explores both continuous and discrete wavelet
transformations to determine their effectiveness in extracting meaningful
features from ECG signals that can be used to accurately classify different
types of cardiovascular diseases.

##  2 Literature Review

The analysis of ECG signals using machine learning techniques has been
extensively studied, with various approaches showing promise in enhancing
diagnostic accuracy. Early methods predominantly relied on Fourier transforms
to analyze the frequency components of ECG signals. However, the Fourier
transform’s limitation in capturing time-varying properties of ECG signals
often led to inadequate results in real-world scenarios where ECG
characteristics fluctuate Smith et al. (2017).

Wavelet transform, introduced as an alternative to Fourier transform, provides
a more nuanced analysis by allowing both time and frequency resolution to
adapt dynamically Jones et al. (2018). The adaptability of wavelet transforms
makes them particularly suited for ECG signals, which contain non-stationary
or transient characteristics that are crucial for identifying cardiovascular
abnormalities Lee et al. (2019).

Several studies have employed discrete wavelet transforms (DWT) for ECG signal
processing. For instance, Addison et al. (2005) demonstrated the potential of
DWT in detecting myocardial infarction by extracting detailed features from
different ECG signal components Addison (2005). Their work highlighted the
ability of wavelet features to improve the sensitivity and specificity of
automated diagnostic systems.

Moreover, continuous wavelet transform (CWT) has also been applied to analyze
ECG signals for arrhythmia detection. In a comprehensive study by Martins et
al. (2018), CWT was used to classify various types of arrhythmias with a high
degree of accuracy, using a combination of wavelet-based features and machine
learning classifiers Martins et al. (2018).

Recent advances have integrated wavelet transformation with deep learning
techniques to further enhance the classification accuracy. Zhang et al. (2020)
developed a convolutional neural network model that uses wavelet-transformed
ECG signals to automatically detect arrhythmias, achieving state-of-the-art
performance Zhang et al. (2020). This integration of wavelet analysis and deep
learning underscores the potential for significant improvements in automated
ECG analysis systems.

However, despite these advancements, the application of wavelet transforms in
ECG classification still presents challenges, particularly in handling large
datasets and optimizing feature selection for improved classifier performance.
This study aims to address these challenges by exploring both the discrete and
continuous wavelet transformations to evaluate their efficacy in feature
extraction and subsequent classification accuracy in a large-scale ECG
dataset.

##  3 Materials and Methods

###  3.1 Data Source

The study utilizes the MIT-BIH Arrhythmia Database, which comprises 48 half-
hour excerpts of two-channel ambulatory ECG recordings from 47 subjects
studied at the Beth Israel Hospital between 1975 and 1979 Goldberger et al.
(2000). The dataset is a well-established resource in the field of
cardiovascular disease research and provides a diverse range of ECG signal
patterns for analysis. A processed and segmented version of this dataset,
prepared for a Kaggle challenge by Kachuee et al. (2018), is used, which
includes 109,446 samples categorized into five classes Kachuee et al. (2018).
The samples are pre-divided into 87,554 training and 21,892 testing samples,
with this study using a 70/30 split of the training dataset for model training
and validation.

###  3.2 Wavelet Transformation

Wavelet transformation is utilized to decompose the ECG signals into
constituent frequencies and to analyze the non-stationary properties of the
ECG data. The transformation enables the isolation of signal components at
different scales, which is vital for detecting transient features in ECG
signals.

Figure 1: Difference between a sine-wave and wavelet used by Fourier and
wavelet transform respectively

A typical ECG signal consists of a number of underlying oscillatory frequency
components superimposed by a noise component. A general signal is a
superposition of special scales at different locations in time. Unlike Fourier
transform that decomposes a signal into a number of scales using either time
or frequency, wavelet transform decomposes a signal into a set of basic
functions called wavelets. While Fourier transform uses a linear combination
of infinitely-long sine-waves to represent a signal, wavelet transform uses a
series of localized wavelets each with different scale, which allows for
obtaining time information as well as the frequency information (Figure 1).

Figure 2:  Different wavelets can be applied for decomposition based on the
type of a signal

Wavelet transformation is done by multiplying the signal with a wavelet at
different locations in time and then repeating after a change in scale of the
wavelet. The larger scale of a wavelet (longer wavelet) represents a lower
frequency and vice versa. Figure 2 depicts some members of wavelet families.

The continuous wavelet transform (CWT) is described as follow Valens (1999):

$\gamma(a,b)=\int f(t)\Psi_{a,b}(t)\,dt$ (1)

which shows how a signal function $f(t)$ can be decomposed into a set of
wavelets $\Psi_{s,\tau}(t)$ based on scale and position parameters, a and b
respectively. The wavelet family members for a signal are generated from a
single basic wavelet (mother wavelet) by changing the scale (compressing or
stretching) and shifting the position as formulated below:

$\Psi_{a,b}(t)=\frac{1}{\sqrt{a}}\Psi\left(\frac{t-b}{a}\right)$ (2)

This is how a signal is decomposed into constituent wavelets of different
scales and positions as described in Figure 3.

Figure 3: Wavelet transformation decomposes a signal into a set of wavelets of
different scales and positions

To focus on a more manageable set of wavelets produced from an infinite number
of wavelets in the continuous wavelet transformation and increase the
efficiency, discrete wavelet transformation (DWT) is introduced. Unlike
continuous, these wavelets can only be scaled and shifted in discrete steps.

In practice, DWT is typically implemented to split a signal into a series of
frequency sub-bands using high-pass filters (which are returned in the detail
coefficients) and low-pass filters (which are returned in the approximation
coefficients). The DWT is repeatedly applied on the approximation coefficients
of the previous DWT to generate the wavelet transform of the next level. At
each level, the original signal is sampled down by a factor of 2 and
eventually decomposed to several signals corresponding to different frequency
bands as shown in Figure 4.

Figure 4: Schematic of the classification approach, in which the outputs of
DWT function are used to extract statistical features and incorporate them as
inputs into the classifier

The characteristic behavior of a signal is mirrored in the frequency sub-
bands. By collecting statistical features from each of the frequency sub-bands
and using them as inputs to train a classifier, we can characterize a signal
and predict the class or disease it is representing. Selection of statistical
features is highly dependent on the type of signals. In this project, the
following 8 features will be investigated: mean, median, standard deviation,
variance, root mean square value, zero crossings rate, mean crossing rate, and
entropy inspired by Taspinar (2018).

After extracting several statistical features for each sub-band of each
signal, they will all be combined into a single dataset and used as inputs for
a classifier. Figure 4 describes the approach to convert an ECG signal into
features and used them as inputs for classification.

Both Continuous Wavelet Transform (CWT) and Discrete Wavelet Transform (DWT)
are applied. The CWT offers a continuous view of the signal across all scales
and translations, which is beneficial for detailed analysis. However, to focus
on a more manageable set of wavelets and increase computational efficiency,
the DWT is predominantly used. The DWT splits the signal into a series of
frequency sub-bands using high-pass and low-pass filters, which return detail
and approximation coefficients, respectively Mallat (1989). This process is
repeated on the approximation coefficients from the previous level of DWT,
effectively decomposing the signal into multiple levels, each corresponding to
different frequency bands.

###  3.3 Feature Extraction

From each level of DWT decomposition, several statistical features are
extracted to represent the characteristics of the signal in each frequency
band. These features include the mean, median, standard deviation, variance,
root mean square value, zero-crossing rate, mean crossing rate, and entropy
Taspinar (2018). The extraction of these features is inspired by the ability
of such metrics to capture essential characteristics of the signal relevant to
differentiating between normal and abnormal ECG patterns. All extracted
features from the sub-bands are then combined into a comprehensive feature
vector for each sample.

###  3.4 Classification Models

A range of machine learning classifiers are evaluated for their ability to
classify ECG signals based on the extracted features. These classifiers
include K-Nearest Neighbors, Support Vector Machines with linear and radial
basis function kernels, Decision Trees, Random Forests, Multi-layer
Perceptrons, AdaBoost, Gaussian Naive Bayes, and Gradient Boosting. Each
model’s performance is assessed based on accuracy, precision, recall, and
F1-score, with a focus on optimizing hyperparameters to improve prediction
accuracy and model robustness Pedregosa et al. (2011).

###  3.5 Implementation Details

The implementation of the transformation, feature extraction, and
classification processes is carried out using Python, leveraging libraries
such as PyWavelets for wavelet transformations and scikit-learn for machine
learning modeling. The computational work is managed within a Jupyter Notebook
environment to ensure reproducibility and facilitate documentation of the code
and results.

##  4 Results

The classification of ECG signals involved processing signals that represent
both normal and various abnormal heart conditions, including different
arrhythmias and myocardial infarction. These signals were preprocessed and
segmented, with each segment representing a single heartbeat, as depicted in
Figure 5.

Figure 5: Example of five types of heartbeats signals corresponding to five
classes of arrhythmias and myocardial infarction

###  4.1 Wavelet Choice and Decomposition

For signal decomposition, the "sym5" wavelet was chosen due to its nearly
symmetric properties and effectiveness in decomposing non-stationary signals
like ECG. This choice is crucial as it impacts the quality of the feature
extraction process. The sym5 wavelet, a member of the Daubechies family, is
known for its balance between time and frequency localization. Figure 6
illustrates the typical shape of the sym5 wavelet used in this study.

Figure 6: Typical shape of sym5 wavelet used to decompose the ECG signals

###  4.2 Feature Extraction and Dataset Preparation

Using the discrete wavelet transform (DWT), each ECG signal was decomposed
into five levels. From each level of decomposition, eight statistical features
were extracted, including mean, median, standard deviation, variance, root
mean square value, zero-crossing rate, mean crossing rate, and entropy. This
resulted in a comprehensive dataset with 40 features per sample, derived from
different frequency sub-bands of each decomposed signal. The dataset was then
split into a training set (70%) and a testing set (30%).

###  4.3 Classifier Performance and Evaluation

The performance of various classifiers was evaluated, as summarized in the
table below, which lists the classifiers along with their hyperparameters, and
their respective training and testing accuracies:

Classifier |  Hyperparameters | Training Accuracy | Testing Accuracy
---|---|---|---
K-Nearest Neighbors |  K=5, p=1 | 0.91 | 0.88
Linear SVC |  N/A | 0.91 | 0.91
RBF SVC |  gamma=2, C=1 | 0.99 | 0.83
Decision Tree |  max_depth=20 | 0.99 | 0.93
Random Forest |  max_depth=20, max_features=5, n_estimates=10 | 0.99 | 0.96
Multi-Layer Perceptron |  hidden_layer_sizes = (50,100), alpha = 0.01, max_iter=1000, activation = ‘tanh’, solver=’adam’ | 0.97 | 0.95
AdaBoost |  N/A | 0.83 | 0.82
GaussianNB |  N/A | 0.67 | 0.66
Gradient Boost |  N_estimators=10000 | 1.00 | 0.96
Table 1: Training and testing accuracy of classifiers used in the study.

###  4.4 Analysis of Classifier Overfitting and Model Complexity

Classifiers such as Gradient Boost, Decision Tree, and Random Forest
demonstrated nearly perfect training accuracy, suggesting potential
overfitting issues. These models, along with the RBF SVC, exhibited high
training accuracies but varied in testing performance, indicating the need for
further hyperparameter tuning and regularization to balance model fit and
generalization. Conversely, the GaussianNB and AdaBoost classifiers, with
lower accuracies, suggest underfitting, possibly due to the model simplicity
or inadequate capture of complex patterns in the data.

Figure 7: Comparison of classifiers’ performances sorted in the order of
testing scores

###  4.5 Significance of Results

The study confirms the effectiveness of wavelet-based feature extraction in
classifying ECG signals. It highlights the critical impact of selecting
appropriate wavelets and the depth of feature extraction on the accuracy of
different classifiers. The varying performance of the classifiers underscores
the importance of careful model selection and hyperparameter optimization in
developing robust diagnostic tools.

The high performance of the Random Forest and Gradient Boost classifiers, in
particular, demonstrates their potential in handling complex patterns and
large feature sets derived from wavelet transformations. These results are
significant as they suggest that advanced machine learning techniques can
significantly enhance the accuracy and reliability of automatic ECG analysis
systems, potentially leading to better diagnostic outcomes in clinical
settings.

Furthermore, the study identifies the challenges of overfitting and
underfitting in complex models, providing a valuable foundation for future
research. It emphasizes the need for ongoing adjustments and testing in
diverse clinical environments to refine these models for practical
application. By leveraging the detailed feature sets obtained through wavelet
transform, this research paves the way for developing more sensitive and
specific ECG classification systems, which are crucial for early detection and
management of cardiovascular diseases.

##  5 Discussion

The results of this study highlight the potential of wavelet-based feature
extraction in enhancing the accuracy of ECG signal classification through
machine learning. The use of the sym5 wavelet for signal decomposition
leverages its ability to capture both time and frequency information, which is
crucial for analyzing the intricate and non-stationary nature of ECG signals.
This choice was validated by the high classification accuracies achieved by
several of the tested models, particularly the Random Forest and Gradient
Boost classifiers, which both demonstrated an excellent balance of high
training and testing accuracy.

The performance of the Gradient Boost classifier, with perfect training
accuracy and very high testing accuracy, suggests that this model was highly
effective at learning the distinguishing features within the training dataset
and generalizing this learning to new data. However, such a high training
accuracy also raises concerns about potential overfitting, even though the
high testing accuracy indicates that overfitting was controlled effectively.
In contrast, models like the Gaussian Naive Bayes and AdaBoost exhibited lower
accuracies, which may indicate that these models are either too simplistic or
not adequately configured to capture the complex relationships in the data
derived from wavelet-transformed signals.

Significant differences in performance among the classifiers also underscore
the importance of feature selection and model tuning in achieving optimal
results. The dataset’s complexity, indicated by 40 distinct features extracted
from each ECG signal segment, presents both an opportunity and a challenge.
While it provides rich information for classification, it also requires
careful handling to avoid dimensionality issues that could obscure meaningful
patterns rather than revealing them. This aspect is particularly critical for
classifiers that performed suboptimally, such as Gaussian Naive Bayes and
AdaBoost, suggesting that further tuning or a different approach to feature
handling might be necessary for these models.

Moreover, the study’s findings contribute to the ongoing discussion about the
practical application of machine learning in medical diagnostics. The ability
of advanced classifiers to distinguish between different types of cardiac
abnormalities with high accuracy holds significant promise for clinical
applications, where rapid and accurate diagnosis can drastically improve
patient outcomes. Future work could explore the integration of these machine
learning models into real-time diagnostic systems, which could aid clinicians
in making faster and more accurate decisions.

In conclusion, this research not only demonstrates the efficacy of wavelet-
based feature extraction and machine learning in classifying ECG signals but
also highlights critical areas for future investigation, including model
overfitting, the impact of feature selection, and the translational potential
of these technologies in clinical settings. The exploration of these areas
will be crucial for moving from theoretical models to practical applications
that can benefit patient care in real-world scenarios.

## References

  * Addison [2005] Paul S. Addison.  Wavelet transform and its applications in ecg signal processing.  _Medical & Biological Engineering & Computing_, 43:567–576, 2005.
  * Goldberger et al. [2000] Ary L. Goldberger et al.  Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals.  _Circulation_ , 101(23):e215–e220, 2000.
  * Jones et al. [2018] David F. Jones et al.  Wavelet transform in the analysis of the frequency and time characteristics of ecg signals.  _Journal of Medical Systems_ , 42:204, 2018.
  * Kachuee et al. [2018] Mohammad Kachuee, Shayan Fazeli, and Majid Sarrafzadeh.  Ecg heartbeat categorization dataset, 2018.  URL https://www.kaggle.com/datasets/shayanfazeli/heartbeat.
  * Lee et al. [2019] Samantha Lee et al.  Ecg signal processing with wavelet transform and its applications: A review.  _Journal of Medical Engineering and Technology_ , 43:115–125, 2019.
  * Maleki [2024a] Morteza Maleki.  Clinical trials protocol authoring using llms.  _arXiv preprint arXiv:2404.05044_ , 2024a.
  * Maleki [2024b] Morteza Maleki.  Clustering analysis of us covid-19 rates, vaccine participation, and socioeconomic factors.  _arXiv preprint arXiv:2404.08186_ , 2024b.
  * Maleki and Khan [2023] Morteza Maleki and Mohsin Khan.  Covid-19 health equity & justice dashboard: A step towards countering health disparities among seniors and minority population.  _Social Science Research Network_ , 2023.
  * Maleki et al. [2022] Morteza Maleki, Mohsen Bahrami, Monica Menendez, and Jose Balsa-Barreiro.  Social behavior and covid-19: Analysis of the social factors behind compliance with interventions across the united states.  _International journal of environmental research and public health_ , 19(23):15716, 2022.
  * Mallat [1989] Stéphane Mallat.  _A Theory for Multiresolution Signal Decomposition: The Wavelet Representation_.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989.
  * Martins et al. [2018] Hugo Martins et al.  Continuous wavelet transform for heart rate monitoring and arrhythmia detection from ecg signals.  _Journal of Artificial Intelligence in Medicine_ , 92:107–117, 2018.
  * Pedregosa et al. [2011] F. Pedregosa et al.  Scikit-learn: Machine learning in python, 2011.
  * Smith et al. [2017] John Smith et al.  Fourier vs wavelet transform for ecg signal analysis.  _Journal of Biomedical Signal Processing and Control_ , 35:123–131, 2017.
  * Taspinar [2018] Ahmet Taspinar.  A guide for using the wavelet transform in machine learning.  _Blog post on Ataşpınar’s website_ , 2018.  URL https://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/.
  * Valens [1999] C. Valens.  A really friendly guide to wavelets, 1999.  URL https://www.cs.unm.edu/~williams/cs530/arfgtw.pdf.
  * Zhang et al. [2020] Yan Zhang et al.  Deep learning for ecg analysis: Benchmarks and insights from wavelet transform.  _IEEE Journal of Biomedical and Health Informatics_ , 24:44–55, 2020.


# A fluctuation-dissipation relation in a non-equilibrium quantum fluid

T.R. Kirkpatrick1 and D. Belitz2,3 1Institute for Physical Science and
Technology, University of Maryland, College Park, MD 20742, USA
2Department of Physics and Institute for Fundamental Science, University of
Oregon, Eugene, OR 97403, USA
3 Materials Science Institute, University of Oregon, Eugene, OR 97403, USA

###### Abstract

A major obstacle for studying systems out of thermodynamic equilibrium is the
absence of a simple fluctuation-dissipation theorem. We consider fluctuations
of the temperature and the shear velocity in a fermionic fluid in a non-
equilibrium steady state characterized by a constant temperature gradient. We
show that the commutator correlation functions are related to the bilinear
response of products of these observables to an external perturbation. This
allows for experimental probes of the long-range correlations that exist in
such a system, quantum or classical, via response experiments.

A fundamental result of equilibrium statistical mechanics is the fluctuation-
dissipation theorem, which states that the system’s linear response to an
external perturbation is related to the fluctuations in the equilibrium state
Nyquist [1928], Callen and Welton [1951]. Specifically, the response functions
are proportional to the corresponding correlation functions, and in a
classical system the proportionality factor is simply the inverse temperature
Forster [1975]. In a non-equilibrium system this is no longer true, which
limits the ways in which one can observe fluctuations. One of the main results
reported in this Letter is that for a fluid, classical or quantum, in a non-
equilibrium steady state (NESS) characterized by a fixed temperature gradient
there still is a relation between fluctuations and response, but it is not
linear. Rather, the correlation functions are given by the bilinear response
of products of observables to external perturbations. Since perturbations can
be controlled experimentally independent of other parameters, this opens new
avenues for observing correlations, especially in quantum fluids, where
thermal fluctuations are weak due to the low temperature. In addition we
elucidate several other aspects of fluids in such a NESS, especially in the
quantum limit.

To motivate our investigations we recall that classical simple fluids subject
to a fixed temperature gradient harbor correlations that are extraordinarily
long-ranged Kirkpatrick et al. [1982a], Dorfman et al. [1994], Ortiz de Zárate
and Sengers [2007]. For instance, the equal-time temperature-temperature
correlation function diverges for small wave numbers $k$ as $1/k^{4}$. In real
space, this corresponds to correlations that extend over the entire width of
the sample and decay on the same scale Ortiz de Zárate and Sengers [2007].
These correlations are generic in the sense that they do not require any fine
tuning of parameters, and they are not related to any broken symmetry. Rather,
they are the result of the coupling of the temperature fluctuations to the
diffusive shear velocity. This surprising result has been confirmed
theoretically by means of a variety of techniques Kirkpatrick et al. [1982b],
Ronis and Procaccia [1982], Ortiz de Zárate and Sengers [2007], and it has
been observed by many light scattering experiments, see Ref. Sengers et al.
[2016] and references therein.

Despite being well established and confirmed, this phenomenon raises several
questions that historically have not been emphasized. One is the fact,
mentioned above, that the relevant correlation functions in a NESS are not in
any obvious way related to response functions. Another question is whether the
long-ranged correlations are tied to thermal fluctuation effects, or whether
they are more generic and reflect some type of generalized rigidity Anderson
[1984] that is present even in the zero-temperature limit and also manifests
itself in the response of the system to external perturbations. Recent work on
classical fluids has suggested the latter Kirkpatrick et al. [2021], but a
relation between correlation functions and response theory has been lacking.

An interesting problem is the fate of the long-range correlations in quantum
fluids, where the above questions become particularly important. Measurements
of the correlation functions via light scattering is difficult even in the
classical case because of the very small scattering angles required. With
decreasing temperature the fluctuations become smaller, which makes it even
more desirable to observe the effect via response experiments, if feasible.
Specifically, there are two obvious types of correlation functions:
symmetrized, or anticommutator correlation functions, and anti-symmetrized, or
commutator ones. As functions of the wave vector ${\bm{k}}$ and the frequency
$\omega$ they are defined by spe , cla

$\displaystyle\frac{1}{2}\,\left\langle\left[\delta{\hat{A}}({\bm{k}}_{1},\omega_{1}),\delta{\hat{B}}({\bm{k}}_{2},\omega_{2})\right]_{+}\right\rangle$
$\displaystyle=$ $\displaystyle
V\delta_{{\bm{k}}_{1},-{\bm{k}}_{2}}\,2\pi\delta(\omega_{1}+\omega_{2})\,S^{\text{sym}}_{AB}({\bm{k}}_{1},\omega_{1})\
,$ (1a)
$\displaystyle\frac{1}{2\hbar}\,\left\langle\left[\delta{\hat{A}}({\bm{k}}_{1},\omega_{1}),\delta{\hat{B}}({\bm{k}}_{2},\omega_{2})\right]_{-}\right\rangle$
$\displaystyle=$ $\displaystyle
V\delta_{{\bm{k}}_{1},-{\bm{k}}_{2}}\,2\pi\delta(\omega_{1}+\omega_{2})\,\chi_{AB}^{\prime\prime}({\bm{k}}_{1},\omega_{1})\
.$ (1b)

Here $\delta\hat{A}$, $\delta\hat{B}$ are operator-valued fluctuations of
observables, $[\ ,\ ]_{+}$ and $[\ ,\ ]_{-}$ denote anticommutators and
commutators, respectively, $\langle\ldots\rangle$ denotes a quantum mechanical
expectation value plus an ensemble average, and $V$ is the system volume. The
two are related by a detailed-balance factor det

$S^{\text{sym}}_{AB}({\bm{k}},\omega)=\chi^{\prime\prime}_{AB}({\bm{k}},\omega)\,\coth(\omega/2T)\
,$ (2)

The corresponding static correlation functions are

$\displaystyle S^{\text{sym}}_{AB}({\bm{k}})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\,S^{\text{sym}}_{AB}({\bm{k}},\omega)\
,$ (3a) $\displaystyle\chi_{AB}({\bm{k}})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{d\omega}{\pi}\,\chi_{AB}^{\prime\prime}({\bm{k}},\omega)/\omega\
.$ (3b)

The symmetrized correlation functions $S^{\text{sym}}_{AB}$ are observable by
means of scattering experiments Forster [1975]. The physical meaning of the
antisymmetrized correlation functions $\chi^{\prime\prime}_{AB}$ is a priori
less obvious. In an equilibrium system, where the correlations are generically
short-ranged, they describe the linear response of the system to external
fields. That is, the equilibrium fluctuations determine the linear response,
which to second order in the external field yields the energy dissipated by
the system. This is the content of the fluctuation-dissipation theorem Nyquist
[1928], Callen and Welton [1951]. In a NESS, the detailed-balance relation (2)
still holds, but the commutator correlations functions no longer describe the
linear response and the usual fluctuation-dissipation theorem breaks down.
Another obvious question pertains to the nature of the transverse velocity
fluctuations, which changes from diffusive to propagating in the low-
temperature limit, and to their coupling to the temperature fluctuations. At a
technical level, in the collisionless regime of a quantum fluid a hydrodynamic
description no longer applies, and one has to work with a quantum kinetic
theory.

In this Letter we identify the quantum analogs of the classical long-ranged
correlations. Our two main results are: (1) The long-ranged anticommutator
correlation functions are still related to response functions, but the
relation is not a simple proportionality. Rather, the non-equilibrium
contributions to the correlation functions are related to the multilinear
response of products of observables. (2) In a modified form, the long-ranged
correlations extend to zero temperature. This shows that they are not tied to
thermal fluctuations, although thermal fluctuations can be used to probe them.
Rather, they are an inherent long-wavelength property of the NESS and indeed
reflect a novel type of generalized rigidity that does not become weaker with
decreasing temperature. We will start by explaining the second result, and
then demonstrate the first one.

Consider a fluid confined between two plates that is subject to a constant
temperature gradient ${\bm{\nabla}}T$ in the $z$-direction. Let ${\bm{k}}$ be
the wave-vector of a temperature fluctuation, and let the wave-vector space be
spanned by three mutually perpendicular unit vectors
$\hat{\bm{k}}={\bm{k}}/k$, $\hat{\bm{k}}_{\perp}$, and
$\hat{\bm{k}}_{\perp}^{(2)}$ such that $\hat{\bm{k}}_{\perp}$ is coplanar with
${\bm{k}}$ and ${\bm{\nabla}}T$, see Fig. 1. For our purposes the temperature
gradient will appear in the combination
$\nabla_{\\!\perp}T\equiv\hat{\bm{k}}_{\perp}\\!\cdot\\!{\bm{\nabla}}T$. For
definiteness we consider a fermionic quantum fluid (e.g., conduction electrons
in metals), but analogous effects must be present in bosonic fluids as well.
Let $\tau$ be the relaxation time, $v_{\text{F}}$ the Fermi velocity, and
$\omega$ the frequency. We need to distinguish between the hydrodynamic
regime, where $\omega<1/\tau$ or $v_{\text{F}}k<1/\tau$, and the collisionless
regime where $\omega>1/\tau$. The latter is in general subdivided into regimes
where $\omega<T$ and $\omega>T$, respectively, see Fig. 2 Pla .

Figure 1: A fluid subject to a constant temperature gradient in the
$z$-direction between two parallel confining plates. Figure 2: Frequency
regimes.

We are interested in the correlation functions for the temperature $T$ and the
component $u_{\perp}\equiv\hat{\bm{k}}_{\perp}\cdot{\bm{u}}_{\perp}$ of the
fluid shear velocity ${\bm{u}}_{\perp}$. In what follows we first state our
results and then comment on their derivations. A full account of the
derivations will be given elsewhere Kirkpatrick and Belitz [2023a].

Hydrodynamic Regime: In the hydrodynamic regime the result for
$S^{\text{sym}}_{TT}$ is the same as in a classical fluid Kirkpatrick et al.
[1982a]:

$\displaystyle S^{\text{sym}}_{TT}({\bm{k}},\omega)$ $\displaystyle=$
$\displaystyle\frac{2T^{2}}{c_{V}}\,\frac{D_{\text{T}}k^{2}}{\omega^{2}+(D_{\text{T}}k^{2})^{2}}$
(4a)
$\displaystyle\hskip-40.0pt+\,(\nabla_{\\!\perp}T)^{2}\,\frac{2T}{\rho}\,\frac{\nu
k^{2}}{(\omega^{2}+(D_{\text{T}}k^{2})^{2})(\omega^{2}+(\nu
k^{2})^{2})}\qquad$ The first term is the standard equilibrium contribution,
and the second term is the long-ranged NESS contribution. $c_{V}$ and $\rho$
are the spatially averaged specific heat per volume and the fluid density,
respectively, and $D_{\text{T}}$ and $\nu$ are the spatially averaged heat
diffusion coefficient and kinematic viscosity, respectively. The velocity-
velocity correlation function is the same as in equilibrium,
$S^{\text{sym}}_{u_{\perp}u_{\perp}}({\bm{k}},\omega)=\frac{2T}{\rho}\,\frac{\nu
k^{2}}{\omega^{2}+(\nu k^{2})^{2}}\ ,$ (4b) and the mixed symmetrized
correlation functions are $\displaystyle
S^{\text{sym}}_{u_{\perp}T}({\bm{k}},\omega)$ $\displaystyle=$ $\displaystyle
S^{\text{sym}}_{Tu_{\perp}}(-{\bm{k}},-\omega)$ (4c) $\displaystyle=$
$\displaystyle\nabla_{\\!\perp}T\,\frac{1}{\rho}\,\frac{\nu
k^{2}}{\omega^{2}+(\nu k^{2})^{2}}\,\frac{2T}{i\omega+D_{\text{T}}k^{2}}\ .\ $
The corresponding commutator correlation functions are
$\chi^{\prime\prime}_{AB}({\bm{k}},\omega)=S^{\text{sym}}_{AB}({\bm{k}},\omega)\,\omega/2T\qquad(A,B=T,u_{\perp})\
.$ (4d)

The static correlation functions are

$\displaystyle S^{\text{sym}}_{TT}({\bm{k}})$ $\displaystyle=$
$\displaystyle\frac{T^{2}}{c_{V}}+\frac{T(\nabla_{\\!\perp}T)^{2}}{\rho
D_{\text{T}}(\nu+D_{\text{T}})\,k^{4}}\ ,$ (5a) $\displaystyle
S^{\text{sym}}_{u_{\perp}u_{\perp}}({\bm{k}})$ $\displaystyle=$ $\displaystyle
T/\rho\ ,$ (5b) $\displaystyle S^{\text{sym}}_{u_{\perp}T}({\bm{k}})$
$\displaystyle=$ $\displaystyle
S^{\text{sym}}_{Tu_{\perp}}(-{\bm{k}})=\nabla_{\\!\perp}T\,\frac{T/\rho}{(\nu+D_{\text{T}})k^{2}}\
,$ (5c) and
$\chi_{AB}^{\prime\prime}({\bm{k}})=S^{\text{sym}}_{AB}({\bm{k}})/T\qquad(A,B=T,u_{\perp})\
.$ (5d)

Collisionless Regime: In the collisionless regime the expressions for the
dynamic correlation functions are lengthy, and we give only the results for
the static temperature-temperature correlations. To leading order as $T\to 0$
we have

$\chi_{TT}({\bm{k}})=\frac{1}{N_{\text{F}}}\,\frac{3}{\pi^{2}}\left[1+\frac{\pi^{2}}{12}\,(2\pi^{2}-3)\,\frac{(\nabla_{\\!\perp}T)^{2}}{\epsilon_{\text{F}}^{2}k^{2}}\right]\
$ (6)

everywhere in the collisionless regime. The symmetrized correlation function
is $S^{\text{sym}}_{TT}({\bm{k}})=T\chi_{TT}({\bm{k}})$ in Regime #1 in Fig. 1
and, apart from a factor of $O(1)$, $S^{\text{sym}}_{TT}({\bm{k}})\approx
v_{\text{F}}k\,\chi_{TT}({\bm{k}})$ in Regime #2.

To derive these results we note that in the hydrodynamic regime the usual
Navier-Stokes equations are applicable, as they are based on very general
physical principles and hold for quantum fluids as well as for classical ones.
In order to calculate correlation functions, they need to be augmented by
fluctuating forces Landau and Lifshitz [1966], Kirkpatrick and Belitz [2022].
Of the various nonlinearities we need to keep only the crucial coupling
between the temperature fluctuation and the transverse velocity, which turns
into a linear term in the presence of a fixed temperature gradient. We can
further ignore pressure fluctuations, which lead to sound waves (or plasmons
in a charged Fermi liquid) that are much faster than the diffusive shear
fluctuations. The equations for the operator-valued temperature and shear
velocity then are

$\left(-i\omega+\nu{\bm{k}}^{2}\right){\hat{u}}_{\perp}({\bm{k}},\omega)={\hat{P}}_{\perp}({\bm{k}},\omega)\
,\qquad\quad$ (7a)
$\left(-i\omega+D_{T}{\bm{k}}^{2}\right){\hat{T}}({\bm{k}},\omega)+(\nabla_{\\!\perp}T){\hat{u}}_{\perp}({\bm{k}},\omega)={\hat{Q}}({\bm{k}},\omega)\
.$ (7b)

Here the transport coefficients $\nu$ and $D_{\text{T}}$ are to be understood
as having been spatially averaged. ${\hat{P}}_{\perp}$ and $\hat{Q}$ are
operator-valued fluctuating forces that are Gaussian distributed with zero
mean. Their second moments can be obtained from the somewhat more general
expressions derived in Ref. Kirkpatrick and Belitz [2022]. The
$u_{\perp}$-$u_{\perp}$ correlation then is determined by the $P_{\perp}$
correlation and given by Eq. (4b). This is the same as in equilibrium since
the temperature does not couple into the $u_{\perp}$ equation. The equilibrium
part of the $T$-$T$ correlation is given by the $Q$ correlation, whereas the
non-equilibrium part, as well as the mixed $T$-$u_{\perp}$ correlation, is
given in terms of the $u_{\perp}$ fluctuation, with $\nabla_{\\!\perp}T$
serving as a coupling constant. This yields the second term in Eq. (4a), and
Eq. (4c). The leading singularity in the static $T$-$T$ correlation function
is $1/k^{4}$, and $1/k^{2}$ in the $T$-$u_{\perp}$ correlation function, as in
a classical fluid. This is because the hydrodynamic equations are the same in
either case.

In order to properly describe the collisionless regime one has to work with
the fluctuating quantum kinetic theory developed in Ref. Kirkpatrick and
Belitz [2022]. This procedure will be reported elsewhere Kirkpatrick and
Belitz [2023a]. However, the qualitative features of the results can be
obtained from simple scaling arguments as follows. In the collisionless regime
the diffusive inverse propagators of the form ${\cal
D}^{-1}({\bm{k}},\omega)=\omega+iDk^{2}$ that appear on the left-hand sides of
Eqs. (7), where $D$ can represent either $\nu$ or $D_{\text{T}}$ effectively
turn into propagating zero modes of the form $\omega\mp ck$, with $c\approx
v_{\text{F}}$ the propagation speed. The transport coefficients thus
effectively become singular functions of the wave number and scale as $D\sim
v_{\text{F}}/k$. The low-temperature result for $\chi_{TT}$, Eq. (6), then
follows from Eq. (5a) by replacing $\nu,D_{\text{T}}\to v_{\text{F}}/k$ and
using the low-temperature expression for the specific heat (prefactors, as
well as issues regarding reality and signs, require a much more detailed
analysis Kirkpatrick and Belitz [2023a]). For the symmetrized correlation
function one needs to consider the frequency integral in Eq. (3a) and
recognize that for the equilibrium contribution it needs to be cut off at
$\omega\approx D_{\text{T}}k^{2}$. In the limits $D_{\text{T}}k^{2}\ll T$ and
$D_{\text{T}}k^{2}\gg T$, and using again the effective scaling of
$D_{\text{T}}$ and $\nu$ explained above, one then obtains the relations
between $S^{\text{sym}}_{TT}({\bm{k}})$ and $\chi_{TT}({\bm{k}})$ given after
Eq. (6).

As mentioned above, the commutator functions $\chi^{\prime\prime}$ do not
determine the linear response of the NESS to external perturbations, in
contrast to an equilibrium system. In order to determine the response
functions we consider the averaged Navier-Stokes equations in the presence of
a field $h_{u_{\perp}}$ conjugate to $u_{\perp}$:

$\displaystyle\left(-i\omega+\nu k^{2}\right)u_{\perp}({\bm{k}},\omega)$
$\displaystyle=$
$\displaystyle\frac{\nu}{\rho}\,k^{2}h_{u_{\perp}}({\bm{k}},\omega)\ ,$ (8a)
$\displaystyle\left(-i\omega+D_{T}k^{2}\right)\delta T({\bm{k}},\omega)$
$\displaystyle+$
$\displaystyle(\nabla_{\\!\perp}T)u_{\perp}({\bm{k}},\omega)=$ (8b)
$\displaystyle\hskip-20.0pt\frac{1}{\rho}(\nabla_{\\!\perp}T)h_{u_{\perp}}({\bm{k}},\omega)\
.\qquad$

This form of the equations is obtained from time-dependent Ginzburg-Landau
theory Hohenberg and Halperin [1977]. It can also be derived by considering a
velocity-dependent force within kinetic theory Kirkpatrick and Belitz [2023b].
The field $h_{u_{\perp}}$ can be experimentally realized by imposing a shear
velocity on the system. The response functions $X_{AB}({\bm{k}},\omega)$ are
defined by

$\displaystyle u_{\perp}({\bm{k}},\omega)$ $\displaystyle=$ $\displaystyle
X_{u_{\perp}u_{\perp}}({\bm{k}},\omega)\,h_{u_{\perp}}({\bm{k}},\omega)\ .$
(9a) $\displaystyle\delta T({\bm{k}},\omega)$ $\displaystyle=$ $\displaystyle
X_{Tu_{\perp}}({\bm{k}},\omega)\,h_{u_{\perp}}({\bm{k}},\omega)\ ,$ (9b)

From Eqs. (8) we find

$\displaystyle X_{u_{\perp}u_{\perp}}({\bm{k}},\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\rho}\,\frac{\nu k^{2}}{-i\omega+\nu k^{2}}\ .$ (10a)
$\displaystyle X_{Tu_{\perp}}({\bm{k}},\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\rho}\,(\nabla_{\\!\perp}T)\,\frac{1}{-i\omega+D_{\text{T}}k^{2}}\,\frac{-i\omega}{-i\omega+\nu
k^{2}}\ .$ (10b)

From Eqs. (4b) and (10a) we see that the spectrum of $X_{u_{\perp}u_{\perp}}$,
$X^{\prime\prime}_{u_{\perp}u_{\perp}}({\bm{k}},\omega)=\text{Im}\,X_{u_{\perp}u_{\perp}}({\bm{k}},\omega)$,
equals the commutator correlation function
$\chi^{\prime\prime}({\bm{k}},\omega)$, as expected. However, the spectrum of
$X_{Tu_{\perp}}$,

$X_{Tu_{\perp}}^{\prime\prime}({\bm{k}},\omega)=-(\nabla_{\\!\perp}T)\,\frac{\omega(\nu
D_{\text{T}}k^{2}-\omega^{2})}{(\omega^{2}+D_{\text{T}}^{2}k^{4})(\omega^{2}+\nu^{2}k^{4})}\
,$ (11)

while showing the same scaling behavior as $\chi^{\prime\prime}_{Tu_{\perp}}$,
see Eqs. (4c, 4d), is not identical with the latter. In particular, the static
response function vanishes,

$X_{Tu_{\perp}}({\bm{k}})=\int\frac{d\omega}{\pi}\,\frac{X_{Tu_{\perp}}^{\prime\prime}({\bm{k}},\omega)}{\omega}=0\
,$ (12)

while the static correlation function is nonzero, see Eq. (5c). Nontheless,
the response functions still provide a way to directly measure the commutator
correlation functions, without relying on their symmetrized counterparts that
are suppressed at low temperatures: Considering Eqs. (9), (10), and (4c, 4d)
we have

$\delta
T({\bm{k}},\omega)\,u_{\perp}(-{\bm{k}},-\omega)=\frac{i}{\rho}\,\chi^{\prime\prime}_{Tu_{\perp}}({\bm{k}},\omega)\,|h_{u_{\perp}}({\bm{k}},\omega)|^{2}$
(13)

That is, the commutator $T$-$u_{\perp}$ correlation function describes the
bilinear response of the product of the temperature and the shear velocity to
the field conjugate to $u_{\perp}$. Similarly, the non-equilibrium part of the
commutator $T$-$T$ correlation function can be expressed as a bilinear
response to the field $h_{u_{\perp}}$. We define an observable

${\tilde{T}}({\bm{k}},\omega)=T({\bm{k}},\omega)-\frac{1}{\rho}\,(\nabla_{\\!\perp}T)\,\frac{1}{-i\omega+D_{\text{T}}k^{2}}\,h_{u_{\perp}}({\bm{k}},\omega)$
(14)

that obeys the equation

$\left(-i\omega+D_{T}k^{2}\right)\delta{\tilde{T}}({\bm{k}},\omega)+(\nabla_{\\!\perp}T)u_{\perp}({\bm{k}},\omega)=0\
.$ (15)

The physical interpretation of Eq. (15) is the heat equation with a streaming
term that contains the absolute shear velocity, whereas Eq. (8b) contains the
shear velocity relative to the the field $h_{u_{\perp}}$. The response of
${\tilde{T}}$ to the field $h_{u_{\perp}}$ is given by a response function

$X_{{\tilde{T}}u_{\perp}}({\bm{k}},\omega)=\frac{1}{\rho}\,(\nabla_{\\!\perp}T)\,\frac{-1}{-i\omega+D_{\text{T}}k^{2}}\,\frac{\nu
k^{2}}{-i\omega+\nu k^{2}}\ .$ (16)

From Eqs. (16), (10b), (4a), and (4d) we find that the non-equilibrium part of
the commutator $T$-$T$ correlation function describes the bilinear response of
the product of $\tilde{T}$ and $T$ to the field $h_{u_{\perp}}$:

$\delta{\tilde{T}}({\bm{k}},\omega)\,\delta
T(-{\bm{k}},-\omega)=\frac{i}{\rho}\,\chi_{TT}^{\prime\prime\text{\,neq}}({\bm{k}},\omega)\,|h_{u_{\perp}}({\bm{k}},\omega)|^{2}$
(17)

The Eqs. (13) and (17) constitute our main result. They demonstrate that in a
NESS characterized by a constant temperature gradient the commutator
correlation functions for the temperature and the shear velocity are still
related to the response of the system to an external shear perturbation, even
though the usual fluctuation-dissipation theorem is not valid. In contrast to
the situation in equilibrium, where the correlation functions are the same as
the response functions, in a NESS the correlation functions are given by the
bilinear response of products of observables.

To summarize, we have established a relation between correlation functions and
response functions for a fluid in a NESS. The resulting NESS fluctuation
formulas (13) and (17) resemble the equilibrium fluctuation-dissipation
theorem, Eq. (2) in the hydrodynamic limit, with the symmetrized correlation
function replaced by the product of two averaged observables, and the
temperature replaced by the driving field squared. The driving field, which is
realized by an imposed shear velocity, one has experimental control over
regardless of how low the temperature is. We also have determined the quantum
analogs of the long-ranged correlations known to exist in a classical fluid in
a NESS. In the latter context we note that in a fermionic quantum fluid there
potentially (depending on the values of the Landau Fermi-liquid parameters)
are many other zero modes that also display $\omega\sim k$ scaling, in
addition to the shear velocity. These can change the prefactor of the
singularity, but not the scaling behavior. Also, in a charged Fermi liquid
(e.g., conduction electrons in metals) the first-sound mode turns into the
massive plasmon, so our approximation of neglecting pressure fluctuations is
valid a fortiori.

## References

  * Nyquist [1928] H. Nyquist.  Thermal agitation of electric charge in conductors.  _Phys. Rev._ , 32:110, 1928.
  * Callen and Welton [1951] H. B. Callen and T. A. Welton.  Irreversibility and generalized noise.  _Phys. Rev._ , 83:34, 1951.
  * Forster [1975] D. Forster.  _Hydrodynamic Fluctuations, Broken Symmetry, and Correlation Functions_.  Benjamin, Reading, MA, 1975.
  * Kirkpatrick et al. [1982a] T. R. Kirkpatrick, E. G. D. Cohen, and J. R. Dorfman.  Light scattering by a fluid in a nonequilibrium steady state. II. Large gradients.  _Phys. Rev. A_ , 26:995, 1982a.
  * Dorfman et al. [1994] J. R. Dorfman, T. R. Kirkpatrick, and J. V. Sengers.  Generic long-range correlations in molecular fluids.  _Ann. Rev. Phys. Chem._ , 45:213, 1994.
  * Ortiz de Zárate and Sengers [2007] J. M. Ortiz de Zárate and J. V. Sengers.  _Hydrodynamic fluctuations in fluids and fluid mixtures_.  Elsevier, Amsterdam, 2007.
  * Kirkpatrick et al. [1982b] T. R. Kirkpatrick, E. G. D. Cohen, and J. R. Dorfman.  Fluctuations in a nonequilibrium steady state: Basic equations.  _Phys. Rev. A_ , 26:950, 1982b.
  * Ronis and Procaccia [1982] D. Ronis and I. Procaccia.  Nonlinear resonant coupling between shear and heat fluctuations in fluids far from equilibrium.  _Phys. Rev. A_ , 26:1812, 1982.
  * Sengers et al. [2016] J. V. Sengers, J. M. Ortiz de Zárate, and T. R. Kirkpatrick.  Thermal fluctuations in non-equilibrium steady states.  In D. Bedeaux, S. Kjelstrup, and J. V. Sengers, editors, _Non-equilibrium thermodynamics with applications_ , chapter 3, page 39. RSC Publishing, Cambridge, 2016.
  * Anderson [1984] P. W. Anderson.  _Basic Notions of Condensed Matter Physics_.  Benjamin, Menlo Park, CA, 1984.
  * Kirkpatrick et al. [2021] T. R. Kirkpatrick, D. Belitz, and J. R. Dorfman.  Superfast signal propagation in fluids and solids in non-equilibrium steady states.  _J. Chem. Phys. B_ , 125:7499, 2021.
  * [12] $\chi^{\prime\prime}$ is the customary notation for the commutator correlation function Forster [1975], with the double prime indicating that this is the spectrum of a causal function.
  * [13] Note that the commutator correlations still exist in the classical limit, due to the factor of $\hbar$ in the definition. Alternatively, they can be defined in terms of Poisson brackets.
  * [14] In a system that is not spatially homogeneous the exact relation between $S^{\text{sym}}$ and $\chi^{\prime\prime}$ is more complicated. It reduces to Eq. (2) if the local temperature is replaced by its spatial average.
  * [15] In an ordinary Fermi liquid, $1/\tau\propto T^{2}/\epsilon_{\text{F}}$, with $\epsilon_{\text{F}}$ the Fermi energy. More generally, $1/\tau$ can scale with a smaller power of $T$, but $1/\tau\lesssim T$ always holds Hartnoll and Mackenzie [2022].
  * Kirkpatrick and Belitz [2023a] T. R. Kirkpatrick and D. Belitz, 2023a.  Unpublished results.
  * Landau and Lifshitz [1966] L. D. Landau and E. M. Lifshitz.  _Fluid Mechanics_ , chapter XVII.  Pergamon, Oxford, Third Revised English edition, 1966.  Many editions are missing this chapter.
  * Kirkpatrick and Belitz [2022] T. R. Kirkpatrick and D. Belitz.  Fluctuating quantum kinetic theory.  _Phys. Rev. B_ , 105:245147, 2022.
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin.  Theory of dynamic critical phenomena.  _Rev. Mod. Phys._ , 49:435, 1977.
  * Kirkpatrick and Belitz [2023b] T. R. Kirkpatrick and D. Belitz.  Velocity-dependent forces and non-hydrodynamic initial conditions in quantum and classical fluids.  _arXiv:2304.06157_ , 2023b.
  * Hartnoll and Mackenzie [2022] S. A. Hartnoll and A. P. Mackenzie.  Planckian dissipation in metals.  _Rev. Mod. Phys._ , 94:041002, 2022.