# Investigating the Suitability of Concept Drift Detection for Detecting
Leakages in Water Distribution Networks

Valerie Vaquet***Autors contributed equally, Fabian Hinder,∗ and Barbara
Hammer

(Machine Learning Group, Bielefeld University, Germany)

###### Abstract

Leakages are a major risk in water distribution networks as they cause water
loss and increase contamination risks. Leakage detection is a difficult task
due to the complex dynamics of water distribution networks. In particular,
small leakages are hard to detect. From a machine-learning perspective,
leakages can be modeled as concept drift. Thus, a wide variety of drift
detection schemes seems to be a suitable choice for detecting leakages. In
this work, we explore the potential of model-loss-based and distribution-based
drift detection methods to tackle leakage detection. We additionally discuss
the issue of temporal dependencies in the data and propose a way to cope with
it when applying distribution-based detection. We evaluate different methods
systematically for leakages of different sizes and detection times.
Additionally, we propose a first drift-detection-based technique for
localizing leakages.

##  1 Introduction

Clean and safe drinking water is a scarce resource in many areas. Almost 80%
of the world’s population is classified as having high levels of threat in
water security [Vörösmarty et al., 2010]. This will aggravate in the future as
due to climate change the already limited water resources will become more
restricted [Rodell et al., 2018]. Currently, across Europe, considerable
amounts of drinking water are lost due to leakages in the
system222https://www.eureau.org/resources/publications/1460-eureau-data-
report-2017-1/file. The issue has been identified by the European Union, which
put the topic on the political agenda recently333https://eur-
lex.europa.eu/legal-
content/EN/TXT/HTML/?uri=CELEX:32020L2184&from=ES##d1e40-1-1.

In order to ensure a reliable drinking water supply, there is a need for
reliable, safe, and efficient water distribution networks (WDNs). In addition
to avoiding water losses, a crucial requirement is to ensure the quality of
the drinking water, e.g. to avoid the spreading and growth of bacteria and
other contaminants. Leakages pose a major risk to water quality as they make
it possible for unwanted substances to enter the water system. Thus,
monitoring the system for leakages is an efficient tool to avoid water loss
and contamination [Eliades and Polycarpou, 2010, Lambert, 1994].

Due to complex network dynamics and changing demand patterns detecting both
small and large leakages is a challenging task. This is aggravated by the fact
that the available data is very limited. Usually, the precise network topology
remains unknown or the documentation contains errors. As smart meter
technologies are not widely distributed there is no real-time demand
information [Cardell-Oliver and Carter-Turner, 2021]. In realistic settings,
this leaves a set of scarce pressure and possibly flow measurements.

Commonly, existing leakage detection methodologies rely on replicating the
system of interest by means of hydraulic models and monitoring the
discrepancies between observations and modeled values. While these approaches
can provide reasonable detection when considering larger leakages, the
approaches struggle when facing smaller leakages [Vrachimis et al., 2022].
Besides, there are major problems hindering the usage of these applications in
real-world applications. First of all, rebuilding a hydraulic system requires
detailed information about the network layout, e.g. exact pipe diameters, pipe
roughness, and elevation levels. Frequently, this kind of information is not
available to water supply companies, as the pipe infrastructure has been in
place for a long time. Second, hydraulic simulations require real-time demands
as an input for the computations. This is not a realistic requirement as the
availability of smart metering systems is still limited due to the associated
cost and privacy concerns [Cardell-Oliver and Carter-Turner, 2021]. Finally,
the need to adapt the simulation arises whenever the system is changing. Next
to the hydraulic approaches, there are also a few machine learning-based
approaches that implement a similar strategy.

In this work, we focus on the problem of leakage detection from the
perspective of handling data streams containing temporal dependencies. More
precisely, we formalize leakages as concept drift and the problem of leakage
detection as drift detection. Our aim is to investigate the suitability of
drift detection for reliable leakage detection, whereby we focus on leakage of
all practically relevant sizes. Our approach is independent of the specific
WDN and requires only real-time pressure measurements. Thus, it is more
flexible and efficient in comparison to the family of hydraulic simulation-
based approaches.

This paper is structured as follows. First, we introduce WDNs and summarize
the main specifics of this domain (section 2). Afterward, we briefly summarize
the body of related work on leakage detection (section 3). In section 4, we
define concept drift and cover both model-loss-based and distribution-based
drift detection. Before evaluating the suitability of these methodologies for
leakage detection in section 6, we discuss the issue of temporal dependencies
in the data collected from WDNs and propose ways to account for those (section
5). Finally, we conclude our paper in section 7.

##  2 Water Distribution Networks

WDNs can be modeled as graphs consisting of nodes representing junctions and
undirected edges representing pipes as the flow direction of the water is not
pre-defined and can change over time due to different inputs and demands in
the network. As usually real-world water distribution systems are extended
over time, it is possible that multiple pipes lie parallel to each other
[Mala-Jetmarova et al., 2018]. As the systems are observed over time,
hydraulic quantities like pressure and flow describe the network graph at each
time step. They can be described by hydraulic formulas given that the network
is known in great detail. Next to the exact pipe layout, different parameters
like elevations, pipe diameters, and pipe roughness are required. Assuming
sensors are installed in $n$ nodes across the graph, for each time step $t$,
measurements $x_{t}=[x_{1}^{t}\dots x_{n}^{t}]$ are collected where
$x_{i}^{t}$ is the value at node $i$.

Real-world measurements of the complete hydraulic state of WDNs are not
available as it is not possible to measure the entire system. Even a precise
topology alongside measurements in many positions is usually not available
since making this information public would be a weakness of the system which
is part of the critical infrastructure. Even if this data were available, for
the task of leakage detection, no ground truth of when exactly a leakage
occurred at which position would be available. Thus, when developing and
evaluating methodologies monitoring WDNs one usually relies on simulated data.
Given the key parameters of a system (layout, elevations, pipe information)
alongside demand patterns realistic network states containing anomalies like
leakages can be simulated using simulation tools like EPANET [Rossman, 2000].

As WDNs are part of the critical infrastructure and are required to work
robustly and safely to ensure the health and well-being of the population,
additional requirements are put on monitoring tools, especially those using AI
technologies. Besides requirements concerning robustness, safety, and fairness
as formulated in the European AI-Act [European Commission, 2021], some
technical attributes of WDNs pose additional challenges to ML approaches. When
working with WDNs, only limited knowledge about the pipe system is available.
Usually, the exact properties of the pipes, e.g. their diameters and
roughness, and the different elevation levels are unknown. In some cases, even
the precise pipe layout is not documented reliably [Eggimann et al., 2017].
Note that these are available for a few benchmark networks, and thus benchmark
scenarios can be generated. However, when designing monitoring systems relying
on this kind of information strongly limits the applicability in practical
applications. As we will discuss in the next section, the lack of this kind of
knowledge limits the application of many state-of-the-art solutions.

Besides limited information about the system setup, the system is also
relatively opaque with respect to the real-time dynamics. Due to installation
costs and challenges regarding the power supply, the availability of pressure
and flow sensors in WDNs is very limited yielding readings at a fraction of
the nodes in the system. Data availability is even more limited for real-time
demand measurements, as households are very rarely equipped with smart meters
for drinking water due to costs and data privacy.

Another property of WDNs is the presence of cycling patterns in demands,
flows, and pressures. When working on ML approaches, one needs to account for
the presence of temporal dependencies. As the water demands of both industrial
and private customers are not constant overall the demands fluctuate over the
course of a day. Additionally, the demand patterns vary between work days and
weekends, and changes caused by different seasons and long-term developments,
e.g. climate change or the COVID pandemic, need to be expected. These
fluctuations increase the difficulty of leakage detection as especially
smaller leakages might be lost in the signals.

##  3 Related Work

The body of related work on leakage detection can be divided into methods
relying on a hydraulic model and very few ML-based approaches. Hydraulic
model-based methods generally aim to replicate the real-world system with a
hydraulic model [Hu et al., 2018]. Usually, the simulation results of the
hydraulic model are then compared to the observations. An anomaly is reported
if the residual of these methods is too large, which is determined either by a
threshold [Romero-Ben et al., 2022], a CUSUM approach [Steffelbauer et al.,
2022], or visual inspection [Marzola et al., 2022]. All these methods share
the downside that they require real-time demands and more information on the
network topology than is usually available [Vrachimis et al., 2022]. Besides,
they lack generalizability across WDNs as the hydraulic model is specifically
designed for one network and even needs adaptation if something changes within
this particular network. While these hydraulic-based approaches yield good
results considering large leakages they usually miss smaller ones [Vrachimis
et al., 2022].

