# Towards Interactively Improving ML Data Preparation Code via “Shadow
Pipelines”

Stefan Grafberger<EMAIL_ADDRESS>0000-0002-9884-9517 AIRLab, University
of AmsterdamThe Netherlands ,  Paul Groth<EMAIL_ADDRESS>0000-0003-0183-6910
University of AmsterdamThe Netherlands and  Sebastian Schelter  schelter@tu-
berlin.de 0000-0003-4722-5840 BIFOLD & TU BerlinGermany

###### Abstract.

Data scientists develop ML pipelines in an iterative manner: they repeatedly
screen a pipeline for potential issues, debug it, and then revise and improve
its code according to their findings. However, this manual process is tedious
and error-prone. Therefore, we propose to support data scientists during this
development cycle with automatically derived interactive suggestions for
pipeline improvements. We discuss our vision to generate these suggestions
with so-called shadow pipelines, hidden variants of the original pipeline that
modify it to auto-detect potential issues, try out modifications for
improvements, and suggest and explain these modifications to the user. We
envision to apply incremental view maintenance-based optimisations to ensure
low-latency computation and maintenance of the shadow pipelines. We conduct
preliminary experiments to showcase the feasibility of our envisioned approach
and the potential benefits of our proposed optimisations.

##  1\. Introduction

Software systems that learn from user data with machine learning (ML) have
become ubiquitous over the last years and participate in many critical
decision-making processes. Unfortunately, real-world experience shows that the
pipelines for data preparation, feature encoding, and model training in ML
systems are often brittle, especially with respect to issues in the data they
process (Breck et al., 2019; Polyzotis et al., 2018; Stoyanovich et al., 2022;
Schelter et al., 2018a). As a consequence, several data-centric techniques are
being developed to detect, quantify, and improve ML applications with respect
to reliability, fairness, and prediction quality (Jia et al., 2019; Schelter
et al., 2021; Li et al., 2019; Schelter et al., 2018b; Chung et al., 2019;
Tramer et al., 2017; Bird et al., 2020). Applying these techniques to ML
pipelines currently requires a high level of expertise, as existing approaches
such as mlinspect (Grafberger et al., 2022, 2021), DataScope (Karlaš et al.,
2023), mlwhatif (Grafberger et al., 2023a, b), Rain (Wu et al., 2020), Gopher
(Pradhan et al., 2022) or ArgusEyes (Schelter et al., 2022, 2023) assume that
data scientists know in advance what kind of errors they are looking for.

The need for interactively improving ML pipelines. In reality, data scientists
typically do not know in advance what pipeline issues to look for, and often
“discover serious issues only after deploying their systems in the real world”
(Holstein et al., 2019). At development time, data scientists currently have
to iteratively screen their pipeline for potential issues, debug these issues,
and then revise and improve the pipelines according to their findings. This
process is tedious, as it requires repeated manual code re-organisation and
re-execution in an environment like a Jupyter notebook.

We argue that ML pipeline development should be accompanied by interactive
suggestions to improve the pipeline code, similar to code inspections in
modern IDEs like IntelliJ (Jetbrains, [n. d.]) or text corrections in writing
assistants like Grammarly (Grammarly, [n. d.]). For that, we can re-use some
techniques from previous work, but are still faced with a

set of challenges: $(i)$ We need low-latency auto-detection of pipeline
improvement opportunities, to seamlessly integrate into the development
workflow; $(ii)$ we should identify pipeline problems spanning several
operators, instead of being artificially limited to screening individual
operators one-at-a-time as in e.g., mlinspect (Grafberger et al., 2022, 2021);
$(iii)$ users should receive provenance-enabled explanations (Chothia et al.,
2016) for detected problems and suggested improvements.

Figure 1. Our vision – several automatically maintained ”shadow pipelines”
give actionable suggestions on how to improve a user’s ML pipeline code at
development time.

Towards low-latency suggestions for pipeline improvement. We envision a system
that instruments a data scientist’s ML pipeline code and creates and maintains
so-called “shadow pipelines” with low-latency to generate suggestions for
improvements, as illustrated in Figure 1. Such a shadow pipeline is a hidden
variant of the original pipeline, which modifies it to auto-detect potential
issues and tries out different pipeline modifications for improvement
opportunities. Subsequently, each shadow pipeline provides the user with code
suggestions to improve the pipeline, accompanied by a provenance-based
explanation and a quantification of the expected impact on the pipeline
outputs. From a technical perspective, the main challenge is to conduct the
required computations with low latency by reusing and updating intermediates
via incremental view maintenance. We introduce the general problem using a
running example (Section 2), present our envisioned approach based on shadow
pipelines (Section 3), and finally conduct preliminary experiments for the
feasibility of our approach (Section 4).

In summary, this paper provides the following contributions:

  * •

We introduce the problem of interactively generating suggestions during ML
pipeline development (Section 2).

  * •

We propose shadow pipelines as an effective approach to computing such
suggestions and discuss ideas for efficiently maintaining them based on
incremental view maintenance (Section 3).

  * •

We conduct preliminary experiments with several shadow pipelines to show the
feasibility of computing low-latency suggestions, where we find that our
proposed optimisations potentially reduce the runtime by up to two orders of
magnitude (Section 4).

  * •

We provide the source code for our experiments at https://github.com/stefan-
grafberger/shadow-pipeline-experiments

##  2\. Problem Statement

Running example. Imagine a data scientist who is developing an ML pipeline for
a social media platform to detect user posts that show signs of the anhedonia
symptom for depression.

⬇

def load_data(countries_to_include):

users = pd.read_csv(’users.csv’)

posts = pd.read_parquet(’posts.parquet’)

users = users[users.country.isin(include_countries)]

return users.merge(posts, on=’user_id’)

# Weak labeling with expert rules

def signs_of_anhedonia(post):

return (post.contains(’(0|no|zero) (motivation)’)

| post.contains(’lost (interest|motivation)’))

& ~(posts.contains(’recover from (0|no|zero) interest’))

# Vectorstore retriever for RAG

def setup_retrieval_corpus(train_data, weak_labels):

vector_store = Chroma.from_texts(texts=train_data[’posts’],

metadatas=weak_labels, embedding=…))

return vectorstore.as_retriever()

# LLM + RAG pipeline using Langchain

def setup_llm_chain():

return ({”context”: retriever | format_docs,

”question”: RunnablePassthrough()}

| ChatPromptTemplate.from_template(…)

| MyCustomLLM()| JsonOutputParser() | format_response)

countries_to_include = [’US’, ’CAN’]

train_data = load_train_data(countries_to_include)

weak_labels = signs_of_anhedonia(train_data[’posts’])

retriever = setup_retrieval_corpus(train_data, weak_labels)

test_data = pd.read_parquet(’expert_labeled_data.parquet’)

rag_chain = setup_llm_chain(retriever)

y_predicted = rag_chain.batch(test[’posts’])

y_test_binarized = label_binarize(test[’signs_of_anhedonia’])

accuracy = accuracy_score(y_predicted, y_test_binarized)

Their pipeline (Listing LABEL:lst:example) loads and integrates data (Lines
1-5) and performs weak labeling of the training data, based on regular
expressions designed by experts (Lines 7-11). Next, the pipeline computes
embeddings of the weakly labeled samples, inserts them into a vector store
(Lines 13-17), and leverages these embeddings to classify samples using an LLM
with retrieval augmentation (Lines 19-24). In lines 26-35, the final pipeline
is then evaluated on a test set (which was manually labeled by experts). Note
that this example is inspired by an existing pipeline from NLP research
(Nguyen et al., 2022).