There are few ML-based approaches for leakage detection [Daniel et al., 2022,
Laucelli et al., 2016, Romano et al., 2014]. However, many are only evaluated
on very small networks and/or lack realistic demands as input for the
simulation data they are tested on. Most of these approaches replace the
hydraulic model with some ML model following the same general idea of
residual-based anomaly detection, for example by using a threshold [Daniel et
al., 2022, Laucelli et al., 2016].

##  4 Detecting Concept Drift

Deploying machine learning-based systems in real-world scenarios, one needs to
account for all kinds of changes and ensure that the models reliably work even
if the observed environment changes. Thus, considerable research focuses on ML
in the presence of changes in the data-generating process, which are called
concept drift or drift for shorthand. In order to obtain a formal definition
of drift, we first need to define a so-called drift process[Hinder et al.,
2023c]:

######  Definition 1.

Let $\mathcal{T}=[0,1]$ and $\mathcal{X}=\mathbb{R}^{d}$. A _drift process_
$(P_{T},\mathcal{D}_{t})$ from the _time domain_ $\mathcal{T}$ to the _data
space_ $\mathcal{X}$ is a probability measure $P_{T}$ on $\mathcal{T}$
together with a Markov kernel $\mathcal{D}_{t}$ from $\mathcal{T}$ to
$\mathcal{X}$, i.e., for all $t\in\mathcal{T}$ $\mathcal{D}_{t}$ is a
probability measure on $\mathcal{X}$ and for all measurable
$A\subset\mathcal{X}$ the map $t\mapsto\mathcal{D}_{t}(A)$ is measurable. We
will just write $\mathcal{D}_{t}$ instead of $(P_{T},\mathcal{D}_{t})$ if this
does not lead to confusion.

Based on this a definition of drift can be obtained:

######  Definition 2.

Let $(P_{T},\mathcal{D}_{t})$ be a drift process. We say that
$\mathcal{D}_{t}$ has _drift_ iff

$\displaystyle\mathbb{P}_{T,S\sim
P_{T}}[\mathcal{D}_{T}\neq\mathcal{D}_{S}]=P_{T}^{2}(\\{(t,s)\in\mathcal{T}^{2}\mid\mathcal{D}_{t}\neq\mathcal{D}_{s}\\})>0.$

Next to a branch of research focusing on the task of keeping an ML model
accurate in the presence of drift (_online or stream learning_), another
branch focuses on _monitoring and analyzing data streams_ with respect to
drift. In this work, we will focus on the latter as we are interested in the
monitoring task of leakage detection.

Drift detection schemes can be generally categorized into _model-loss-based_
and _distribution-based_ approaches. While the latter directly investigates
the observed data, model-loss-based approaches first train a model and then
analyze its loss as a proxy for change in the data distribution. Here, the
rationale is that a drift event changes the data in a way that the model
cannot approximate well anymore, causing a decline in the model loss. This
approach is valid if used in online learning where the goal is to keep a model
which performs well at its prediction task. Yet, as argued by [Hinder et al.,
2023a, Hinder et al., 2023b] the relation between model-loss and drift is
rather loose – in case the model does not provide sufficient complexity to
approximate the data distribution well (i) the drift might stay undetected as
it is smoothed out by the model or in converse (ii) the model might change
because of irrelevant drifts, e.g. a change in the ratio of classes, causing a
change in the model-loss even though no real drift occurred. For these
reasons, when working on monitoring tasks, one should apply distribution-based
approaches.

However, frequently model-loss-based approaches are also used in monitoring
tasks. For instance, the residual-based strategy described in section 3 which
is used across most state-of-the-art approaches published on leakage detection
can be categorized as model-loss-based approach. Therefore, we will
investigate the suitability of both types of drift detection methods in this
work.

###  4.1 Model-loss-based Drift detection

Performing model loss-based drift detection the idea is to analyze the model
loss as a proxy for drift in the data distribution. In case of WDNs, there are
two reasonable inference tasks a model can perform as a proxy for the
prediction: Either one performs a _forecasting_ task where the goal is to
predict the measurement of next time step $x_{t+1}$ based on the sensor
measurements collected up to time $t$, or one performs an _interpolation_ task
where the goal is to predict one sensor by the measurements of all other
sensors, i.e. for each node position $i$, a model
$f_{i}:\mathbb{R}^{n-1}\rightarrow\mathbb{R},\,f_{i}(x^{t}_{\setminus
i})=\hat{x}^{t}_{i}$ is trained. $x^{t}_{\setminus i}$ only includes the
measurements of the other sensors $\\{s_{j}\\}_{j\in\\{1,\dots n\\},\,j\neq
i}$ at time $t$. The latter strategy has been employed as a virtual sensor
imputation strategy in case of sensor faults. Even very simple ML models could
successfully perform the interpolation task [Vaquet et al., 2022].

###  4.2 Distribution-based Drift detection

Most distribution-based approaches follow the strategy of _comparing two
samples_. This can be done by _statistical testing_ : Considering two windows
some test is applied to determine whether they are generated by the same
distribution. One particular popular test is the _Kolmogorov-Smirnow (KS)
test_ [Kolomogorov, 1933]. However, since this test is only working on one-
dimensional data, a feature-wise testing scheme is required. In contrast, the
_kernel two-sample test_ which relies on the maximum mean discrepancy (MMD)
can test distributions of a higher degree by computing a kernel matrix as a
descriptor [Gretton et al., 2006].

Another option to include vectorial data is to use a _virtual classifier_ as a
first step. This strategy follows the fact that machine learning classifiers
perform better than random guessing if they can exploit some properties of the
data. The idea is to consider the two windows as different classes in a
classification task. If the obtained model performs better than it would be by
chance, the distributions of the windows differ, i.e. a drift occurred. For
this methodology, any classifier can be chosen. However, the model class
determines the detection capabilities. We will consider the D3 detection
scheme [Gözüaçık et al., 2019] in our experiments.

A methodology not focusing on two samples but instead relying on _meta-
statistics_ is the shape drift detector [Hinder et al., 2021]. It computes the
discrepancy of two consecutive time windows which was introduced as the drift
magnitude that is sliding through the stream. Plotting those discrepancies
over time, in case of drift, this graph will take a particular shape that can
be used to detect and pinpoint it. We will not use the ShapeDD directly in our
experiments as its main purpose is to increase the performance of the MMD, but
we will make use of it for analysis purposes.

Finally, a _block-based_ detection scheme can be used. In contrast to all
aforementioned methods, it does not rely on comparing two data samples
obtained at different time points to test for drift but rather searches
directly for a dependency of data and time which was identified to be an
equivalent description of drift by [Hinder et al., 2020]. This task can be
performed by a standard independence test; in this work, we will make use of
the HSIC-test [Gretton et al., 2007] which is another kernel-based method.

As discussed in section 2, different kinds of daily, weekly, and seasonal
patterns have to be expected. These patterns introduce certain temporal
dependencies to the data. As already discussed, these patterns might increase
the difficulty of detecting leakages. Considering this from a theoretical
viewpoint, this problem can be summarized by the need to account for the
temporal dependencies when performing drift detection. Thus before, describing
the drift detection schemes used in our experiments, in the next section we
will analyze the temporal patterns in the data.

##  5 Temporal dependencies in the data

We already raised the issue of temporal dependencies in data collected from
WDNs. In this section, we will analyze the dataset which we will use in our
experiments later on. For this purpose, we will first briefly introduce the
dataset and provide our analysis.

###  5.1 L-Town Benchmark Data

Figure 1: Visualization of sensor data for one year (no leak).

In this work, we will consider the L-Town network since it is relatively
complex in comparison to other benchmarks and one year of realistic real-time
demands are available for this system allowing us to simulate realistic data
for our experimental evaluation. The L-Town network resembles parts of the old
town of Limassol, Cyprus. In our experiments, we consider area A consisting of
661 nodes and 764 edges with 29 optimally placed pressure sensors [Vrachimis
et al., 2022]. We run simulations with four different leakage sizes ranging
from 7mm to 19mm at all pipes using the ATMN package which builds on EPANET.
Each scenario contains data of 364 days with a measuring frequency of 15
minutes. We always consider one leakage per scenario which starts at some
point of the scenario and stays present until the scenario ends.

###  5.2 Analysis

Analyzing the data, as expected we observe daily, weekly, and seasonal
patterns. As visualized in fig. 1, the pressure follows a clear weekly pattern
as can be seen in the zoom-in subplot. In order to control those dependencies
we perform two analysis strategies: 1) subtracting the “standard week”, 2)
subtracting the values of the previous week.