Iterative improvements for the pipeline. Next, imagine that the data scientist
finds that this initial pipeline version provides unsatisfactory accuracy and
starts to manually debug and revise the pipeline to figure out how to improve
it. They might, for example, write debugging code to look for mislabeled
training samples via Shapley Values (Karlaš et al., 2023) and detect that some
wrongly labeled samples are fed to the LLM from the retrieval corpus. After
mapping the embedded samples back to the original posts, the data scientist
would pinpoint the problem to the regular expressions used for weak labeling.
They would inspect the problematic posts and extend the regular expressions to
cover the cases that were not handled well. Next, the data scientist would re-
run the updated pipeline and observe that the accuracy improved a little but
that there are still instances of likely mislabeled samples. They again
isolate the samples to which the label error detection points and realise that
some users’ posts contain too many typos to be reliably matched. As a result,
the data scientist would once again rewrite the pipeline to integrate a
spelling correction step. Once no more likely mislabeled samples are detected,
the data scientist would try other evaluation techniques, e.g., to find small
slices in the test data where the pipeline does not perform well (Chung et
al., 2019; Sagadeeva and Boehm, 2021) and start another cycle of iterative
debugging and rewriting.

Towards automated, low-latency suggestions for improvements. As already
mentioned in the introduction, data scientists spend a lot of time in this
iterative process of “guessing” what a potential problem might be, rewriting
the pipeline to determine if that is indeed the case and then testing and
integrating improvements. Furthermore, this tedious process requires a high
level of expertise, as the data scientists need to identify the actual
improvement opportunities and concrete steps to fix the detected problems.

We argue that we should automate the underlying steps in the outlined
iterative improvement cycle as follows. First, common techniques for
identifying problematic portions of the data or the data preparation
operations (like label error detection or the detection of problematic data
slices) should be automatically run on the pipeline in the background. Once
potential issues are detected, data scientists should be presented with
provenance-based explanations for the locations in the data or pipeline code
that cause the issues. Furthermore, the data scientist should also be
presented with a list of suggestions on how to improve their pipeline, which
ideally already quantify the expected impact on a metric of interest such as
accuracy or fairness.

Existing systems such as mlwhatif (Grafberger et al., 2023b) or ArgusEyes
(Schelter et al., 2022) fall short to address this problem, as they assume
that the data scientist already knows exactly what issues to look for. They
are designed as one-off approaches with no consideration of the iterative
development cycle, and in general incur rather high execution times.

##  3\. Envisioned Approach

In the following, we describe our ideas for automatically maintained shadow
pipelines to assist the user with ML pipeline development.

Shadow pipelines. Figure 1 gives a high-level overview of our envisioned
approach, in reference to the running example from Section 2. On the left, we
see an automatically extracted “logical

plan” of a user’s ML pipeline, consisting of both relational and ML-specific
operators. We showed in previous work (Grafberger et al., 2022, 2023a) how to
extract such plans by instrumenting Python code. For the user pipeline, we
generate so-called shadow pipelines (shown on the right). A shadow pipeline is
a hidden variant of the original pipeline, which modifies it to auto-detect
potential issues and tries out different pipeline modifications to detect
improvement opportunities. Each shadow pipeline consists of several general
steps. The first step is issue detection, for which the shadow pipeline
introduces operators that take intermediates from the original pipeline as
input to screen for potential problems. For detected issues, a shadow pipeline
continues with root cause determination, aiming to localise specific pipeline
operators and/or input tuples responsible for the identified issue. Finally,
the shadow pipeline integrates and evaluates potential fixes, generates the
corresponding provenance-based change explanations, and provides a
quantification of the expected impact on the pipeline outputs. To prevent
information overload, shadow pipelines will have to use various techniques to
carefully select tuples for user explanations, e.g., via stratified sampling,
slice-finding (Chung et al., 2019; Sagadeeva and Boehm, 2021), based on
Shapley values (Karlaš et al., 2023; Jia et al., 2019), and by focusing on
test samples with label changes.

Low-latency computation of shadow pipelines. A key technical challenge is
ensuring low-latency shadow pipeline computations for an interactive user
experience. For that, we envision to reuse intermediates from the original
pipeline, sometimes even on a tuple-level, via incremental view maintenance
(IVM) techniques (McSherry et al., 2013). Every shadow pipeline must only
compute the minimum required difference compared to the original user pipeline
and should ideally avoid costly operations like model re-training. In cases
where such re-training is unavoidable, we plan to explore cheap proxy models
(Jia et al., 2019) to estimate the impact of a change on expensive models like
neural networks. Moreover, certain issue detection techniques and potential
fixes can be evaluated on selected subsets of the data only, when faced with
costly operations like LLM inference APIs or embedding computations for text
or image data. Some shadow pipelines may themselves involve expensive
operators like slice finding, Shapley value computation, spelling correction,
or text translation. Efficient implementations of these operators will be
crucial, and we may even need to work with samples only in extreme cases.
Furthermore, the parallel execution of different shadow pipelines and the
prioritisation of different analyses will be crucial for performance.

Maintenance of shadow pipelines. During the iterative development cycle of
data scientists, they continuously rewrite and re-run their pipelines. In
light of such rewrites, we again plan to update the original pipeline and its
shadow pipelines with IVM techniques, based on detected changes in operators
or input tuples. If the user followed a concrete code suggestion from our
shadow pipeline, we may even be able to update the original pipeline with
results previously computed in the shadow pipeline. However, in general,
updating ML pipelines presents some tough challenges, as IVM approaches
(McSherry et al., 2013) typically assume a fixed query with changing input
data, while in our scenario, both data and operators can change between
iterations. Therefore, we plan to apply a best-effort approach, where we only
apply IVM if certain simplifying assumptions are met, e.g. that only a single
operator changes.

Application to the running example. Consider the exemplary “label errors”
shadow pipeline shown in Figure 1, which uses Shapley values (Karlaš et al.,
2023; Jia et al., 2019) to identify likely label errors in the user pipeline’s
weakly labeled samples. Next, the shadow pipeline simulates fixes in the
regexes by flipping the labels of samples with negative Shapley values in the
train data, and re-evaluates the pipeline to observe the impact on its
accuracy. In the initial shadow pipeline run, we can directly apply the
Shapley value algorithm (Karlaš et al., 2023) on intermediates captured from
the original pipeline. Next, we can update the labels of these problematic
samples directly in the vector store metadata, and thereby reuse the
previously computed train embeddings, which do not change. Finally, the shadow
pipeline only needs to rerun the expensive inference for test set predictions
that relied on train samples with flipped labels. Similar optimisations can
quickly update the shadow pipeline with low latency in case of user changes,
and an incremental version of the Shapley value algorithm could even decrease
the latency further.

##  4\. Preliminary Experiments

We present preliminary experiments on computing and maintaining shadow
pipelines with hardcoded optimisations. Our goal is to validate the
feasibility of continuously managing these pipelines in the background as a
user works on their ML pipeline.

User pipelines. We experiment with two versions of the NLP pipeline from our
running example in Section 2. The first version, called rag, uses retrieval
augmentation and a trained LLM as detailed in Listing LABEL:lst:example. The
second version, called train, does not use an LLM but trains a simple neural
network classifier on the embeddings.

Shadow pipelines. We evaluate three shadow pipelines inspired by the issues
outlined in Section 2: slices uses SliceLine (Sagadeeva and Boehm, 2021) to
detect subsets of the test set where the user pipeline does not work well.
These samples correspond to posts in non-English languages; therefore, slices
integrates an additional translation step for posts from these languages as a
fix. The data-errors shadow pipeline determines how robust the user pipeline
is against data quality issues such as typos by evaluating synthetic errors,
and integrates an additional spell-checking step to improve the robustness.
Label-errors uses Shapley values (Jia et al., 2019; Karlaš et al., 2023) to
identify potentially mislabeled training samples from weak labeling,
subsequently flips the label of the most likely mislabeled samples and re-
evaluates the pipeline. All shadow pipelines also compute provenance-based
explanations.