Figure 2: Visualization of standard week (line), shaded area shows standard
deviation.

The standard week is the average over all weeks. We have illustrated it in
fig. 2. As can be seen, the standard deviation is comparably small implying
that the weekly patterns are relatively consistent throughout the year.
Furthermore, we observe the same strong daily patterns we have observed before
as well as a difference in the shape of two days, probably corresponding to
weekends.

Figure 3: Visualization of sensor residuals after subtracting standard (no
leak). Orange line marks mean trend across all sensors.

By subtracting the standard week from the original signals we obtain the
signals shown in fig. 3. The plot shows one examplary sensor reading (blue
line) as well as the minimal and maximal sensor reading at each given point in
time as shaded area. As can be seen, the feature runs across the entire range
implying very strong fluctuations. Furthermore, as can be seen in the zoomed-
in plot there is a change in fluctuation that follows a daily pattern. In
addition, we also added a trend line (orange) which follows a cosine shape.
This is a plausible finding as we expect a cyclic pattern across several years
that correlates with the seasons. However, this pattern may render change
detection schemes useless as it induces changes that are not caused by leaks.

Figure 4: Visualization of sensor residuals after subtracting value of last
week (no leak). Orange line marks mean trend across all sensors.

As an alternative, we considered subtracting the value of the last week rather
than a standard week. The results are illustrated in fig. 4. Due to the small
variance in the computation of the standard week, this is already a good proxy
for the standard week. However, it is better suited to cope with long-term
changes as can be seen from the trend line. Furthermore, we again observe
strong oscillations whose intensities follow a daily pattern. We will find
that this strategy is actually quite efficient in section 6.1.

Figure 5: Plot of shape curve for different windows sizes. Red line marks time
point of drift, orange crosses mark candidate time points (shape heuristic)
with transparency indicating value of MMD at that point.

From both analyses, we expect that we can easily cope with the periodic
patterns if we only compare the data on a by-week basis. This is due to the
fact that there is a noticeable difference between the values of weekends and
weekdays so that day-wise is too short to resolve this dependency.
Furthermore, longer periods of time will be strongly affected by the seasonal
trends. We will further discuss those ideas in the next section.

###  5.3 Coping with temporal dependencies in the data

We observed substantial temporal patterns in the data which we need to account
for when utilizing drift detection schemes for the task of detecting leakages.
For model-loss-based drift detection approaches we assume that the models can
generalize well. Thus, in this setting, no additional actions need to be
taken. In contrast, when using distribution-based schemes, we need to
carefully incorporate our knowledge of the different temporal cycles in the
data to successfully detect leakages.

In preliminary experiments, we used a pre-processing technique, which
subtracted a standard week to eliminate cycles in the data. However, this
strategy assumes that we can model this standard week successfully which
requires some leakage-free historical data. Since we aim to develop a
methodology that requires as little information as possible to generalize to
new networks, we additionally experimented with choosing the window sizes such
that the detection schemes do not suffer from seasonalities. Here, as
discussed before our idea is to eliminate the cyclic patterns by choosing
exactly one week per window. Thereby daily and weekly patterns are eliminated
while the windows are still small enough to not be affected by long-term
dependencies. Since this strategy resulted in better results while requiring
no additional information, we will use this option in our experimental
evaluation instead of performing a preprocessing step subtracting the standard
or the previously observed week.

We analyzed the efficiency of this approach by making use of the Shape Drift
Detector [Hinder et al., 2021] which postprocesses the MMD of two consecutive
sliding windows in order to find candidate drift points. The result of this
postprocessing is the shape curve that indicates a candidate drift point when
it changes sign from positive to negative. We illustrated the shape curve for
different window lengths and one leakage in fig. 5. As can be seen, a short
window of just one day indicates drifts at the weekends, while for a window
length of two weeks (or more), the drift induced by seasonal effects becomes
stronger than the effect of the leak. Furthermore, for slightly less than a
week the shape curve shows strong oscillations that correspond to the fact
that we in fact have a weekly reoccurring pattern.

##  6 Experiments

For all our experiments, we use the data benchmark which we described in
section 5.1. We will first investigate the suitability of model-loss-based
drift detection and focus on data distribution-based drift detection
later.444The experimental code is available at
https://github.com/FabianHinder/Drift-and-Water

###  6.1 Model-Loss-Based Drift Detection

Figure 6: Mean squarred error of the forecasting task for different leakage
sizes. Setting without leakage is reported as baselines. Line marks mean
value, shaded area is minimal value to mean+standard deviation across folds.
Figure 7: ROC-AUC for the model-loss-based drift detection for increasing
leakagesizes. Forecasting task. Solid line marks mean value, shaded area
mean$\pm$standard deviation, dashed line median.

To evaluate the model-loss-based detection schemes relying on forecasting and
interpolation schemes for different regression models. Namely, we consider
$k$NN, polynomial ridge regression, random forests, and linear ridge
regression. RBF-ridge and RBF-/Poly-/Linear-SVR were considered but discarded
after initial considerations due to weak performances in the regression task.
In our experiments, we first analyze how well the models perform on the
forecasting and interpolation task, and analyze their generalization
capabilities to out-of-sample examples, e.g. to scenarios containing leakages.
In the second step, we analyze how well the schemes are suited for detecting
leakages. To do so we checked the underlying assumption that the model would
perform better on the original training data (without leak) compared to the
leaky data. In order for this to facilitate a useful strategy we need to be
able to define a threshold $\theta$ such that $\text{MSE}(x_{t})>\theta$
indicates a leakage at time $t$ and vice versa. Considering this as a
classification problem with the classes “no leak” and “leak” we can apply the
ROC AUC score to evaluate the performance of our models for leakage detection.
In order to be more resilient to slowly growing leakages we do not consider
model updates. Thus, we end up with the following procedure:

  1. 1.

Select one fold: Extract two consecutive weeks from the baseline dataset

  2. 2.

Train the model on the data for the respective task

  3. 3.

Compute the errors of the model for the remaining year for each data point
$E_{0}$

  4. 4.

Compute the errors of the model for the entire year for each leakagelocation
and size $E_{1}$

  5. 5.

Compute the detection performance for this fold $\text{ROC-
AUC}([0,\dots,0\;,\;1\dots,1],E_{0}+E_{1})$

Recall that the ROC-AUC measures how well the obtained scores separate the
leaky and non-leaky setups. The score is 1 if the largest error without
leakage is smaller than the smallest error with leakage, it is 0.5 if the
assignment is random. Thus, the ROC-AUC provides a scale-invariant upper bound
on the performance of every concrete threshold. Furthermore, the ROC-AUC is
not affected by class imbalance.

#### Forecasting

Figure 8: Mean absolute weight vector (and standard deviation) of forecasting
task using ElasticNet. Features correspond to time points relative to the time
point that is to be predicted (red dot).

The results of our experiment evaluating the generalization ability of
forecasting models are summarized in fig. 6. We observe increasing errors for
increasing leakage diameters. We find that our models generalize to small
leakages, i.e. the MSE reported for the smallest leakage size is very similar
to the baseline error on scenarios not containing any leakage. This might pose
a problem for the detection of small leakages as the drift in the system is
not reflected in the model loss. The small errors observed in the leakage-free
scenarios indicate that the trained models can generalize well to the temporal
patterns in the data.

Figure 9: Mean squarred error of the interpolation task for different leakage
sizes. Setting without leakage is reported as baselines. Line marks mean
value, shaded area is minimal value to mean+standard deviation. Figure 10:
ROC-AUC for the model-loss-based drift detection for increasing leakage sizes.
Interpolation task. Solid line marks mean value, shaded area mean$\pm$standard
deviation, dashed line median.

Concerning the experiments on leakage detection, we report very low scores
which align with random chance for small leakage sizes across all regression
models. While the detection capabilities increase with increasing leakage
size, the scores are still very weak.

As can be seen in fig. 6 linear models achieve comparably good performances.
To get a better understanding of which features are important for the
prediction we train an ElasicNet, a sparse linear model, and analyze its
weight vector which is illustrated in fig. 8. As can be seen, besides the
value right before the value that is to be predicted, the values one and two
weeks prior to it are very important. Furthermore, there is a slight relevance
for one day prior to the event. All other values are zero, indicating that
only those four time points are relevant for the prediction. This aligns well
with our initial analysis of the temporal patterns in the data (section 5.2).
We also repeated the same analysis with Ridge regression and Random Forest
leading to similar results.

#### Interpolation