Optimisations. For slices, SliceLine (Chung et al., 2019; Sagadeeva and Boehm,
2021) is applied directly on the original pipeline’s intermediates (parts of
the unfeaturised test data and the corresponding labels) to identify the most
problematic slice. In data-errors, we randomly corrupt 10% of the data with
typos, and recompute embeddings and run inference only for corrupted rows to
measure the impact. For potential fixes, we analogously restrict expensive
translation and spelling correction steps, embedding computation, and
inference to relevant test subsets. Label-errors also changes the train data.
For rag, we ensure that the expensive inference is re-run only for test set
predictions based on train samples with flipped labels. The train variant of
the user pipeline, unfortunately, requires neural network retraining; for
that, we evaluate a special optimisation referred to as opt-proxy, which
substitutes the neural network with a cheaper model and uses its performance
difference as a proxy for the label update’s impact. For subsequent shadow
pipeline updates after the user’s pipeline change, we re-use intermediates
from the original shadow pipeline run wherever possible. In particular, we
partially re-use translation, typo corruption, spell correction, and inference
results.

Data and APIs. We use a synthetic dataset with 1,000 rows, inspired by real
depression detection data (Nguyen et al., 2022), since the original dataset
cannot be shared due to its sensitive nature. For reproducibility, we fix
random seeds and replay cached answers from external APIs for the LLM and
translation step with artificial wait times. We repeat each experiment seven
times, preceded by one warm-up run, and report the median runtimes.

Figure 2. Benefits of our proposed optimisations for computing shadow
pipelines. Our shadow pipeline optimisations decrease the runtime by up to a
factor of 38.

Computing shadow pipelines. We measure the execution time for computing the
shadow pipelines after first executing the original pipeline. In each case, we
compare the runtime of the naive baseline (repeated execution of manually
rewritten versions of the original pipeline) to the runtime of the optimised
execution, where we implement hardcoded shadow pipelines with our proposed
optimisations.

Results and discussion. Figure 2 illustrates that the optimised shadow
pipelines decrease the runtime by up to a factor of 38 (55 with the proxy
model). The optimised variant consistently runs at least 4.7 times faster in
all scenarios. We find that the runtime reduction primarily depends on three
key factors: inference time, costly operations within shadow pipelines (e.g.,
spell correction, translation), and the need for retraining (not applicable
for rag). Substituting neural network retraining with a proxy model results in
a noticeable runtime improvement of 55, whereas the variant without a proxy
model only achieves a speedup of nine. Notably, for the scenarios with the
highest runtime (slices on train and data-errors on rag), the initial shadow
pipeline execution accounts for the majority of the runtime due to expensive
operators, where we can reuse intermediate results for shadow pipeline
maintenance. It is noteworthy that code inspections in modern IDEs like
IntelliJ can also incur several seconds of latency, which seems generally
acceptable for users, with the option to disable specific inspections.

Maintaining shadow pipelines. In the second experiment, we simulate that a
user makes a small change to their pipeline code (the update of a regex), and
then measure how long it takes to update the results of both the original
pipeline and its shadow pipelines.

Results and discussion. The incremental shadow pipeline maintenance is up to
626 times faster than the full naive re-execution, with original pipeline
runtimes decreasing up to fivefold. The additional re-use of intermediates
from previous shadow pipeline executions leads to bigger runtime decreases
than their initial optimised execution in the previous experiment (where we
only saw decreases by up to a factor of 55). Additionally, shadow pipeline
maintenance takes less than one second in all but one scenarios. The label
flipping for the most likely label errors in label-errors on rag affects train
rows retrieved for many test predictions, which leads to high inference costs,
which we cannot avoid. Notably, we now observe particularly low latencies for
the scenarios where we observed the highest runtimes in the previous
experiment (slices on train and data-errors on rag), as we can reuse expensive
intermediates from the initial execution of the respective shadow pipelines.

Figure 3. Benefits of our proposed optimisations for incrementally updating
the pipelines. The runtime (shown on a log scale) is less than one second in
all but one scenarios.

##  5\. Next Steps & Open Questions

We plan to explore how to create and efficiently maintain shadow pipelines for
user-defined pipelines, by automatically applying our outlined optimisations.
In particular, we aim to develop a runtime for incremental view maintenance
and shadow pipeline parallelisation (e.g., by building on DuckDB).
Additionally, we plan to conduct a user study to evaluate the usefulness of
our envisioned system and its provenance-based explanations for data
scientists. Several questions remain open, including whether to directly build
certain optimisations into the shadow pipelines or build a custom optimiser,
the extent to which we will need custom incremental operators, the reliance on
sampling, and approximations such as proxy models. Furthermore, we aim to
explore the integration of LLM-based code suggestions (GitHub, 2021) for
detected data-related problems into our system.

## References

  * (1)
  * Bird et al. (2020) Sarah Bird, Miro Dudík, Richard Edgar, Brandon Horn, Roman Lutz, Vanessa Milan, Mehrnoosh Sameki, Hanna Wallach, and Kathleen Walker. 2020.  Fairlearn: A toolkit for assessing and improving fairness in AI.  _Microsoft, Tech. Rep. MSR-TR-2020-32_ (2020).
  * Breck et al. (2019) Eric Breck, Neoklis Polyzotis, Sudip Roy, Steven Whang, and Martin Zinkevich. 2019.  Data Validation for Machine Learning.  _MLSys_ (2019).
  * Chothia et al. (2016) Zaheer Chothia, John Liagouris, Frank McSherry, and Timothy Roscoe. 2016.  _Explaining outputs in modern data analytics_.  Technical Report. ETH Zurich.
  * Chung et al. (2019) Yeounoh Chung, Tim Kraska, Neoklis Polyzotis, Ki Hyun Tae, and Steven Euijong Whang. 2019.  Slice finder: Automated data slicing for model validation. In _2019 IEEE 35th International Conference on Data Engineering (ICDE)_. IEEE, 1550–1553.
  * GitHub (2021) GitHub. 2021.  GitHub Copilot · Your AI pair programmer.  https://copilot.github.com/.
  * Grafberger et al. (2023a) Stefan Grafberger, Paul Groth, and Sebastian Schelter. 2023a.  Automating and Optimizing Data-Centric What-If Analyses on Native Machine Learning Pipelines.  _SIGMOD_ (2023).
  * Grafberger et al. (2022) Stefan Grafberger, Paul Groth, Julia Stoyanovich, and Sebastian Schelter. 2022.  Data distribution debugging in machine learning pipelines.  _VLDBJ_ (2022).
  * Grafberger et al. (2023b) Stefan Grafberger, Shubha Guha, Paul Groth, and Sebastian Schelter. 2023b.  mlwhatif: What If You Could Stop Re-Implementing Your Machine Learning Pipeline Analyses over and over?  _Proc. VLDB Endow._ 16, 12 (aug 2023), 4002–4005.  https://doi.org/10.14778/3611540.3611606
  * Grafberger et al. (2021) Stefan Grafberger, Shubha Guha, Julia Stoyanovich, and Sebastian Schelter. 2021.  MLINSPECT: A Data Distribution Debugger for Machine Learning Pipelines.  _SIGMOD_ (2021).
  * Grammarly ([n. d.]) Grammarly. [n. d.].  Demo.  https://demo.grammarly.com/.
  * Holstein et al. (2019) Kenneth Holstein, Jennifer Wortman Vaughan, Hal Daumé III, Miro Dudik, and Hanna Wallach. 2019.  Improving fairness in machine learning systems: What do industry practitioners need?. In _Proceedings of the 2019 CHI conference on human factors in computing systems_. 1–16.
  * Jetbrains ([n. d.]) Jetbrains. [n. d.].  Code inspections.  https://www.jetbrains.com/help/idea/code-inspection.html#access-inspections-and-settings.
  * Jia et al. (2019) Ruoxi Jia, David Dao, Boxin Wang, Frances Ann Hubis, Nezihe Merve Gurel, Bo Li, Ce Zhang, Costas J Spanos, and Dawn Song. 2019.  Efficient task-specific data valuation for nearest neighbor algorithms.  _PVLDB_ (2019).
  * Karlaš et al. (2023) Bojan Karlaš, David Dao, Matteo Interlandi, Sebastian Schelter, Wentao Wu, and Ce Zhang. 2023.  Data Debugging with Shapley Importance over Machine Learning Pipelines. In _The Twelfth International Conference on Learning Representations_.
  * Li et al. (2019) Peng Li, Xi Rao, Jennifer Blase, Yue Zhang, Xu Chu, and Ce Zhang. 2019.  Cleanml: A benchmark for joint data cleaning and machine learning.  _ICDE_ (2019).
  * McSherry et al. (2013) Frank McSherry, Derek Gordon Murray, Rebecca Isaacs, and Michael Isard. 2013.  Differential Dataflow.. In _CIDR_.
  * Nguyen et al. (2022) Thong Nguyen, Andrew Yates, Ayah Zirikly, Bart Desmet, and Arman Cohan. 2022.  Improving the Generalizability of Depression Detection by Leveraging Clinical Questionnaires. In _Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , Smaranda Muresan, Preslav Nakov, and Aline Villavicencio (Eds.). Association for Computational Linguistics, Dublin, Ireland, 8446–8459.  https://doi.org/10.18653/v1/2022.acl-long.578
  * Polyzotis et al. (2018) Neoklis Polyzotis, Sudip Roy, Steven Euijong Whang, and Martin Zinkevich. 2018.  Data lifecycle challenges in production machine learning: a survey.  _SIGMOD Record_ 47, 2 (2018).
  * Pradhan et al. (2022) Romila Pradhan, Jiongli Zhu, Boris Glavic, and Babak Salimi. 2022.  Interpretable data-based explanations for fairness debugging. In _Proceedings of the 2022 International Conference on Management of Data_. 247–261.
  * Sagadeeva and Boehm (2021) Svetlana Sagadeeva and Matthias Boehm. 2021.  Sliceline: Fast, linear-algebra-based slice finding for ml model debugging. In _Proceedings of the 2021 International Conference on Management of Data_. 2290–2299.
  * Schelter et al. (2018a) Sebastian Schelter, Felix Biessmann, Tim Januschowski, David Salinas, Stephan Seufert, and Gyuri Szarvas. 2018a.  On challenges in machine learning model management.  _IEEE Data Engineering Bulletin_ (2018).
  * Schelter et al. (2023) Sebastian Schelter, Stefan Grafberger, Shubha Guha, Bojan Karlaš, and Ce Zhang. 2023.  Proactively Screening Machine Learning Pipelines with ArgusEyes.  _SIGMOD_ (2023).
  * Schelter et al. (2022) Sebastian Schelter, Stefan Grafberger, Shubha Guha, Olivier Sprangers, Bojan Karlaš, and Ce Zhang. 2022.  Screening Native ML Pipelines with “ArgusEyes”.  _CIDR_ (2022).
  * Schelter et al. (2018b) Sebastian Schelter, Dustin Lange, Philipp Schmidt, Meltem Celikel, Felix Biessmann, and Andreas Grafberger. 2018b.  Automating Large-Scale Data Quality Verification.  _Proc. VLDB Endow._ 11, 12 (aug 2018), 1781–1794.  https://doi.org/10.14778/3229863.3229867
  * Schelter et al. (2021) Sebastian Schelter, Tammo Rukat, and Felix Biessmann. 2021.  JENGA - A Framework to Study the Impact of Data Errors on the Predictions of Machine Learning Models.  _EDBT_ (2021).
  * Stoyanovich et al. (2022) Julia Stoyanovich, Bill Howe, Serge Abiteboul, H.V. Jagadish, and Sebastian Schelter. 2022.  Responsible Data Management. In _Communications of the ACM_.
  * Tramer et al. (2017) Florian Tramer, Vaggelis Atlidakis, Roxana Geambasu, Daniel Hsu, Jean-Pierre Hubaux, Mathias Humbert, Ari Juels, and Huang Lin. 2017.  Fairtest: Discovering unwarranted associations in data-driven applications. In _2017 IEEE European Symposium on Security and Privacy (EuroS &P)_. IEEE, 401–416.
  * Wu et al. (2020) Weiyuan Wu, Lampros Flokas, Eugene Wu, and Jiannan Wang. 2020.  Complaint-driven training data debugging for query 2.0. In _Proceedings of the 2020 ACM SIGMOD International Conference on Management of Data_. 1317–1334.


billiard potential destroys the hyperbolic structure of the dispersive
billiard, and there are several known mechanisms for creating elliptic, KAM
non-degenerate periodic motions of a single particle in the billiard-like
Hamiltonian (21) at arbitrary small $\delta$ for the case of dispersive
$d(\geqslant 2)$-dimensional containers [51, 52, 36, 38]. It is conjectured
that such islands appear for dispersive billiard-like Hamiltonians generically
[51]. The phase space volume of the islands in these cases vanishes with
$\delta$, yet, its scaling with $\delta$ is known, and depends on the
asymptotic behavior of the container smoothing potential, $V$, near the
boundary. Our techniques imply that such stability islands of the single
particle motion can produce also choreographic solutions of the $N$ particle
system provided the perturbations induced by the interaction potential, $W$,
are much smaller than the islands size. Hence, the results regarding the
existence of KAM-stable choreographic solutions in containers that support
stable motion can probably be extended to any container.

Physical relevance of the coherent states. Under which conditions can the
constructed choreographic solutions be observed in realistic multi-particle
systems? General estimates on the probability for an initial condition to
belong to a KAM-torus of a Hamiltonian system with $N$-degrees of freedom are
quite pessimistic even for small $N$ [12, 13], but examples of such systems
where stability islands are well-noticeable are also known [2, 45], e.g. in
numerical experiments of [38] the islands are seen for $N$ as large as 20.
Therefore, Nekhoreshev-type estimates on the life time of the coherent states
are, probably, most relevant for the physical realizability question.
Sufficiently long living (i.e., effectively stable) coherent states can be of
direct interest when $N$ is not very large. For example, for a gas in a three-
dimensional rectangular box, the construction of $N<(l_{1}l_{2})/\rho^{2}$
particles moving vertically in synchrony, as in Theorem 5, corresponds to
pulsating fronts. For a gas in a convex container, the states of
$N<|L^{*}|/\rho$ particles constructed in Theorem 3 corresponds to rings of
current of a specific non-trivial spatial form. For each such state, the
effective stability imposes limitations on $N$ and the parameters of the
system. Additional classes of choreographic solutions may allow to study
coherent states with a larger number of particles.

Multi-path choreographic solutions provide such a class. When the single-
particle system is non-integrable and has an elliptic periodic orbit, it has,
typically, infinitely many elliptic orbits (e.g. around an elliptic periodic
orbit there are typically many resonant elliptic orbits, around which there
are secondary resonances, and so on [3]). Similarly, near a homoclinic loop to
a saddle center of a non-integrable system many stable periodic orbits co-
exist [29, 28]. When the corresponding periodic paths in the configuration
space do not intersect, one can obtain KAM-stable motions of repelling
particles along several such paths: the particles on the same path must have
the same frequency to avoid collisions, but the particles on different paths
may have different frequencies. When the paths do intersect, one needs the
frequencies of the motion along the different paths to be in resonance. This
condition is not actually restrictive: given any finite set of elliptic orbits
with arbitrary periods, one can tune the partial energies such that for every
two paths the ratio of the periods would become rational. For particles in the
billiard this is done just by normalizing the motion speed for each path by
the path length.