The results of our experiment evaluating the generalization ability of the
interpolation models are summarized in fig. 9. We observe considerably lower
errors than for the forecasting task. Again the errors increase with
increasing leakage diameters and the models generalize to small leakages. In
this setting, we find that the simple linear models perform much better than
the $k$NN and the random forest. This aligns with findings published in
[Vaquet et al., 2022] using simple linear models as virtual sensors.

The evaluation of the detection experiments is summarized in fig. 10. In
contrast to the forecasting setting, we observe reasonable ROC-AUC scores for
leakage sizes of about 15mm-19mm. However, again, small leakages pose a
difficult problem for model-loss-based drift detection schemes. Similar to the
previous experiment, we observe scores that are only marginally above those
that would be obtained by random guessing.

In summary, in our experiments, we find that model-loss-based detection
schemes are only suitable for detecting large leakages. This aligns with both
results of other works and the theory we briefly discussed in section 4. In
this use case the models generalize too well to out-of-distribution samples to
reliably detect small leakages.

###  6.2 Distribution-based Drift Detection

Figure 11: Performance of unsupervised drift detectors for different
displacements (discrepancy between split point (assumed time-point of drift)
and actual drift). Solid line marks mean value, shaded area mean$\pm$standard
deviation, dashed line median. Figure 12: Performance of unsupervised drift
detectors for different leakage sizes. Solid line marks mean value, shaded
area mean$\pm$standard deviation, dashed line median.

In this section, we describe our experiments with the different distribution-
based drift detection schemes introduced in section 4.2. As argued before,
choosing a suitable window size is crucial for the successful usage of
distribution-based drift detection schemes. Analyzing the system at hand we
decided to rely on two windows of one week each to eliminate the temporal
patterns from the data. In addition to the choice of window size, the position
of the split point, i.e. the border of the two considered windows, affects the
performance of the detection scheme. The larger the displacement from the
actual leakage onset time, the harder the detection will be. However, as it is
desirable to detect leakages as soon as possible, it would be desirable to
detect leakages even if the displacement is still large, e.g. if the split
point still lies before the leakage occurred. If we get a reasonably high
detection score, some kind of warning mechanism can be implemented where
warnings can be confirmed when more data arrives.

In fig. 11 we summarize the results of our distribution-based detection
schemes with the different models. Assuming the split point lies exactly on
the leakage, we report much better detection results than for the model-loss-
based detection strategies. We report smaller scores for DAWIDD which is to be
expected in this case as the model is not benefiting from our window size
choice due to its block-based nature. Thus, it is affected by the temporal
dependencies in the oscillations. While we obtain smaller scores for the
statistical test-based methods for small leakages, the virtual classifier-
based methods reliably detect leakages of a diameter of 7mm.

As assumed, concerning the position of the split point, we observe a decline
in the score with increasing split points. However, even for a displacement of
4 days, we obtain better scores than for the best model-based detection
schemes. In the presence of a displacement of 6 days, we still get reasonable
scores for large leakages when using the D3 detection scheme with a linear
model. Thus detecting large leakages can be realized very fast, while for
smaller leakages it takes a little more time to obtain a reliable warning.

These findings can be confirmed when analyzing fig. 12. One can see that apart
from the smallest leakage size, with a displacement of 4 days, i.e. 3 days
after the leakage occurred, most detection schemes yield a reasonable score,
with the D3 variants performing best. Using the linear version of D3 even the
smallest considered leakages can be detected with a delay of 5 days.

In conclusion, we found that distribution-based detection schemes outperform
model-based detection if the window size is chosen in a suitable way. The
methods are capable of detecting all leakage sizes. Large leakages can be
detected with a reasonably small delay. In practical applications implementing
a warning mechanism early on could be realized to react early on.

###  6.3 Leakage localization

Besides the task of leakage detection, i.e. determining whether or not there
is a leakage somewhere in the WDN, there is also the more specific task of
leakage localization, i.e. assuming there is a leakage somewhere determine the
pipe where it occurrs [Vrachimis et al., 2022]. This task is usually
considered much harder and is commonly approached by trying to find the
additional demand that is necessary to make the predicted pressures match the
observed pressures[Li et al., 2022, Daniel et al., 2022, Wang et al., 2022,
Marzola et al., 2022]. This of course again requires a lot of data usually not
available to us. As our drift detection approach performed quite well on the
detection task we consider the possibility to extend the methodology for
leakage detection. Here the idea is quite simple the closer a sensor is to the
leakage’s actual position, the stronger the influence and thus the drift,
leading to a particularly small $p$-value for that feature. As only the
Kolmogorov-Smirnov test operates feature-wise we only take it into account. We
make use of the same setup as before and narrow down the number of potential
leakage locations by searching the sensor node closest to it. The closest
sensor node is selected by choosing the one that has the smallest $p$-value,
i.e. is considered to be particularly drifting by the test.

In the following let $S$ be all sensor nodes, $s^{*}$ be the selected sensor
node, and $v$ be the node where there leakage actually occurred, even if it is
not the sensor node. Furthermore, $d$ denotes the graph distance in the WDN,
i.e. $d(a,b)$ is the length of the shortest path connecting $a$ and $b$. We
make use of three metrics: distance between selected and actual node (Dist.;
$d(s^{*},v)$), number of sensor nodes closer to the actual node (#Cls.;
$|\\{s\in S\mid d(s,v)<d(s^{*},v)\\}|$), and relative distance between actual
node, selected and optimal node (rel.D.; $d(s^{*},v)/\min_{s\in S}d(s,v)$)
which is normalized in contrast to the simple distance and smooth in contrast
to the closer node metric.

Table 1: Results of leakage localization size | Dist. | rel.D. | #Cls.
---|---|---|---
(mm) |  $\mu$$\pm$ | $\sigma$ |  $\mu$$\pm$ | $\sigma$ |  $\mu$$\pm$ | $\sigma$
7 |  10.1$\pm$ | 13.1 |  2.6$\pm$ | 4.9 |  3.3$\pm$ | 7.0
11 |  5.5$\pm$ | 4.7 |  1.3$\pm$ | 1.4 |  0.6$\pm$ | 2.0
15 |  5.1$\pm$ | 3.9 |  1.2$\pm$ | 1.2 |  0.5$\pm$ | 1.5
19 |  5.0$\pm$ | 3.6 |  1.2$\pm$ | 1.1 |  0.4$\pm$ | 1.4

The results are shown in table 1. As can be seen, the results are quite
promising. Furthermore, we observe that the precision is decreasing for
smaller leakage sizes, which is to be expected considering the results from
the last experiment.

##  7 Conclusion

In this work, we investigated the suitability of model-loss-based and
distribution-based drift detection methods. Combining distribution-based
detection with knowledge of WDNs, we provide detection schemes that
successfully detect leakages of all sizes with reasonable detection delays.
Analyzing model-loss-based techniques that are widely implemented in the water
domain, we confirmed theoretical results that raise the issue of the loose
connection between model loss and drift.

We assume that our work is not limited to WDNs but can also be realized for
anomaly detection in other critical infrastructure systems like gas or
electrical grids. In practical applications a further analysis of the leakages
is necessary – solely detecting leakages is not sufficient to take appropriate
actions. We proposed a first localization strategy that is directly from
detection efforts. Considering these follow-up tasks in more detail through
the lens of concept drift is an interesting path for further research.
Especially leakage localization requires further investigation as our initial
results still leave space for improvement and require an extension to the
actual leakage position not just the closest sensor node if applied in
practice. A theoretical justification for this approach would be desirable.

## Acknowledgements

We gratefully acknowledge funding from the European Research Council (ERC)
under the ERC Synergy Grant Water-Futures (Grant agreement No. 951424).