Choreographies in the box. The billiard in the box is integrable, and
therefore we build choreographic multi-particle regimes based on families of
parabolic periodic solutions (instead of elliptic orbits). We have considered
only one of such families in this paper, of parallel vertical motion, but
there are many types of them, for example, diamond shaped orbits. In our
construction of Section 5, to avoid collisions, a single particle sits on each
periodic path. However, for other types of parabolic families, such as diamond
shape orbits in a $d\geqslant 3$ box, many particles may occupy the same path
and parallel paths without any collisions. In particular, it may be possible
to create KAM-stable choreographic motions of $N\propto(l/\rho)^{3}$ particles
in a three-dimensional box of characteristic size $l$. Similar solutions can
be built for particles in ellipsoidal billiard, and for other systems of
weakly-interacting particles whose individual dynamics is integrable.

Dynamics of the averaged system. We have shown that the stable choreographic
motions are controlled by effective potentials defined on the torus
corresponding to the set of phases of the individual particles. Our main
result only refers to the fact that near the minimum of such potential the
dynamics are, generically, KAM stable. However, one may also ask a question of
the global dynamics: can the averaged system be completely integrable, or can
one find additional stable motions far from the minimum of the potential?
Various types of KAM-stable solutions of the averaged system should generate
new types of non-trivial coherent states which may depend differently on
physical parameters and be relevant for a larger variety of physical settings
than the dynamically simplest types of choreographies we found here. Also, as
explained in Appendix A, in the particular case of the equidistant particles’
phases, our effective potentials share the same symmetries as the potentials
of the classical Fermi-Pasta-Ulam chains, yet they form a larger class, and it
may be interesting to study this broader class of systems.

Solid coherent states in a high temperature gas. In this paper, we operate in
the limit where the motion along the periodic orbit is faster than the
oscillations of the phase differences between the particles. For highly
energetic particles in a container, if the non-averaged interaction potential
has a minimum (like the Lennard-Jones potential), one can think of an opposite
limit where the frequencies $\omega$ of the small oscillations of the
particles near this minimum are much faster than the frequency of the periodic
billiard-like motion of the center of mass: $\omega\gg\frac{\sqrt{2h}}{L^{*}}$
(the kinetic energy of the center-of-mass motion can still be much larger than
the energy of the fast but small oscillations of the particles: $I\ll
h/\omega$). It may be interesting to study whether this can also lead to
stable choreographic motions of molecules.

Large $N$ limit. As one can see, different types of coherent states can
correspond to different scaling of the number of particles as a function of
the size of the container. In general, when considering the limit
$N\to\infty$, one should also decide how the parameters of the system (the
system size $l$, the energy per particle $h$, effective particle diameter
$\rho$, etc.) scale with $N$: different types of scalings correspond to
different physical situations. Determining which type of scalings correspond
to various stable coherent states is the key for resolving the question of the
realizability of such states in physically relevant settings.

## Acknowledgements:

VRK work was funded by Israel Science Foundation (grant 787/22). DT work was
funded by the Leverhulme Trust (grant RPG-2021-072). We thank V. Gelfreich and
L. Chierchia for useful discussions.

## Appendix

##  Appendix A  Choreographic solutions with equidistant phases.

The evolution of the phases of the choreographic solutions is described, in
the first-order approximation, by the averaged Hamiltonian

$H=\frac{1}{2}(I_{0})^{2}+U(\theta),\quad(\theta,I_{0})\in\mathbb{T}^{N}\times\mathbb{R}^{N},$
(183)

see Lemma 3.1 (here we scale $I_{0}$ such that $a=1$ in the matrix $A$ of
(46), see (7)). The averaged potential $U$ is given by

$U(\theta^{(1)},\dots,\theta^{(N)})=\mathop{\sum_{n,m=1,\dots,N}}_{n\neq
m}W_{avg}(\theta^{(n)}-\theta^{(m)}),$

where $W_{avg}$ is an even, $2\pi$-periodic function defined by (11).

System (183) is a generalization of the classical Fermi-Pasta-Ulam chain: the
difference is that in the FPU we have only $m=n+1\mod N$ in the sum describing
the potential $U$, i.e., the interaction between the phases $\theta$ is short-
range, while in our case all phases typically interact with each other. Like
in the FPU, any uniformly distributed particle configuration
$\theta=\theta_{eq}$, where

$\theta_{eq}^{(2)}-\theta_{eq}^{(1)}=\ldots=\theta_{eq}^{(N)}-\theta_{eq}^{(N-1)}=\theta_{eq}^{(1)}+2\pi-\theta_{eq}^{(N)}=\frac{2\pi}{N}$

is an equilibrium of system (183). Indeed, it is always an extremum of $U$
since $\frac{\partial U}{\partial\theta^{(n)}}|_{\theta_{eq}}=0$ for each $n$:

$\frac{\partial U}{\partial\theta^{(n)}}|_{\theta_{eq}}=\sum_{m\neq
n}W_{avg}^{\prime}(\theta_{eq}^{(n)}-\theta_{eq}^{(m)})=\sum^{N-1}_{k=1}W_{avg}^{\prime}(\frac{2\pi}{N}k)=0,$

where the last equality follows because $W_{avg}^{\prime}$ is an odd and
$2\pi$-periodic function.

Such configurations form a line of extrema of $U$ (parameterized by the choice
of $\theta_{eq}^{(1)}$). It is a line of minima of $U$ when the Hessian matrix
at $\theta=\theta_{eq}$ is positive semi-definite. We have

$\frac{\partial^{2}U}{\partial\theta^{(n)}\partial\theta^{(m)}}(\theta_{eq})=\begin{cases}-W_{avg}^{\prime\prime}(\frac{2\pi(n-m)}{N})=u_{n-m}&\mbox{
for }\;n\neq m\\\
\sum^{N-1}_{j=1}W_{avg}^{\prime\prime}(\frac{2\pi}{N}j)=u_{0}=-\sum^{N-1}_{k=1}u_{k}&\mbox{
for }\;n=m.\\\ \end{cases}$ (184)

Note that the numbers $u_{k}$ satisfy $u_{k}=u_{-k}=u_{N-k}$ because
$W_{avg}^{\prime\prime}$ is even and $2\pi$-periodic.

The Hessian matrix (184) is a circulant matrix, so its eigenvectors are the
Fourier modes

$v_{j}=(1,\exp(i\frac{2\pi j}{N}),\exp(i\frac{4\pi
j}{N}),\ldots,\exp(i\frac{2(N-1)j\pi}{N})),\quad j=0,\ldots,N-1.$

The corresponding eigenvalues are:

$\lambda_{j}=u_{0}+u_{1}\exp(i\frac{2\pi j}{N})+u_{2}\exp(i\frac{4\pi
j}{N})+\ldots+u_{N-2}\exp(i\frac{2\pi j(N-2)}{N})+u_{N-1}\exp(i\frac{2\pi
j(N-1)}{N}),$

i.e.,

$\lambda_{j}=2\sum_{1\leqslant k<N/2}u_{k}(\cos(k\frac{2\pi
j}{N})-1)+((-1)^{j}-1)u_{N/2}$ (185)

(we take $u_{N/2}=0$ in this formula when $N$ is odd).

As we see, $\lambda_{0}=0$, which is due to the translational symmetry of the
Hamiltonian (183). When all other $\lambda_{j}$ are strictly positive, the
line of uniformly distributed particles’ configurations consists of minima of
the potential $U$. This happens, for example, when $u_{k}<0$ for all $k\neq
0$, i.e., when $W_{avg}$ is a convex function, which is consistent with the
repelling nature of the interaction.

The frequencies of small oscillations around the line of minima are equal to
$\sqrt{\lambda_{j}}$. It follows from (185) that

$\lambda_{j}=\lambda_{N-j}$

for all $1\leqslant j\leqslant N/2$. Hence, the standard non-resonance
assumption on the interaction potential breaks at $\theta=\theta_{eq}$. In
particular, adding a small perturbation to (183) (without breaking the
translational symmetry) would result, in general, in the destruction of the
ellipticity.

Yet, we show next that the KAM Assumption IP1 still holds generically. Notice
that the resonance relations are due to the discrete symmetries of the
potential near $\theta=\theta_{eq}$. Namely, the system in a small
neighborhood of the minima line is symmetric with respect to the
transformations $T:\theta^{(n)}\rightarrow\theta^{(n+1)}$ (where $n$ is taken
$\mod N$) and $S:\theta^{(n)}\rightarrow\theta^{(N+1-n)}$. The maps $S$ and
$T$ generate the so-called $N$-th dihedral group. This is the same group of
symmetries as in the FPU. The normal form theory for the FPU was built by Rink
[39]. In fact, he derived the normal form near an equilibrium of a general
$S,T$-symmetric Hamiltonian provided the Hamiltonian has no additional
resonances. His work applies to our Hamiltonian (even though the interactions
here are for all particle pairs whereas in [39] the FPU chain with only
nearest neighbors interactions was considered). Indeed, notice that for any
prescribed sequence $u_{k},k=1,\dots,N-1$ there is a smooth even and
$2\pi$-periodic potential $W_{avg}$ satisfying
$W_{avg}^{\prime\prime}(\frac{2\pi}{N}k)=u_{k}$, so, generically, for our
system, no additional independent resonance relations appear.

The symmetric Rink normal form is given by Theorem 8.2 in [39] and for odd $N$
is written as follows:

$H_{Rink}=\sum_{1\leqslant
j<\frac{N}{2}}\sqrt{\lambda_{j}}\;a_{j}+\frac{1}{2}\sum_{1\leqslant
j,k<\frac{N}{2}}(C^{a}_{jk}a_{j}a_{k}+C^{b}_{jk}b_{j}b_{k})$ (186)

whereas for even $N$

$H_{Rink}=\sum_{1\leqslant
j\leqslant\frac{N}{2}}\sqrt{\lambda_{j}}\;a_{j}+\frac{1}{2}\sum_{1\leqslant
j,k\leqslant\frac{N}{2}}C^{a}_{jk}a_{j}a_{k}+\frac{1}{2}\sum_{1\leqslant
j,k<\frac{N}{2}}C^{b}_{jk}b_{j}b_{k}+\sum_{1\leqslant
j\leqslant\frac{N}{4}}C^{d}_{j}(d_{j}d_{\frac{N}{2}-j}-c_{j}c_{\frac{N}{2}-j}).$
(187)

The terms $a,b,c,d$ are quadratic functions, so the normal form is of order
$4$ (the $C$’s are constant coefficients such that
$C^{*}_{jk}=C^{*}_{kj},*\in\\{a,b,c,d\\}$). The transition to the normal form
is done as follows. First, one makes a linear symplectic coordinate
transformation which diagonalizes the quadratic part of the Hamiltonian near
$(\theta=\theta_{eq},I_{0}=0)$. One defines

$\begin{array}[]{l}\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{for
}\;1\leqslant j<N/2:\\\ \displaystyle\qquad
z_{j}=\frac{1}{\sqrt{2N}}\sum_{n=1}^{N}e^{-2\pi
i\frac{j}{N}n}(I_{0}^{(n)}+i\sqrt{\lambda_{j}}\;(\theta^{(n)}-\theta_{eq}^{(n)})),\\\
\displaystyle\qquad\zeta_{j}=\frac{1}{\sqrt{2N\lambda_{j}}}\sum_{n=1}^{N}e^{2\pi
i\frac{j}{N}n}(I_{0}^{(n)}-i\sqrt{\lambda_{j}}\;(\theta^{(n)}-\theta_{eq}^{(n)})),\\\
\displaystyle\qquad z_{N-j}=-\;\frac{1}{\sqrt{2N}}\sum_{n=1}^{N}e^{-2\pi
i\frac{j}{N}n}(I_{0}^{(n)}-i\sqrt{\lambda_{j}}\;(\theta^{(n)}-\theta_{eq}^{(n)})),\\\
\displaystyle\qquad\zeta_{N-j}=\frac{1}{\sqrt{2N\lambda_{j}}}\sum_{n=1}^{N}e^{2\pi
i\frac{j}{N}n}(I_{0}^{(n)}+i\sqrt{\lambda_{j}}\;(\theta^{(n)}-\theta_{eq}^{(n)})),\\\
\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{and,
if $N$ is even}:\\\ \displaystyle\qquad
z_{\frac{N}{2}}=\frac{1}{\sqrt{2N\lambda_{\frac{N}{2}}}}\sum_{n=1}^{N}(-1)^{n}(I_{0}^{(n)}+i\sqrt{\lambda_{j}}\;(\theta^{(n)}-\theta_{eq}^{(n)})),\\\
\displaystyle\qquad\zeta_{\frac{N}{2}}=-\;\frac{1}{\sqrt{2N}}\sum_{n=1}^{N}(-1)^{n}(I_{0}^{(n)}-i\sqrt{\lambda_{j}}\;(\theta^{(n)}-\theta_{eq}^{(n)})),\end{array}$

cf. [39], formulas (7.2),(7.3). By the translational invariance of the average
potential $U$, the Hamiltonian (183) in these coordinates is the sum of the
term $\frac{1}{2N}\sum_{n=1}^{N}(I^{(n)}_{0})^{2}$ and a function which
depends only on $(z_{j},\zeta_{j})$ with $j=1,\ldots,N-1$. This function is
the Hamiltonian of the system reduced by the translation symmetry group
($\theta^{(n)}\mapsto\theta^{(n)}+c,\;c\in\mathbb{R}^{1},n=1,\ldots,N$). Since
$(z,\zeta)=0$ is an equilibrium of the reduced system, the Taylor expansion of
the reduced Hamiltonian at zero starts with quadratic terms – these terms
coincide with the quadratic part of $H_{Rink}$ in (186) and (187), i.e., with
$\sum_{1\leqslant j\leqslant\frac{N}{2}}\sqrt{\lambda_{j}}\;a_{j}$. Next, one
does a symplectic transformation (identity plus terms of the second order and
higher) of the variables $(z_{j},\zeta_{j})$, $j=1,\ldots,N-1$, which brings
the reduced Hamiltonian to the form which coincides with the normal forms
(186) or (187) up to terms of order $5$. It is a standard fact that this
system has KAM-tori if the normal form has KAM-tori in an arbitrarily small
neighborhood of zero, so we further focus on finding KAM-tori in the normal
forms.

In the new variables $(z,\zeta)$ the functions $a,b,c,d$ in $H_{Rink}$ are
given by

$a_{j}=i(z_{j}\zeta_{j}-z_{N-j}\zeta_{N-j}),\qquad
b_{j}=i(z_{j}\zeta_{j}+z_{N-j}\zeta_{N-j}),$
$c_{j}=\frac{1}{\sqrt{\lambda_{j}}}(z_{j}z_{N-j}+\lambda_{j}\zeta_{j}\zeta_{N-j}),\qquad
d_{j}=\frac{-i}{\sqrt{\lambda_{j}}}(z_{j}z_{N-j}-\lambda_{j}\zeta_{j}\zeta_{N-j}),$

for $1\leqslant j<\frac{N}{2}$; we also have, for even $N$,

$a_{\frac{N}{2}}=iz_{\frac{N}{2}}\zeta_{\frac{N}{2}}.$

Note that $a,b,c,d$ are real and satisfy

$a_{j}^{2}=b_{j}^{2}+c_{j}^{2}+d_{j}^{2},$ (188)

see [39], formulas (8.6),(8.7). The coefficients $C_{jk}^{*}$ in $H_{Rink}$
are polynomials of the coefficients of the original Hamiltonian depending
rationally on $\sqrt{\lambda_{j}}$. One can check that a generic potential $U$
corresponds to a generic choice of $C$’s.