## References

  * European Commission, 2021 European Commission (2021).  Artificial Intelligence Act.
  * Cardell-Oliver and Carter-Turner, 2021 Cardell-Oliver, R. and Carter-Turner, H. (2021).  Activity-aware privacy protection for smart water meters.  In 8th ACM BuildSys, BuildSys ’21, page 31–40. Association for Computing Machinery.
  * Daniel et al., 2022 Daniel, I., Pesantez, J., Letzgus, S., Khaksar Fasaee, M. A., Alghamdi, F., Berglund, E., Mahinthakumar, G., and Cominola, A. (2022).  A Sequential Pressure-Based Algorithm for Data-Driven Leakage Identification and Model-Based Localization in Water Distribution Networks.  J. Water Resour. Plan. Manag., 148.
  * Eggimann et al., 2017 Eggimann, S., Mutzner, L., Wani, O., Schneider, M. Y., Spuhler, D., Moy de Vitry, M., Beutler, P., and Maurer, M. (2017).  The Potential of Knowing More: A Review of Data-Driven Urban Water Management.  Environmental Science & Technology, 51(5):2538–2553.
  * Eliades and Polycarpou, 2010 Eliades, D. G. and Polycarpou, M. M. (2010).  A Fault Diagnosis and Security Framework for Water Systems.  IEEE Transactions on Control Systems Technology, 18(6):1254–1265.
  * Gözüaçık et al., 2019 Gözüaçık, Ö., Büyükçakır, A., Bonab, H., and Can, F. (2019).  Unsupervised concept drift detection with a discriminative classifier.  In Proceedings of the 28th ACM international conference on information and knowledge management, pages 2365–2368.
  * Gretton et al., 2006 Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., and Smola, A. J. (2006).  A kernel method for the two-sample-problem.  In NIPS 2006, pages 513–520.
  * Gretton et al., 2007 Gretton, A., Fukumizu, K., Teo, C. H., Song, L., Schölkopf, B., and Smola, A. J. (2007).  A kernel statistical test of independence.  In NIPS 2007, pages 585–592.
  * Hinder et al., 2020 Hinder, F., Artelt, A., and Hammer, B. (2020).  Towards non-parametric drift detection via dynamic adapting window independence drift detection (DAWIDD).  In ICML 2020, volume 119, pages 4249–4259. PMLR.
  * Hinder et al., 2021 Hinder, F., Brinkrolf, J., Vaquet, V., and Hammer, B. (2021).  A shape-based method for concept drift detection and signal denoising.  In 2021 IEEE Symposium Series on Computational Intelligence (SSCI), pages 01–08. IEEE.
  * Hinder et al., 2023a Hinder, F., Vaquet, V., Brinkrolf, J., and Hammer, B. (2023a).  On the Change of Decision Boundary and Loss in Learning with Concept Drift.  In IDA XXI, volume 13876, pages 182–194. Springer Nature Switzerland, Cham.
  * Hinder et al., 2023b Hinder, F., Vaquet, V., Brinkrolf, J., and Hammer, B. (2023b).  On the Hardness and Necessity of Supervised Concept Drift Detection:.  In 12th ICPRAM, pages 164–175, Lisbon, Portugal. SCITEPRESS.
  * Hinder et al., 2023c Hinder, F., Vaquet, V., and Hammer, B. (2023c).  One or two things we know about concept drift – a survey on monitoring evolving environments.
  * Hu et al., 2018 Hu, C., Li, M., Zeng, D., and Guo, S. (2018).  A survey on sensor placement for contamination detection in water distribution systems.  Wireless Networks, 24(2):647–661.
  * Kolomogorov, 1933 Kolomogorov, A. (1933).  Sulla determinazione empirica di una legge didistribuzione.  Giorn Dell’inst Ital Degli Att, 4:89–91.
  * Lambert, 1994 Lambert, A. (1994).  Accounting for Losses: The Bursts and Background Concept.  Water and Environment Journal, 8(2):205–214.
  * Laucelli et al., 2016 Laucelli, D., Romano, M., Savić, D., and Giustolisi, O. (2016).  Detecting anomalies in water distribution networks using EPR modelling paradigm.  Journal of Hydroinformatics, 18(3):409–427.
  * Li et al., 2022 Li, Z., Wang, J., Yan, H., Li, S., Tao, T., and Xin, K. (2022).  Fast Detection and Localization of Multiple Leaks in Water Distribution Network Jointly Driven by Simulation and Machine Learning.  J. Water Resour. Plan. Manag., 148(9).
  * Mala-Jetmarova et al., 2018 Mala-Jetmarova, H., Sultanova, N., and Savic, D. (2018).  Lost in Optimisation of Water Distribution Systems? A Literature Review of System Design.  Water, 10(3):307.
  * Marzola et al., 2022 Marzola, I., Mazzoni, F., Alvisi, S., and Franchini, M. (2022).  Leakage Detection and Localization in a Water Distribution Network through Comparison of Observed and Simulated Pressure Data.  J. Water Resour. Plan. Manag., 148(1):04021096.
  * Rodell et al., 2018 Rodell, M., Famiglietti, J. S., Wiese, D. N., Reager, J., Beaudoing, H. K., Landerer, F. W., and Lo, M.-H. (2018).  Emerging trends in global freshwater availability.  Nature, 557(7707):651–659.
  * Romano et al., 2014 Romano, M., Kapelan, Z., and Savić, D. A. (2014).  Automated Detection of Pipe Bursts and Other Events in Water Distribution Systems.  J. Water Resour. Plan. Manag., 140(4):457–467.
  * Romero-Ben et al., 2022 Romero-Ben, L., Alves, D., Blesa, J., Cembrano, G., Puig, V., and Duviella, E. (2022).  Leak Localization in Water Distribution Networks Using Data-Driven and Model-Based Approaches.  J. Water Resour. Plan. Manag., 148(5).
  * Rossman, 2000 Rossman, L. A. (2000).  EPANET 2: users manual.  US Environmental Protection Agency. Office of Research and Development.
  * Steffelbauer et al., 2022 Steffelbauer, D. B., Deuerlein, J., Gilbert, D., Abraham, E., and Piller, O. (2022).  Pressure-Leak Duality for Leak Detection and Localization in Water Distribution Systems.  J. Water Resour. Plan. Manag., 148(3).
  * Vaquet et al., 2022 Vaquet, V., Artelt, A., Brinkrolf, J., and Hammer, B. (2022).  Taking Care of Our Drinking Water: Dealing with Sensor Faults in Water Distribution Networks.  In ICANN 2022, volume 13530, pages 682–693. Springer Nature Switzerland, Cham.
  * Vörösmarty et al., 2010 Vörösmarty, C. J., McIntyre, P. B., Gessner, M. O., Dudgeon, D., Prusevich, A., Green, P., Glidden, S., Bunn, S. E., Sullivan, C. A., Liermann, C. R., et al. (2010).  Global threats to human water security and river biodiversity.  nature, 467(7315):555–561.
  * Vrachimis et al., 2022 Vrachimis, S. G., Eliades, D. G., Taormina, R., Kapelan, Z., Ostfeld, A., Liu, S., Kyriakou, M., Pavlou, P., Qiu, M., and Polycarpou, M. M. (2022).  Battle of the Leakage Detection and Isolation Methods.  Journal of Water Resources Planning and Management, 148(12).
  * Wang et al., 2022 Wang, X., Li, J., Liu, S., Yu, X., and Ma, Z. (2022).  Multiple Leakage Detection and Isolation in District Metering Areas Using a Multistage Approach.  J. Water Resour. Plan. Manag., 148(6).


# Observation of the massive Lee-Fukuyama phason in a charge density wave
insulator

Soyeun Kim  Department of Physics, University of Illinois at Urbana-Champaign,
Urbana, IL 61801, USA  Materials Research Laboratory, University of Illinois
at Urbana-Champaign, Urbana, IL 61801, USA  Yinchuan Lv  Department of
Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
Materials Research Laboratory, University of Illinois at Urbana-Champaign,
Urbana, IL 61801, USA  Xiao-Qi Sun  Department of Physics, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Chengxi Zhao  Materials
Research Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL
61801, USA  Materials Science and Engineering Department, University of
Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Nina Bielinski
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL
61801, USA  Materials Research Laboratory, University of Illinois at Urbana-
Champaign, Urbana, IL 61801, USA  Azel Murzabekova  Department of Physics,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Materials
Research Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL
61801, USA  Kejian Qu  Department of Physics, University of Illinois at
Urbana-Champaign, Urbana, IL 61801, USA  Materials Research Laboratory,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Ryan A.
Duncan  Stanford PULSE Institute, SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, United States  Stanford Institute for Materials and Energy
Sciences, SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
Quynh L. D. Nguyen  Stanford PULSE Institute, SLAC National Accelerator
Laboratory, Menlo Park, CA 94025, United States  Stanford Institute for
Materials and Energy Sciences, SLAC National Accelerator Laboratory, Menlo
Park, CA 94025, USA  Mariano Trigo  Stanford PULSE Institute, SLAC National
Accelerator Laboratory, Menlo Park, CA 94025, United States  Stanford
Institute for Materials and Energy Sciences, SLAC National Accelerator
Laboratory, Menlo Park, CA 94025, USA  Daniel P. Shoemaker  Materials Research
Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
Materials Science and Engineering Department, University of Illinois at
Urbana-Champaign, Urbana, IL 61801, USA  Barry Bradlyn  Department of Physics,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  Fahad
Mahmood<EMAIL_ADDRESS>Department of Physics, University of Illinois at
Urbana-Champaign, Urbana, IL 61801, USA  Materials Research Laboratory,
University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA

The lowest-lying fundamental excitation of an incommensurate charge density
wave (CDW) material is widely believed to be a massless phason – a collective
modulation of the phase of the CDW order parameter. However, as first pointed
out by Lee and Fukuyama [1], long-range Coulomb interactions should push the
phason energy up to the plasma energy of the CDW condensate, resulting in a
massive phason and a fully gapped spectrum. Whether such behavior occurs in a
CDW system has been unresolved for more than four decades. Using time-domain
THz emission spectroscopy, we investigate this issue in the material
(TaSe4)2I, a classical example of a quasi-one-dimensional CDW insulator. Upon
transient photoexcitation at low temperatures, we find the material strikingly
emits coherent, narrow-band THz radiation. The frequency, polarization and
temperature-dependence of the emitted radiation imply the existence of a
phason that acquires mass by coupling to long-range Coulomb interaction. Our
observations constitute the first direct evidence of the massive “Lee-
Fukuyama” phason and highlight the potential applicability of fundamental
collective modes of correlated materials as compact and robust sources of THz
radiation.

The fundamental collective modes (amplitudon and phason) of a broken symmetry
ordered state (Fig. 1a) have been key in establishing foundational theories
across various fields of physics, including gauge theories in particle
physics; and superconductors, antiferromagnets and charge-density wave (CDW)
materials in condensed matter physics [2, 3, 4]. The phason is typically
massless, in accordance with Goldstone’s theorem, which necessitates the
emergence of a massless boson for a broken symmetry in systems in which the
ground or vacuum state is continuously degenerate. A prominent exception
occurs in superconducting systems. Here, even though the ground state is
continuously degenerate, the long-range Coulomb interaction pushes the
longitudinal phason up to the plasma frequency [5, 6], and so a massless
Goldstone boson does not exist and the low-lying excitation spectrum is fully
gapped. This behavior is the celebrated Anderson-Higgs mechanism which
established the deep connection between symmetry breaking and gauge fields,
and ultimately explained how all fundamental particles acquire mass from
interactions with the Higgs field.

Unlike superconductivity, an incommensurate CDW is believed to have a massless
phason, typically understood in terms of softening of a longitudinal acoustic
phonon branch around the CDW wavevector $\vec{q}_{\text{CDW}}$ (Fig. 1b).
Below the CDW transition temperature ($T_{\text{CDW}}$), this mode softening
results in the linear-in-wavevector, zero-gap dispersion of the phason,
implying that the CDW can freely slide for excitation wavevector $\vec{q}=0$.
In any real material, however, random impurities and disorder restrict this
sliding motion, leading to a small gap in the phason dispersion (the pinning
frequency) (Fig. 1b). Thus, the sliding CDW motion can only be observed if a
strong enough electric field is applied to first depin the phason. The
resulting sliding motion of the CDW can then be measured in DC transport
experiments as has been done in various systems [7].

Figure 1: Collective modes of an incommensurate charge-density-wave (CDW)
phase. (a) Ginzburg-Landau free energy ($\mathcal{F}$) and real space
representation of the amplitudon (dashed arrow) and the phason (solid line
arrow) for a CDW order parameter. (b-c) Dispersion relations of the CDW
collective modes. Below $T_{\text{CDW}}$, the otherwise undistorted acoustic
phonon (LA) softens near $\vec{q}_{\text{CDW}}$ and renormalizes into the
amplitudon and phason bands. (b) At moderate $T$ ($\lesssim T_{\text{CDW}}$)
the long-range Coulomb repulsion $U$ is screened by quasiparticles. (c) When
the system is cooled well below $T_{\text{CDW}}$, the screening weakens and
the spectral weight from the massless (acoustic) phason is transferred to the
massive (IR-active) phason of frequency $\omega_{\text{LF}}$.

However, the above phenomenology of a massless phason (or disorder pinned
phason at low frequency) assumes the absence of long-range Coulomb
interactions ($U$). This assumption is believed to be valid because the
presence of normal electrons at a non-zero temperature can screen $U$.
However, if $U$ were sufficiently strong, or if the density of normal
electrons were sufficiently low, then the CDW phason at $\vec{q}=0$ should be
pushed to higher energies (even above the amplitudon energy) (Fig. 1c). This
behavior was highlighted in the seminal work of Lee, Rice and Anderson (LRA)
on CDW dynamics [8, 9] and soon formalized by Lee and Fukuyama [1] who noted
the similarity of the emergence of the massive CDW phason with the Anderson-
Higgs mechanism in a superconductor. Later works [10, 11] predicted that the
massive (optical) phason could indeed dominate over the massless (acoustic)
phason at sufficiently low temperatures where charged quasiparticles cannot
sufficiently screen $U$ [10, 11]. We note that in superconductors the plasma
frequency is much larger than the single particle gap, rendering the phase
mode unobservable deep in the superconducting phase. In a CDW system, however,
the relevant scale is the plasma frequency of the condensate, which can lie
far below the single-particle gap. To date, direct experimental evidence of
the massive “Lee-Fukuyama” phason in CDW systems has been elusive.

Here, we present evidence for the generation and detection of a massive Lee-
Fukuyama phason in the CDW insulator (TaSe4)2I, a quasi-1D material [12] that
undergoes an incommensurate CDW transition below $T_{\text{CDW}}\sim$ 260 K
[13, 14, 15, 16] with a CDW gap for single particle excitations of
2$\Delta_{\text{CDW}}\approx 250-300$meV [17, 18, 19, 20]. (TaSe4)2I is unique
among quasi-1D CDW systems due to its unusually high resistivity in the low
temperature insulating state [21], such that the long-range Coulomb
interaction can remain unscreened and the Lee-Fukuyama phason can have have
significant spectral weight. To investigate the collective modes of the CDW
order in (TaSe4)2I, we performed time-domain THz emission spectroscopy using
an ultrafast photoexcitation ‘pump’ pulse with an energy of 1.2 eV ($\lambda$
= 1030 nm) (See Fig. 2a for the experimental geometry and [21] for further
details of the experimental setup). Note that the photoexcitation energy here
is greater than 2$\Delta_{\text{CDW}}$ and so the pump pulse initially creates
single-particle excitations across the CDW gap. As (TaSe4)2I is structurally
chiral (SG. 97) and lacks inversion symmetry, this photoexciation generates a
transient current with a duration of a few picoseconds (ps). Such phenomena is
well-known as a photo-galvanic or a photo-Dember effect, both of which can
occur in systems lacking inversion symmetry [22, 23]. The transient current
then results in a short ps-duration burst of THz radiation in the far-field
which we measure in the time-domain using standard electro-optical sampling
(EOS) [21].

Figure 2: Terahertz (THz) emission from (TaSe4)2I. (a) The geometry of the
sample orientation, incident light (pump) and emitted THz polarizations. The
TaSe4 chains ($c$-axis) are oriented along the $z$-axis, and the pump beam is
${p}$-polarized with incident angle $\theta_{i}=45^{\circ}$. The plane of
incidence is represented as the gray shaded region. (b) The THz emission
signal ($E_{\text{THz}}$) as a function of time ($t_{\text{delay}}$), measured
at 7 K. The Fourier transform amplitude ($A_{\text{FT}}$) is plotted in the
inset. The 0.23 THz mode is marked with an asterisk (*). (c) The $p$\- and
$s$-polarized components of the THz emission signal, $E_{\text{THz}}^{p}$
(upper) and $E_{\text{THz}}^{s}$ (lower), respectively. $E_{\text{THz}}^{p}$
is shown with an offset.

Figure 2b shows the measured time profile of the THz electric field
(${E}_{\text{THz}}(t)$) emitted from (TaSe4)2I upon photoexcitation at $T$ = 7
K $\ll$ $T_{\text{CDW}}$. Here the pump is $p$-polarized, with an electric
field component along the quasi-1D chains of (TaSe4)2I. Two features are
immediately evident in $E_{\text{THz}}(t)$: a transient peak around
$t_{\text{delay}}=0$ ps followed by a long-lived coherent oscillation that
lasts over $\sim$ 80 ps. In the frequency domain (Fig. 2b inset), the
transient peak around $t_{\text{delay}}\sim 0$ ps manifests as a broad
distribution over frequencies from 0.1 to 2 THz, while the long-lived coherent
oscillation manifests as a sharp peak centered at 0.23 THz. For the rest of
this work, we refer to the transient ($t_{\text{delay}}\sim 0$ ps) peak as
$E_{\text{tr}}$ and the coherent oscillation as $E_{\text{osc}}$. As noted
above, $E_{\text{tr}}$ is what we typically expect from such an experiment due
to a photo-galvanic or a photo-Dember effect. The transient current can be
estimated from the measured $E_{\text{tr}}$, and is greater than the current
necessary to depin the dynamics of CDW order in (TaSe4)2I [21].