For odd $N$, the normal form (186) is completely integrable and generically
has KAM tori [39]. This gives us that the the KAM assumption IP1 is fulfilled
generically for system (183) for odd $N$. Let us consider the case of even
$N$. By [39], Corollary 9.3, the $(N-1)$ degrees of freedom system (187) has a
number of quadratic integrals in involution: $a_{j}$, $j=1,\ldots,\frac{N}{2}$
and $b_{j}-b_{\frac{N}{2}-j}$, $1\leqslant j<\frac{N}{4}$. This set is
incomplete; while the normal form $H_{Rink}$ for the nearest-neighbor FPU
chain with an even $N$ has additional integrals and is completely integrable
[19, 20], it is not known whether this normal form is completely integrable
for a generic choice of coefficients $C$ and $\lambda$ in (187), or for a
general choice of the potential $U$ in (183).

However, the restriction of (187) to the invariant subspace
$\\{a_{s}=0,1\leqslant s\leqslant\frac{N}{4}\\}$ is completely integrable:
since $b_{s}$, $c_{s}$ and $d_{s}$ all vanish for $1\leqslant
s\leqslant\frac{N}{4}$ by (188), the restricted Hamiltonian is given by

$H=\sum_{\frac{N}{4}<j\leqslant\frac{N}{2}}\sqrt{\lambda_{j}}\;a_{j}+\frac{1}{2}\sum_{\frac{N}{4}<j,k\leqslant\frac{N}{2}}C^{a}_{jk}a_{j}a_{k}+\frac{1}{2}\sum_{\frac{N}{4}<j,k<\frac{N}{2}}C^{b}_{jk}b_{j}b_{k},$

and the quadratic functions $a_{j}$ ($\frac{N}{4}<j\leqslant\frac{N}{2}$) and
$b_{j}$ ($\frac{N}{4}<j<\frac{N}{2}$) give a complete set of its integrals.
The restricted Hamiltonian has the same structure as the completely integrable
normal form (186). So, in the same way as it was done in [39] for system
(186), one establishes that for a generic choice of non-zero values of the
integrals $a_{j}$ ($\frac{N}{4}<j\leqslant\frac{N}{2}$) and $b_{j}$
($\frac{N}{4}<j<\frac{N}{2}$) the corresponding joint level set of these
integrals is a KAM-torus $\mathcal{T}$ of the restricted system. Moreover,
these integrals $a_{j}$ and $b_{j}$ are the action variables (see formula
(9.3) in [39]). The dynamics in a small neighborhood of such torus are
described, in the main approximation, by the Hamiltonian (187) averaged over
the angle variables conjugate to the actions. The only terms in (187) that
depend on these angle variables are $d_{\frac{N}{2}-j}$ and
$c_{\frac{N}{2}-j}$ with $j<\frac{N}{4}$; let us show that their averaged
values are zero.

Indeed, by the ergodicity of the flow on the invariant torus $\mathcal{T}$, we
can replace the averaging over the angle variables by the time averaging. Due
to the commutation relations (see formula (9.1) in [39])

$\\{b_{k},c_{m}\\}=2d_{m}\delta_{km},\qquad\\{b_{k},d_{m}\\}=-2d_{m},\qquad\\{a_{k},c_{m}\\}=\\{a_{k},d_{m}\\}=0,$

where $\\{\cdot,\cdot\\}$ is the Poisson bracket and $\delta_{km}$ is the
Kronecker delta, we have that for $\frac{N}{4}<m<\frac{N}{2}$

$\frac{d}{dt}c_{m}=\\{c_{m},H_{Rink}\\}=-\Omega_{m}\;d_{m},\qquad\frac{d}{dt}d_{m}=\\{d_{m},H_{Rink}\\}=\Omega_{m}\;c_{m},$

where

$\Omega_{m}=2\sum_{\frac{N}{4}<k<\frac{N}{2}}C^{b}_{km}b_{k}.$

For a generic choice of the actions $b_{k}$ ($\frac{N}{4}<k<\frac{N}{2}$), the
frequencies $\Omega_{m}$ are all non-zero, hence $c_{m}$ and $d_{m}$ perform
harmonic oscillations and their time-average is zero. Thus, the averaged
system (187) is

$H=\sum_{1\leqslant
j\leqslant\frac{N}{2}}\sqrt{\lambda_{j}}\;a_{j}+\frac{1}{2}\sum_{1\leqslant
j,k\leqslant\frac{N}{2}}C^{a}_{jk}a_{j}a_{k}+\frac{1}{2}\sum_{1\leqslant
j,k<\frac{N}{2}}C^{b}_{jk}b_{j}b_{k}.$

This system has the same structure as (186) and is completely integrable (the
integrals are $a_{j}$ and $b_{j}$). It follows that the invariant torus
$\mathcal{T}$ is normally elliptic. As in [39], one checks that a generic
Liouville torus of this systems satisfies the twist condition, hence
$\mathcal{T}$ is surrounded by KAM-tori and these tori persist when we proceed
from the averaged system to the original system $H_{Rink}$ in a neighborhood
of $\mathcal{T}$.