Figure 3: Temperature ($T$) and pump fluence ($F$) dependence of the 0.23 THz
mode. (a) THz emission signal as a function of delay time at a few selected
temperatures. The signals at different temperatures are offset for clarity.
(b) Fourier transforms of the traces shown in (a). (c) Spectral weight ($SW$)
of the 0.23 THz mode as a function of temperature, normalized to the $SW$ at 7
K. Dashed line indicates 0.3$T_{\text{CDW}}\sim$ 80 K. (d) THz emission signal
at 7 K as a function of delay time at a few selected pump fluences. The
signals at different fluences are shown with offsets for clarity. (e) Fourier
transforms for the traces shown in (d). (f) $SW$ of the 0.23 THz mode as a
function of pump fluence, normalized to the maximum value of the $SW$.

We focus on the observed narrow-band THz emission $E_{\text{osc}}$ since this
aspect of the data is particularly striking. $E_{\text{osc}}$ lasts well over
80 ps – much longer than the typical few ps duration signal expected from
semiconductors [23] and semimetals [24, 25, 26] in THz emission experiments.
Another unusual feature of $E_{\text{osc}}$ is the observed waveform envelope
in the time domain which appears to gradually increase in magnitude starting
at $t_{\text{delay}}$ = 0 ps. In impulsive excitation ultrafast experiments,
the signal is usually peaked at $t_{\text{delay}}$ = 0 from where it
exponentially decays. Additionally, while the measured transient peak
$E_{\text{tr}}$ is both horizontally and vertically polarized, the coherent
oscillation $E_{\text{osc}}$ is only horizontally polarized (Fig. 2c). In our
experimental geometry, this corresponds to $E_{\text{osc}}$ being polarized
along the quasi-1D chains of (TaSe4)2I as shown in Fig. 2a.

To investigate the origin of the radiating mode, we study the $E_{\text{osc}}$
spectra at various different temperatures (Figs. 3a,b). The $E_{\text{osc}}$
signal is strongest at 7 K and gradually decreases with increasing
temperature, showing a sudden drop at $\sim$ 80 K (at about 30% of
$T_{\text{CDW}}$). This feature is also clear in the integrated spectral
weight of the 0.23 THz mode as a function of temperature shown in Fig. 3c.
Note that the initial transient signal, $E_{\text{tr}}$, is similar at both 70
K and 100 K but the oscillating mode $E_{\text{osc}}$ is significantly weaker
at 100 K. To further understand the origin of the radiating mode, we studied
the dependence of $E_{\text{THz}}$ on the incident pump fluence, shown in
Figs. 3d-f. The spectral shape of $E_{\text{osc}}$ does not change with
decreasing pump fluence - only the overall spectral weight decreases linearly
with decreasing fluence (Fig. 3f), indicating that we are in a perturbative
regime.

Based on the observations above, we can rule out several possible sources for
the radiating mode at 0.23 THz, such as optical phonons [27, 23], longitudinal
acoustic phonons [28] or the CDW amplitudon. The lowest-frequency optical
phonons in (TaSe4)2I are at much higher frequencies (1.1 THz [29]) and THz
oscillatory emission from optical or acoustic phonons would be peaked at
$t_{\text{delay}}\sim 0$ ps followed by an exponential decay, unlike the
waveform shape we observe. The amplitudon in (TaSe4)2I cannot have a net IR-
dipole moment (due to the $D_{2}$ (222) point group symmetry of the
incommensurate phase [30]) which is necessary for the far-field coherent
radiation observed here. Furthermore, the spectral weight of THz emission due
to phonons or the amplitude mode is not expected to decrease to zero near
$T\sim$ 80 K. We also note, that the relatively low pump fluences used in our
experiments ($\sim$ 580 $\mu$Jcm-2) and the linear pump fluence dependence
precludes explanations in terms of a complete melting and recovery of CDW
order, which would have resulted in a saturating-type behavior with pump
fluence.

Figure 4: Model of a radiating massive Lee-Fukuyama phason coupled to a non-
radiating acoustic phonon. (a) $p(t)$ and $q(t)$ represent the time-profile of
the massive phason and acoustic phonon, respectively, as obtained from a
solution of a coupled harmonic oscillators model [21]. This particular set
corresponds to $p(t)$ obtained from a fit to the data at 7 K as shown in (b).
$N$ is a constant to rescale $p(t)$.

Having ruled out the possible sources discussed above, we posit that the
massive Lee-Fukuyama phason in (TaSe4)2I gives the coherent radiation at 0.23
THz. Within the LRA model for CDW dynamics [8, 1], the massive Lee-Fukuyama
phason frequency $\omega_{\text{LF}}$ at $\vec{q}=0$ is given as
$\omega_{\text{LF}}^{2}=\frac{3}{2}\lambda\omega_{\text{Q}}^{2}$, where
$\omega_{\text{Q}}$ is the frequency of the normal state longitudinal acoustic
phonon at $\vec{q}_{\text{CDW}}$ (Fig. 1b) and $\lambda$ is the electron-
phonon coupling constant. For (TaSe4)2I, $\omega_{\text{Q}}\approx 0.25$ THz,
as measured using inelastic neutron scattering [31, 32] while $\lambda\sim
0.6$ [33] for the strong electron-phonon coupling in (TaSe4)2I. This gives an
expected $\omega_{\text{LF}}\approx 0.24$ THz, close to our observed mode
frequency of 0.23 THz. Additionally, the temperature dependence of the 0.23
THz mode, which shows a considerable spectral weight only for $T\lesssim$ 0.3
$T_{\text{CDW}}$, agrees quite well with the random phase approximation (RPA)
calculations of Virosztek and Maki [11] for the temperature dependence of a
massive phason. Taken together, our observations strongly suggest that the
narrow-band 0.23 THz radiation originates from the predicted massive Lee-
Fukuyama phason in a CDW insulator.

To determine the excitation mechanism for the massive phason, we further
discuss the temporal waveform of the $E_{\text{osc}}$ highlighted earlier.
Note that the $E_{\text{osc}}$ shown in Fig. 2b is roughly zero at
$t_{\text{delay}}$ = 0 ps and then becomes significant only at later times. We
fit the behavior of $E_{\text{osc}}$ in time using a model of two coupled
modes, one of which is excited at $t_{\text{delay}}$ = 0 ps while only the
other one radiates as in Fig. 4 (See [21] for details on the fitting
procedure). This fitting model is consistent with a picture where the
radiating mode is the $\vec{q}=0$ massive phason with frequency
$\omega_{\mathrm{LF}}$, while the non-radiating mode is an acoustic phonon
mode with wavevector $\vec{q}=\vec{q}_{\text{CDW}}$ and frequency
$\omega_{\mathrm{Q}}\approx\omega_{\mathrm{LF}}$. The presence of such an
acoustic mode in the CDW phase was established in Refs. [32, 28]. Since the
acoustic phonon frequency $\omega_{\mathrm{Q}}\ll|\vec{q}_{\text{CDW}}|c$ is
well outside the light cone for free space radiation, the acoustic phonon
cannot directly radiate. However, since it has wavevector
$\vec{q}=\vec{q}_{\text{CDW}}$, it can scatter off the CDW modulation and
excite the massive phason mode. The massive phason mode at $\vec{q}=0$ will
then radiate, explaining the time profile of $\vec{E}_{\text{osc}}$ in Figs. 3
and 4. We hypothesize that the incident IR pump excites a large number of
single-particle excitations above the CDW gap. These excitations rapidly relax
via the generation of high-energy phonons which in turn eventually decay into
acoustic phonons. The low-lying acoustic phonon with
$\vec{q}=\vec{q}_{\text{CDW}}$ and
$\omega_{\mathrm{Q}}\approx\omega_{\mathrm{LF}}$ will then slowly decay via
scattering into the massive phason.

We note that (TaSe4)2I was recently argued to be an axion insulator (in the
presence of an external magnetic field) based on magnetoconductance
measurements of the sliding phason dynamics [19], though such an
interpretation is under active debate [34]. Since (TaSe4)2I is a good
insulator at low temperatures [21], the reported magnetoconductance
measurements [19] were limited to temperatures above 80 K. Our results imply
that at low temperature long-range Couloumb interaction is crucial for a full
understanding of the CDW order in (TaSe4)2I and related axion insulator
candidates. Contributions to the longitudinal magnetoconductance due to axion
electrodynamics may occur near the massive phason frequency, which can provide
a promising way to manipulate axionic states coherently.