## References

  * [1] R. Abraham and J. E. Marsden.  Foundations of mechanics.  Number 364. American Mathematical Soc., 2008.
  * [2] C. Antonopoulos, T. Bountis, and C. Skokos.  Chaotic dynamics of n-degree of freedom hamiltonian systems.  International Journal of Bifurcation and Chaos, 16(06):1777–1793, 2006.
  * [3] V. Arnold, V. Kozlov, and A. Neishtadt.  Mathematical Aspects of Classical and Celestial Mechanics, volume 3.  Springer Science & Business Media, 2007.
  * [4] P. Bálint, T. Gilbert, D. Szász, and I. P. Tóth.  What mathematical billiards teach us about statistical physics?  Pure and Applied Functional Analysis, 6(1):1–35, 2021.
  * [5] L. Barreira and Y. Pesin.  Smooth ergodic theory and nonuniformly hyperbolic dynamics.  In Handbook of dynamical systems. Vol. 1B, pages 57–263. Elsevier B. V., Amsterdam, 2006.  With an appendix by Omri Sarig.
  * [6] L. Bunimovich.  Kinematics, equilibrium, and shape in Hamiltonian systems: the “LAB” effect.  Chaos. An Interdisciplinary Journal of Nonlinear Science, 13(3):903–912, 2003.
  * [7] L. Bunimovich, C. Liverani, A. Pellegrinotti, and Y. Suhov.  Ergodic systems of $n$ balls in a billiard table.  Communications in Mathematical Physics, 146(2):357–396, 1992.
  * [8] L. A. Bunimovich, Y. G. Sinai, and N. I. Chernov.  Statistical properties of two-dimensional hyperbolic billiards.  Russian Mathematical Surveys, 46(4):47–106, aug 1991.
  * [9] A. Chenciner, J. Gerver, R. Montgomery, and C. Simó.  Simple choreographic motions of $N$ bodies: a preliminary study.  In Geometry, mechanics, and dynamics, pages 287–308. Springer, New York, 2002.
  * [10] A. Chenciner and R. Montgomery.  A remarkable periodic solution of the three-body problem in the case of equal masses.  Annals of Mathematics. Second Series, 152(3):881–901, 2000.
  * [11] N. Chernov and R. Markarian.  Chaotic Billiards.  Number 127. American Mathematical Soc., 2006.
  * [12] L. Chierchia and G. Pinzari.  Properly-degenerate kam theory (following vi arnold).  Discrete & Continuous Dynamical Systems-S, 3(4):545, 2010.
  * [13] L. Chierchia and G. Pinzari.  The planetary n-body problem: symplectic foliation, reductions and invariant tori.  Inventiones mathematicae, 186(1):1–77, 2011.
  * [14] M. M. Dvorin and V. F. Lazutkin.  Existence of an infinite number of elliptic and hyperbolic periodic trajectories for convex billiards.  Akademija Nauk SSSR. Funkcionalnyi Analiz i ego Prilozenija, 7(2):20–27, 1973.
  * [15] J. Fejoz, A. Knauf, and R. Montgomery.  Classical $n$-body scattering with long-range potentials.  Nonlinearity, 34(11):8017–8054, 2021.
  * [16] H. Fukuda, T. Fujiwara, and H. Ozaki.  Figure-eight choreographies of the equal mass three-body problem with lennard-jones-type potentials.  Journal of Physics A: Mathematical and Theoretical, 50(10):105202, feb 2017.
  * [17] G. Fusco, G. F. Gronchi, and P. Negrini.  Platonic polyhedra, topological constraints and periodic solutions of the classical n-body problem.  Inventiones mathematicae, 185(2):283–332, 2011.
  * [18] M. Guardia, P. Martin, J. Paradela, and T. M. Seara.  Hyperbolic dynamics and oscillatory motions in the 3 body problem, 2022\.
  * [19] A. Henrici and T. Kappeler.  Results on normal forms for FPU chains.  Communications in Mathematical Physics, 278(1):145–177, 2008.
  * [20] A. Henrici and T. Kappeler.  Resonant normal form for even periodic FPU chains.  Journal of the European Mathematical Society (JEMS), 11(5):1025–1056, 2009.
  * [21] A. Katok and B. Hasselblatt.  Introduction to the modern theory of dynamical systems, volume 54 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, 1995.  With a supplementary chapter by Katok and Leonardo Mendoza.
  * [22] A. Knauf.  The $n$-centre problem of celestial mechanics for large energies.  Journal of the European Mathematical Society (JEMS), 4(1):1–114, 2002.
  * [23] V. Kozlov and D. Treschev.  A genetic introduction to the dynamics of systems with impacts.  AMS, Providence, 1991.
  * [24] V. V. Kozlov.  Billiards, invariant measures, and equilibrium thermodynamics.  Regular & Chaotic Dynamics. International Scientific Journal, 5(2):129–138, 2000.
  * [25] V. V. Kozlov.  Polynomial conservation laws for the Lorentz and the Boltzmann-Gibbs gases.  Uspekhi Matematicheskikh Nauk, 71(2(428)):81–120, 2016.
  * [26] N. S. Krylov.  Works on the foundations of statistical physics.  Princeton Series in Physics. Princeton University Press, Princeton, N.J., 1979.  Translated from the Russian by A. B. Migdal, Ya. G. Sinai [Ja. G. Sinaĭ] and Yu. L. Zeeman [Ju. L. Zeeman], With a preface by A. S. Wightman, With a biography of Krylov by V. A. Fock [V. A. Fok], With an introductory article “The views of N. S. Krylov on the foundations of statistical physics” by Migdal and Fok, With a supplementary article “Development of Krylov’s ideas” by Sinaĭ.
  * [27] S. Lansel, M. A. Porter, and L. A. Bunimovich.  One-particle and few-particle billiards.  Chaos. An Interdisciplinary Journal of Nonlinear Science, 16(1):013129, 12, 2006.
  * [28] L. Lerman.  Hamiltonian systems with loops of a separatrix of a saddle-center.  Selecta Math. Sov, 10(3):297–306, 1991.
  * [29] L. Lerman and K. Trifonov.  Saddle-center and periodic orbit: Dynamics near symmetric heteroclinic connection.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(2):023113, 2021.
  * [30] K. Meyer, G. Hall, and D. Offin.  Introduction to Hamiltonian Dynamical Systems and the N-Body Problem, volume 90.  Springer Science & Business Media, 2008.
  * [31] R. Montgomery.  N-body choreographies.  Scholarpedia, 5(11):10666, 2010.  revision #137285.
  * [32] R. Montgomery.  The three-body problem and the shape sphere.  American Mathematical Monthly, 122(4):299–321, 2015.
  * [33] C. Moore.  Braids in classical dynamics.  Physical Review Letters, 70(24):3675–3679, 1993.
  * [34] A. Neishtadt, K. Sheng, and V. Sidorenko.  Stability analysis of apsidal alignment in double-averaged restricted elliptic three-body problem.  Celestial Mechanics and Dynamical Astronomy, 133(10):1–23, 2021\.
  * [35] T. Ouyang and Z. Xie.  Star pentagon and many stable choreographic solutions of the newtonian 4-body problem.  Physica D: Nonlinear Phenomena, 307:61–76, 2015.
  * [36] A. Rapoport and V. Rom-Kedar.  Nonergodicity of the motion in three-dimensional steep repelling dispersing potentials.  Chaos, 16(4):043108, 6, 2006.
  * [37] A. Rapoport, V. Rom-Kedar, and D. Turaev.  Approximating multi-dimensional Hamiltonian flows by billiards.  Comm. Math. Phys., 272(3):567–600, 2007.
  * [38] A. Rapoport, V. Rom-Kedar, and D. Turaev.  Stability in high dimensional steep repelling potentials.  Comm. Math. Phys., 279(2):497–534, 2008.
  * [39] B. Rink.  Symmetry and resonance in periodic fpu chains.  Communications in Mathematical Physics, 218(3):665–685, 2001.
  * [40] V. Rom-Kedar and D. Turaev.  Billiards: A singular perturbation limit of smooth Hamiltonian flows.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(2):026102, 2012.
  * [41] N. Simányi.  The complete hyperbolicity of cylindric billiards.  Ergodic Theory and Dynamical Systems, 22(1):281–302, 2002.
  * [42] N. Simányi.  Proof of the Boltzmann-Sinai ergodic hypothesis for typical hard disk systems.  Inventiones Mathematicae, 154(1):123–178, 2003.
  * [43] N. Simányi.  Proof of the ergodic hypothesis for typical hard ball systems.  Annales Henri Poincaré. A Journal of Theoretical and Mathematical Physics, 5(2):203–233, 2004.
  * [44] N. Simányi and D. Szász.  Hard ball systems are completely hyperbolic.  Annals of Mathematics. Second Series, 149(1):35–96, 1999.
  * [45] C. Simó.  Dynamical properties of hamiltonian systems with applications to celestial mechanics.  In Central Configurations, Periodic Orbits, and Hamiltonian Systems, pages 169–232. Springer, 2015.
  * [46] Y. G. Sinai.  On the foundations of the ergodic hypothesis for a dynamical system of statistical mechanics.  In Doklady Akademii Nauk, volume 153, pages 1261–1264. Russian Academy of Sciences, 1963.
  * [47] Y. G. Sinai.  Dynamical systems with elastic reflections.  Russian Mathematical Surveys, 25(2):137, 1970.
  * [48] Y. G. Sinai.  Hyperbolicity and chaos.  In Boltzmann’s legacy 150 years after his birth (Rome, 1994), volume 131 of Atti Convegni Lincei, pages 107–110. Accad. Naz. Lincei, Rome, 1997.
  * [49] J. G. Sinaĭ.  Dynamical systems with elastic reflections. Ergodic properties of dispersing billiards.  Akademiya Nauk SSSR i Moskovskoe Matematicheskoe Obshchestvo. Uspekhi Matematicheskikh Nauk, 25(2 (152)):141–192, 1970.
  * [50] D. Turaev.  Richness of chaos in the absolute newhouse domain.  In Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II–IV: Invited Lectures, pages 1804–1815. World Scientific, 2010.
  * [51] D. Turaev and V. Rom-Kedar.  Elliptic islands appearing in near-ergodic flows.  Nonlinearity, 11(3):575, 1998.
  * [52] D. Turaev and V. Rom-Kedar.  Soft billiards with corners.  J. Statist. Phys., 112(3-4):765–813, 2003.
  * [53] M. P. Wojtkowski.  Design of hyperbolic billiards.  Communications in Mathematical Physics, 273(2):283–304, 2007.
  * [54] G. Yu.  Simple choreographies of the planar newtonian n-body problem.  Archive for Rational Mechanics and Analysis, 225(2):901–935, 2017\.