To conclude, we emphasize that converting an ultrafast infrared pulse into a
multicycle THz pulse with a narrow bandwidth of a few GHz, as demonstrated
here, is in sharp contrast to existing ultrafast-based emitters that produce
single-cycle, broadband THz pulses. Our work thus reveals the possibility of
using fundamental collective modes as sources of robust and narrow-band THz
radiation. Such sources would enable resonant control of light-matter
interaction and narrow-band spectroscopy at the intrinsic energy scales of
correlated materials. Furthermore, our methodology here can provide a direct
dynamical probe of collective excitations and effects of long-range
interactions in materials with modulated order parameters such as spin or pair
density waves.

Acknowledgements: We thank Peter Abbamonte, Peter Armitage, Dipanjan
Chaudhuri, Tai-Chang Chiang, Lance Cooper, Steve Kivelson, Anshul Kogar,
Patrick Lee, Vidya Madhavan, Darius Torchinsky, and Benjamin Wieder for
insightful discussions.
This work was supported by the Quantum Sensing and Quantum Materials, an
Energy Frontier Research Center funded by the U.S. Department of Energy (DOE),
Office of Science, Basic Energy Sciences (BES), under Award No.DE-SC0021238.
F.M. acknowledges support from the EPiQS program of the Gordon and Betty Moore
Foundation, Grant GBMF11069. R.A.D. acknowledges support from the Bloch
Postdoctoral Fellowship in Quantum Science and Engineering of the Stanford
University Quantum Fundamentals, Architectures, and Machines initiative
(Q-FARM), and from the Karel Urbanek and Marvin Chodorow Postdoctoral
Fellowship of the Stanford University Department of Applied Physics. The
authors acknowledge the use of the spectroscopic ellipsometry setup at the
Institute for Basic Science (IBS) in Korea (Grant No. IBS-R009-D1).

Data availability. All relevant data are available on reasonable request.

Competing financial interests: The authors declare no competing financial
interests.

Author contributions: S.K., Y.L., N.B., A.M. and F.M. performed the THz
emission spectroscopy experiments and the corresponding data analysis. X-Q.S
and B.B. developed the theoretical understanding and modelling. C.Z., K.Q. and
D.P.S synthesized and characterized the samples. R.A.D., Q.L.D.N., and M.T.
gave crucial insights to the understanding and analysis of the data. S.K.,
X-Q.S, B.B. and F.M. wrote the manuscript with input from all the authors.
F.M. conceived and supervised this project.

## References

  * [1] Lee, P. A. & Fukuyama, H.  Dynamics of the charge-density wave. II. Long-range Coulomb effects in an array of chains.  _Physical Review B_ 17, 542–548 (1978).
  * [2] Schwartz, M. D.  _Quantum Field Theory and the Standard Model_ (Cambridge University Press, 2013).
  * [3] Fradkin, E.  _Field Theories of Condensed Matter Physics_ (Cambridge University Press, 2013).
  * [4] Coleman, P.  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015).
  * [5] Anderson, P. W.  Coherent excited states in the theory of superconductivity: Gauge invariance and the Meissner effect.  _Physical Review_ 110, 827 (1958).
  * [6] Anderson, P. W.  Random-phase approximation in the theory of superconductivity.  _Physical Review_ 112, 1900 (1958).
  * [7] Grüner, G.  The dynamics of charge-density waves.  _Reviews of Modern Physics_ 60, 1129–1181 (1988).
  * [8] Lee, P. A., Rice, T. M. & Anderson, P. W.  Fluctuation effects at a peierls transition.  _Physical Review Letters_ 31, 462–465 (1973).
  * [9] Lee, P., Rice, T. & Anderson, P.  Conductivity from charge or spin density waves.  _Solid State Communications_ 88, 1001–1007 (1993).
  * [10] Wong, K. Y. M. & Takada, S.  Effects of quasiparticle screening on collective modes: Incommensurate charge-density-wave systems.  _Physical Review B_ 36, 5476–5492 (1987).
  * [11] Virosztek, A. & Maki, K.  Collective modes in charge-density waves and long-range Coulomb interactions.  _Physical Review B_ 48, 1368–1372 (1993).
  * [12] Gressier, P., Guemas, L. & Meerschaut, A.  Preparation and structure of ditantalum iodide octaselenide, Ta2ISe8.  _Acta Crystallographica Section B Structural Crystallography and Crystal Chemistry_ 38, 2877–2879 (1982).
  * [13] Wang, Z. _et al._ Charge density wave transport in (TaSe4)2I.  _Solid State Communications_ 46, 325–328 (1983).
  * [14] Maki, M., Kaiser, M., Zettl, A. & Grüner, G.  Charge density wave transport in a novel inorganic chain compound (TaSe4)2I.  _Solid State Communications_ 46, 497–500 (1983).
  * [15] Fujishita, H., Sato, M. & Hoshino, S.  Incommensurate superlattice reflections in quasi one dimensional conductors (MSe4)2I (M=Ta and Nb).  _Solid State Communications_ 49, 313–316 (1984).
  * [16] Lee, K.-B., Davidov, D. & Heeger, A.  X-ray diffraction study of the CDW phase in (TaSe4)2I: Determination of the CDW modulation amplitude.  _Solid State Communications_ 54, 673–677 (1985).
  * [17] Donovan, S., Kim, Y., Alavi, B., Degiori, L. & Gruner, G.  The optical spectrum of charge density wave condensates.  _Solid State Communications_ 75, 721–724 (1990).
  * [18] Degiorgi, L. & Grüner, G.  Pinned and bound collective-mode state in charge-density-wave condensates.  _Physical Review B_ 44, 7820–7827 (1991).
  * [19] Gooth, J. _et al._ Axionic charge-density wave in the Weyl semimetal (TaSe4)2I.  _Nature_ 575, 315–319 (2019).
  * [20] Kim, S. _et al._ Kramers-Weyl fermions in the chiral charge density wave material (TaSe4)2I.  _arXiv:2108.10874_ (2021).
  * [21] See Supplementary Information .
  * [22] Belinicher, V. I. & Sturman, B. I.  The photogalvanic effect in media lacking a center of symmetry.  _Soviet Physics Uspekhi_ 23, 199–223 (1980).
  * [23] Sakai, K. (ed.) _Terahertz Optoelectronics_ (Springer Berlin Heidelberg, 2005).
  * [24] Rees, D. _et al._ Helicity-dependent photocurrents in the chiral Weyl semimetal RhSi.  _Science Advances_ 6, 29 (2020).
  * [25] Gao, Y. _et al._ Chiral terahertz wave emission from the Weyl semimetal TaAs.  _Nature Communications_ 11, 720 (2020).
  * [26] Luo, L. _et al._ A light-induced phononic symmetry switch and giant dissipationless topological photocurrent in ZrTe5.  _Nature Materials_ 20, 329–334 (2021).
  * [27] Dekorsy, T., Auer, H., Bakker, H. J., Roskos, H. G. & Kurz, H.  THz electromagnetic emission by coherent infrared-active phonons.  _Physical Review B_ 53, 4005–4014 (1996).
  * [28] Schaefer, H. _et al._ Dynamics of charge density wave order in the quasi one dimensional conductor (TaSe4)2I probed by femtosecond optical spectroscopy.  _The European Physical Journal Special Topics_ 222, 1005–1016 (2013).
  * [29] Lorenzo, J. E., Currat, R., Monceau, P., Hennion, B. & Levy, F.  Neutron investigation of optic-phonon branches in the quasi-one-dimensional compound (TaSe4)2I.  _Physical Review B_ 47, 10116–10121 (1993).
  * [30] Shi, W. _et al._ A charge-density-wave topological semimetal.  _Nature Physics_ 17, 381–387 (2021).
  * [31] Fujishita, H., Shapiro, S. M., Sato, M. & Hoshino, S.  A neutron scattering study of the quasi-one-dimensional conductor (TaSe4)2I.  _Journal of Physics C: Solid State Physics_ 19, 3049–3057 (1986).
  * [32] Lorenzo, J. E. _et al._ A neutron scattering study of the quasi-one-dimensional conductor.  _Journal of Physics: Condensed Matter_ 10, 5039–5068 (1998).
  * [33] Kim, T. W., Donovan, S., Grüner, G. & Philipp, A.  Charge-density-wave dynamics in (Ta1-xNbxSe4)2I alloys.  _Physical Review B_ 43, 6315–6325 (1991).
  * [34] Sinchenko, A. A., Ballou, R., Lorenzo, J. E., Grenet, T. & Monceau, P.  Does (TaSe4)2I really harbor an axionic charge density wave?  _Applied Physics Letters_ 120, 063102 (2022).