11institutetext: University of Hamburg, Department of Informatics
11email<EMAIL_ADDRESS>22institutetext: University Medical Center Hamburg-Eppendorf
22email<EMAIL_ADDRESS>
# Segmenting Medical Instruments in Minimally Invasive Surgeries using
AttentionMask

Christian Wilms 11 Alexander Michael Gerlach 11 Rüdiger Schmitz 22
Simone Frintrop 11

###### Abstract

Precisely locating and segmenting medical instruments in images of minimally
invasive surgeries, medical instrument segmentation, is an essential first
step for several tasks in medical image processing. However, image
degradations, small instruments, and the generalization between different
surgery types make medical instrument segmentation challenging. To cope with
these challenges, we adapt the object proposal generation system AttentionMask
and propose a dedicated post-processing to select promising proposals. The
results on the ROBUST-MIS Challenge 2019 show that our adapted AttentionMask
system is a strong foundation for generating state-of-the-art performance. Our
evaluation in an object proposal generation framework shows that our adapted
AttentionMask system is robust to image degradations, generalizes well to
unseen types of surgeries, and copes well with small instruments.

##  1 Introduction

The localization and pixel-precise segmentation of medical instruments in
images of minimally invasive surgeries, medical instrument segmentation, is
important for numerous medical applications. For instance, it allows for skill
assessment [15], procedure identification, or automated camera steering [29],
among other tasks [22]. However, localizing and segmenting medical instruments
is limited by several application-specific characteristics [19, 22]. Images
from minimally invasive surgeries are often degraded by motion blur, blood,
smoke, or poor illumination [22]. Some of these challenging conditions are
visible in the example images in Fig. 1. Additionally, minimally invasive
surgeries are applied in various body parts, instruments, and surrounding
tissues vary. Hence, approaches for medical instrument segmentation have to be
robust and generalizable to support minimally invasive surgeries in clinical
practice.

|  |
---|---|---
|  |
Figure 1: Example images from minimally invasive surgeries with medical
instruments. The challenges induced by motion blur (left column), smoke
(central column), and poor illumination (top right) are well visible. All
images are taken from the ROBUST-MIS Challenge 2019 test set [22, 17].

Multiple formulations for medical instrument segmentation exist [2, 22]. These
include generation of a common binary mask for all instruments [14, 6],
segmentation of individual instances [12, 7, 22, 4] or instrument parts [23,
19]. We focus on the segmentation of individual instances without
classification. Various computer vision approaches based on semantic
segmentation [23, 19, 12] and instance segmentation [7, 22, 4] were proposed
for the different formulations. Besides these kinds of approaches, object
proposal generation methods [20, 21, 11, 24, 25, 26] can also be part of the
armamentarium when focusing on the segmentation of individual instances.
Object proposal generation methods produce a set of class-agnostic object
candidates as boxes or pixel-precise masks. These candidates are typically
used to limit the search space in object detection [16] or instance
segmentation [8]. Since the systems are class-agnostic, they generalize better
to unseen configurations than instance segmentation methods [20, 18, 13].

In this paper, we adapt the object proposal generation system AttentionMask
[24] to medical instrument segmentation in the context of the ROBUST-MIS
Challenge 2019 [22]. The ROBUST-MIS Challenge 2019 aims to assess the state-
of-the-art in medical instrument segmentation and foster more research on this
task by introducing a dataset and an evaluation framework. Utilizing
AttentionMask for medical instrument segmentation is very promising since it
has shown a strong performance in discovering objects of various sizes [24]
and demonstrated good results in other complex real-world applications [27,
28]. To adapt AttentionMask, we modify the training to fit the new data and
propose a dedicated post-processing to select the most relevant object
proposals. The results show that our adapted AttentionMask is a promising
foundation for generating high-quality medical instrument segmentations. We
further investigate the influence of the instruments’ size as well as the
robustness and generalization ability of our system.

##  2 Related Work

This section gives a brief overview of approaches to medical instrument
segmentation and object proposal generation.

###  2.1 Medical Instrument Segmentation

As mentioned above, the different formulations of medical instrument
segmentation lead to various approaches. Tackling the problem of generating
one binary segmentation mask for all instruments, [14] and [6] utilize
encoder-decoder architectures. [23] and [19] apply semantic segmentation
networks to locate and classify different instrument parts.

For segmenting individual instances, the task we address in this paper,
several authors use instance segmentation systems. [4] propose CCAM as a
combination of YOLACT++ [3] with an attention mechanism. Most ROBUST-MIS
Challenge 2019 participants that segment individual instances employ the
instance segmentation system Mask R-CNN [8]. Mostly, few details like the
anchor boxes or the training process are adapted to the new data. Similarly,
[7] use Mask R-CNN with temporal consistency constraints to utilize available
video data. An exception to these instance segmentation-based approaches are
[12]. They first solve the binary segmentation task and extract instances
using connected component analysis.

In contrast to these approaches, we adapt an object proposal generation system
showing strong generalization and robustness abilities.

###  2.2 Object Proposal Generation

Object proposal generation is the class-agnostic discovery of objects with
boxes or masks (proposals) [1]. For each proposal, an objectness score is
estimated to rank the proposals. Besides methods using hand-crafted features
[10], CNN-based approaches were proposed [20, 21, 11, 24, 25, 26]. [20] and
[21] generate foreground-background segmentations and objectness scores for
several overlapping crops from an image pyramid. In contrast, FastMask [11]
generates a feature pyramid inside a CNN for a more efficient processing.
AttentionMask [24] further improves the efficiency by utilizing attention to
only extract windows in relevant parts of the feature pyramid. The improved
efficiency allows a better discovery of small objects utilizing a new module.
We will discuss AttentionMask in more detail in Sec. 3.1 since we adapt the
system to segment medical instruments.

##  3 Method

We adapt the object proposal system AttentionMask [24] (see Fig. 2) to medical
instrument segmentation in images of minimally invasive surgeries. As
discussed above, AttentionMask showed a strong performance in localizing and
segmenting objects [24] as well as a good robustness w.r.t. image degradations
[27]. To adapt AttentionMask, we update the training to fit the medical
instrument segmentation application and integrate a post-processing that
selects the most promising object proposals in three steps. In the following,
we first briefly review AttentionMask and present our updated training in Sec.
3.1. Subsequently, we introduce our dedicated post-processing in Sec. 3.2.

###  3.1 AttentionMask

Figure 2: System figure of AttentionMask. The input image is processed by the
first part of the backbone network (left yellow module), yielding the feature
pyramid’s base level ($\mathcal{S}_{8}$). Further levels are generated by
downsampling up to level $\mathcal{S}_{128}$. At each pyramid level, windows
are extracted in regions with high scale-specific attention (see red areas in
heatmaps) to omit background areas. A pixel-precise mask and an objectness
score are generated for each extracted window. Note that we cropped the input
image for clarity. The figure is based on [24].

AttentionMask locates and segments medical instruments by first extracting a
feature map from the input image using the first part of a ResNet-based
backbone [9] (left yellow module in Fig. 2). The result is a feature map
downsampled by a factor of 8 w.r.t. the input image. This pyramid level
($\mathcal{S}_{8}$) serves as the base for the subsequent feature pyramid,
which is generated using the remainder of the backbone and residual pooling
layers (orange module in Fig. 2) [11]. Overall, the feature pyramid consists
of 8 pyramid levels representing the image downsampled by a factor of 8 for
small instruments to 128 for large instruments.

Across the feature pyramid, several windows are extracted to locate and
segment the instruments. To only extract windows in relevant areas,
AttentionMask utilizes a trained attention mechanism (scale specific
objectness attention in Fig. 2). This attention highlights areas with
instruments that fit the pyramid levels in terms of size. Subsequently,
$10\times 10$ windows are extracted from the feature maps in high attention
areas (green windows in Fig. 2). Note that extracting fixed-size windows from
the different pyramid levels allows AttentionMask to locate and segment
instruments of various sizes. Given the extracted windows, AttentionMask
generates segmentation masks for the instruments and objectness scores to rank
the generated instrument proposals.

To train the system, we adapt the training regime used in [24] by reducing the
learning rate to 0.00001 to fit the medical instrument segmentation data and
train the system for only 12 epochs.

###  3.2 Post-Processing

We apply three post-processing steps to select the most promising from the
initial instrument proposals. This step is necessary since the dense window
sampling in high attention areas within AttentionMask leads to several similar
proposals covering the same object. Although this is common in object proposal
generation, the evaluation frameworks typically used in medical instrument
segmentation penalize redundant proposals.

The first post-processing step removes instrument proposals with an objectness
score below 0.8 to revoke proposals that most likely cover only object parts
or background patches. Second, we combine groups of at least five partially
overlapping proposals and discard all other proposals. This step counters the
effects of the dense proposal generation in AttentionMask, leading to several
proposals covering the same instrument with slightly different segmentations.
From each group of proposals, we generate one final proposal by selecting all
pixels contained in at least 10% of the group’s proposals. These three steps
remove 98.2% of the proposals from the original AttentionMask output to better
fit the evaluation frameworks commonly utilized in medical instrument
segmentation. Note that we do not apply this post-processing when evaluating
our approach in the object proposal generation framework.

##  4 Evaluation

To assess the quality of our adapted AttentionMask system, we evaluate on the
dataset of the ROBUST-MIS Challenge 2019 [22, 17]. The dataset contains 10040
annotated images from 30 minimally invasive surgeries. We divide the official
training set consisting of 5983 images into a training split (80%) and a
validation split (20%). The remaining 4057 images are used as the test set
similar to [22]. This test set was subdivided by [22] into three stages to
test the generalization ability. While the first stage comprises images of
known patients and surgeries, the second stage includes data of known surgery
types but different patients. Finally, the third and most challenging stage
consists of images from an unseen surgery type. Hence, the second and third
stages of the test data are independent of the training data. We follow [22]
and report the stage-specific results.

For measuring the quality of the results, we use two strategies. First, we
follow the ROBUST-MIS Challenge 2019 [22] and calculate the Multi-instance
Dice Similarity Coefficient (MI_DSC) and the Multi-instance Normalized Surface
Dice (MI_NSD). While the MI_DSC assesses the segmentations’ quality by
comparing entire regions with the annotation, MI_NSD only uses the areas
around the contours. To better assess the consistency of the performance and
the applicability in clinical practice, we focus on the worst-case results and
report the 5% quantile for MI_DSC and MI_NSD as in [22].

As a second evaluation strategy, we use a standard object proposal evaluation
protocol [20, 21, 11, 24, 25, 26]. Hence, we evaluate the results using
Average Recall (AR) [10] for different numbers of proposals. AR assesses how
many annotated objects are discovered and how well they are segmented. We also
calculate the recall for annotated objects of different size categories,
similar to [20, 21, 11, 24, 25, 26]. The recall defines how many annotated
objects are located given an Intersection over Union (IoU) value for assessing
the overlap between the proposal and the annotated object. We define our size
categories in terms of relative size (object area divided by image area). The
categories are XS ($<1\%$), S($1\%-2\%$), M ($2\%-5\%$), L ($5\%-10\%$), and
XL ($>10\%$).

We compare our results to the participants of the ROBUST-MIS Challenge 2019
caresyntax, CASIA_SRL, fisensee [12], SQUASH, Uniandes, VIE, and www based on
the results presented in [22]. Additionally, we compare to CCAM [4] in the
evaluation framework of [22]. Note that we evaluate different versions of our
proposed approach. In the ROBUST-MIS Challenge 2019 evaluation framework (Sec.
4.1), we evaluate our adapted AttentionMask with and without the dedicated
post-processing as well as the adapted AttentionMask with an optimal ranking.
The latter version gives an upper limit for our approach and selects the best
fitting proposal per ground truth object. Hence, it assumes an optimal
ranking. In the object proposal generation framework (Sec. 4.2 and Sec. 4.3),
we only evaluate our adapted AttentionMask without the dedicated post-
processing to match the evaluation framework.

###  4.1 Quantitative Results in the ROBUST-MIS Framework

First, we present the results using MI_DSC and MI_NSD on the ROBUST-MIS
Challenge 2019 test set stage 3 in Tab. 1. This test set only features a
surgery type not seen in training. The results indicate that our dedicated
post-processing improves the results of our adapted AttentionMask system, but
only outperforms two and three other systems, respectively. The best systems
in terms of MI_DSC outperform our system with post-processing by 0.17 due to
the impaired initial ranking. If we assume an optimal ranking for our system,
the results substantially improve, outperforming all systems by at least 67.8%
and 74.3% in terms of MI_DSC and MI_NSD. These results show the potential of
our adapted AttentionMask system for medical instrument segmentations, but
also clearly indicate the low quality of the initial ranking.

Note that a result of 0.00 does not imply that a method misses all
instruments, since we only focus on the worst-case scenarios and report 5%
quantiles [22].

Table 1: Quantitative results of the three versions of our adapted AttentionMask, CCAM [4], and the participants of the ROBUST-MIS Challenge 2019 [22] in terms of Multi-instance Dice Similarity Coefficient (MI_DSC) and Multi-instance Normalized Surface Dice (MI_NSD) on stage 3 of the ROBUST-MIS Challenge 2019 test set [22, 17]. All results are 5% quantiles. Bold font highlights the best current results, while italic font denotes the best results given an optimal ranking for our system to showcase its potential. Participant/System | MI_DSC$\uparrow$ | MI_NSD$\uparrow$
---|---|---
VIE | 0.00 | 0.00
caresyntax | 0.00 | 0.00
fisensee [12] | 0.17 | 0.16
CASIA_SRL | 0.19 | 0.27
SQUASH | 0.22 | 0.26
Uniandes | 0.26 | 0.29
www | 0.31 | 0.35
CCAM [4] | 0.31 | 0.34
Ours (without post-processing) | 0.00 | 0.00
Ours (with post-processing) | 0.14 | 0.19
Ours (with optimal ranking) | 0.52 | 0.61

###  4.2 Quantitative Object Proposal Generation Results

In the second part of our evaluation, we analyze our adapted AttentionMask
without the dedicated post-processing in more detail using Average Recall (AR)
measures. Removing the post-processing matches the evaluation framework of
object proposal generation. Table 2 presents the overall results in terms of
AR for the three stages of the ROBUST-MIS Challenge 2019 test set. The results
indicate that the first 10 proposals locate most instruments that our adapted
AttentionMask discovers. Between AR@1 and AR@10, the results improve by 0.280
across all stages, while the improvement from AR@10 to AR@100 is only 0.063 on
average. Hence, the rough ranking of proposals in AttentionMask has a high
quality, while the precise ranking of the most relevant proposals is
suboptimal. Comparing the results across the different stages, no major
difference is visible. For instance, in terms of AR@10 the difference between
the best result (stage 2) and the worst result (stage 3) is only 0.082.
Therefore, AttentionMask shows strong generalization abilities across
different patients (stage 2) and surgery types (stage 3).

Table 2: Quantitative results of our adapted AttentionMask without post-processing in terms of Average Recall (AR) using different numbers of proposals on the three stages of the ROBUST-MIS Challenge 2019 test set [22, 17]. Test split stage | AR@1$\uparrow$ | AR@10$\uparrow$ | AR@100$\uparrow$
---|---|---|---
Stage 1 | 0.182 | 0.471 | 0.533
Stage 2 | 0.214 | 0.502 | 0.554
Stage 3 | 0.156 | 0.420 | 0.497

We also evaluate the quality of our instrument segmentations w.r.t. the
instruments’ size. Figure 3 shows the results per stage and size category in
terms of the located instruments (recall) for different Intersection over
Union (IoU) levels. The results indicate that except for very small
instruments (XS), the results for the different size categories (S-XL) are
very similar. However, even for very small instruments, more than 80% of the
instruments are located and segmented properly using the common IoU threshold
of 0.5 [5]. These results are in line with the findings in [24], showing a
strong performance of AttentionMask on small objects. Only for IoU levels
above 0.7 the results across all sizes start to drop substantially. Overall,
the results of AttentionMask are only weakly related to the instruments’ size,
which is different from the participants of the ROBUST-MIS Challenge 2019
according to [22].

(a) Stage 1 (b) Stage 2 (c) Stage 3

Figure 3: Results of our adapted AttentionMask without post-processing in
terms of recall on the three stages of the ROBUST-MIS Challenge 2019 test set
[22, 17]. The results per stage are divided by the relative size of the
instruments. The categories are XS ($<1\%$), S($1\%-2\%$), M ($2\%-5\%$), L
($5\%-10\%$), and XL ($>10\%$).

###  4.3 Qualitiative Results

Figure 4 presents qualitative results of our adapted AttentionMask on four
images from the challenging stage 3 of the ROBUST-MIS Challenge 2019 test set.
The first three examples show high-quality instrument segmentations, despite
different challenges like multiple objects, smoke, or suboptimal illumination.
These results indicate a strong robustness of AttentionMask and are in line
with the findings in [22] and [27]. However, the final example shows a
challenging scenario where AttentionMask fails to discover two of three
instruments. We attribute this error to the strong occlusion of the green and
blue instruments (mask color in ground truth). Another typical limitation is
visible in the second example. Due to the inherited downsampling process in
CNNs, AttentionMask misses small details of the instruments, like parts of the
tip in this example. This behavior is similar to the findings of [25, 26]
applying AttentionMask on typical computer vision datasets.

|  |
---|---|---
|  |
|  |
|  |
Input image | Ours | Ground truth
Figure 4: Qualitative results of our adapted AttentionMask without post-
processing on four images from the ROBUST-MIS Challenge 2019 test set [22, 17]
covering different challenging conditions. Filled colored contours denote
found objects, while not filled red contours denote missed objects. Note that
only the best fitting proposal (highest IoU) is visualized per annotated
object.

##  5 Conclusion

This paper addressed medical instrument segmentation in minimally invasive
surgery images. We adapted the object proposal generation system AttentionMask
and proposed a dedicated post-processing pipeline to utilize AttentionMask for
medical instrument segmentation. Results on the ROBUST-MIS Challenge 2019
dataset showed that our adapted AttentionMask system is a promising foundation
for precise and robust medical instrument segmentation. Specifically, we
showed robustness of our system w.r.t. image degradations like smoke or
suboptimal illumination. Importantly, AttentionMask-based medical instrument
segmentation generalizes well to unseen scenarios and localizes instruments of
almost all sizes with high quality. The latter behavior is especially
intriguing since it differs from previous systems [22].

A main limitation of our system lies in the non-optimal ranking of the
segmentation mask proposals. As some of these proposals do not accurately
reflect the detailed shape of the corresponding object, non-optimal ranking
results in a considerable performance degradation in terms of segmentation
metrics that focus exactly on the quality of the segmented object outline. We
have shown that falsely high ranking of inaccurate proposals is the central
hurdle for our system to outperform state-of-the-art systems. Superpixel-based
refinement [25, 26] may alleviate this problem through reduction of ill-shaped
object proposals.

## References

  * [1] Alexe, B., Deselaers, T., Ferrari, V.: What is an object? In: Computer Vision and Pattern Recognition (2010)
  * [2] Allan, M., Shvets, A., Kurmann, T., Zhang, Z., Duggal, R., Su, Y.H., Rieke, N., Laina, I., Kalavakonda, N., Bodenstedt, S., et al.: 2017 robotic instrument segmentation challenge. arXiv preprint arXiv:1902.06426 (2019)
  * [3] Bolya, D., Zhou, C., Xiao, F., Lee, Y.J.: YOLACT++: Better real-time instance segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence (2020)
  * [4] Ceron, J.C.A., Chang, L., Ochoa-Ruiz, G., Ali, S.: Assessing YOLACT++ for real time and robust instance segmentation of medical instruments in endoscopic procedures. In: International Conference of the IEEE Engineering in Medicine and Biology Society (2021)
  * [5] Everingham, M., Van Gool, L., Williams, C.K.I., Winn, J., Zisserman, A.: The PASCAL Visual Object Classes (voc) challenge. International Journal of Computer Vision 88(2) (2010)
  * [6] García-Peraza-Herrera, L.C., Li, W., Fidon, L., Gruijthuijsen, C., Devreker, A., Attilakos, G., Deprest, J., Vander Poorten, E., Stoyanov, D., Vercauteren, T., Ourselin, S.: ToolNet: holistically-nested real-time segmentation of robotic surgical tools. In: International Conference on Intelligent Robots and Systems (2017)
  * [7] González, C., Bravo-Sánchez, L., Arbelaez, P.: ISINet: An instance-based approach for surgical instrument segmentation. In: Medical Image Computing and Computer Assisted Intervention (2020)
  * [8] He, K., Gkioxari, G., Dollár, P., Girshick, R.: Mask R-CNN. In: International Conference on Computer Vision (2017)
  * [9] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Computer Vision and Pattern Recognition (2016)
  * [10] Hosang, J., Benenson, R., Dollár, P., Schiele, B.: What makes for effective detection proposals? IEEE Transactions on Pattern Analysis and Machine Intelligence 38(4) (2015)
  * [11] Hu, H., Lan, S., Jiang, Y., Cao, Z., Sha, F.: FastMask: Segment multi-scale object candidates in one shot. In: Computer Vision and Pattern Recognition (2017)
  * [12] Isensee, F., Maier-Hein, K.H.: OR-UNet: an optimized robust residual U-Net for instrument segmentation in endoscopic images. arXiv preprint arXiv:2004.12668 (2020)
  * [13] Kim, D., Lin, T.Y., Angelova, A., Kweon, I.S., Kuo, W.: Learning open-world object proposals without learning to classify. arXiv preprint arXiv:2108.06753 (2021)
  * [14] Laina, I., Rieke, N., Rupprecht, C., Vizcaíno, J.P., Eslami, A., Tombari, F., Navab, N.: Concurrent segmentation and localization for tracking of surgical instruments. In: Medical Image Computing and Computer Assisted Intervention (2017)
  * [15] Lin, S., Qin, F., Bly, R.A., Moe, K.S., Hannaford, B.: Automatic sinus surgery skill assessment based on instrument segmentation and tracking in endoscopic video. In: MICCAI International Workshop on Multiscale Multimodal Medical Imaging (2019)
  * [16] Mai, X., Zhang, H., Meng, M.Q.H.: Faster R-CNN with classifier fusion for small fruit detection. In: International Conference on Robotics and Automation (2018)
  * [17] Maier-Hein, L., Wagner, M., Ross, T., Reinke, A., Bodenstedt, S., Full, P.M., Hempe, H., Mindroc-Filimon, D., Scholz, P., Tran, T.N., et al.: Heidelberg colorectal data set for surgical data science in the sensor operating room. Scientific Data 8(1) (2021)
  * [18] Ošep, A., Voigtlaender, P., Weber, M., Luiten, J., Leibe, B.: 4D generic video object proposals. In: International Conference on Robotics and Automation (2020)
  * [19] Pakhomov, D., Premachandran, V., Allan, M., Azizian, M., Navab, N.: Deep residual learning for instrument segmentation in robotic surgery. In: MICCAI International Workshop on Multiscale Multimodal Medical Imaging (2019)
  * [20] Pinheiro, P.O., Collobert, R., Dollár, P.: Learning to segment object candidates. In: Advances in Neural Information Processing Systems (2015)
  * [21] Pinheiro, P.O., Lin, T.Y., Collobert, R., Dollár, P.: Learning to refine object segments. In: European Conference on Computer Vision (2016)
  * [22] Roß, T., Reinke, A., Full, P.M., Wagner, M., Kenngott, H., Apitz, M., Hempe, H., Mindroc-Filimon, D., Scholz, P., Tran, T.N., et al.: Comparative validation of multi-instance instrument segmentation in endoscopy: Results of the ROBUST-MIS 2019 challenge. Medical Image Analysis 70 (2021)
  * [23] Shvets, A.A., Rakhlin, A., Kalinin, A.A., Iglovikov, V.I.: Automatic instrument segmentation in robot-assisted surgery using deep learning. In: International Conference on Machine Learning and Applications (2018)
  * [24] Wilms, C., Frintrop, S.: AttentionMask: Attentive, efficient object proposal generation focusing on small objects. In: Asian Conference on Computer Vision (2018)
  * [25] Wilms, C., Frintrop, S.: Superpixel-based refinement for object proposal generation. In: International Conference on Pattern Recognition (2020)
  * [26] Wilms, C., Frintrop, S.: DeepFH segmentations for superpixel-based object proposal refinement. Image and Vision Computing 114 (2021)
  * [27] Wilms, C., Heid, R., Sadeghi, M.A., Ribbrock, A., Frintrop, S.: Which airline is this? Airline logo detection in real-world weather conditions. In: International Conference on Pattern Recognition (2020)
  * [28] Wilms, C., Johanson, R., Frintrop, S.: Localizing small apples in complex apple orchard environments. arXiv preprint arXiv:2202.11372 (2022)
  * [29] Zhang, J., Gao, X.: Object extraction via deep learning-based marker-free tracking framework of surgical instruments for laparoscope-holder robots. International Journal of Computer Assisted Radiology and Surgery 15(8) (2020)


# Ergodic decay laws in Newtonian and relativistic chaotic scattering

Diego S. Fernández<EMAIL_ADDRESS>Álvaro G. López  Jesús M. Seoane
Miguel A. F. Sanjuán  Nonlinear Dynamics, Chaos and Complex Systems Group,
Departamento de Física, Universidad Rey Juan Carlos, Tulipán s/n, 28933
Móstoles, Madrid, Spain

###### Abstract

In open Hamiltonian systems, the escape from a bounded region of phase space
according to an exponential decay law is frequently associated with the
existence of hyperbolic dynamics in such a region. Furthermore, exponential
decay laws based on the ergodic hypothesis are used to describe escapes in
these systems. However, we uncover that the presence of the set that governs
the hyperbolic dynamics, commonly known as the chaotic saddle, invalidates the
assumption of ergodicity. For the paradigmatic Hénon-Heiles system, we use
both theoretical and numerical arguments to show that the escaping dynamics is
non-ergodic independently of the existence of KAM tori, since the chaotic
saddle, in whose vicinity trajectories are more likely to spend a finite
amount of time evolving before escaping forever, is not utterly spread over
the energy shell. Taking this into consideration, we provide a clarifying
discussion about ergodicity in open Hamiltonian systems and explore the
limitations of ergodic decay laws when describing escapes in this kind of
systems. Finally, we generalize our claims by deriving a new decay law in the
relativistic regime for an inertial and a non-inertial reference frames under
the assumption of ergodicity, and suggest another approach to the description
of escape laws in open Hamiltonian systems.

###### pacs:

05.45.Ac,05.45.Df,05.45.Pq

##  I Introduction

When a point mass travels freely and meets a region of interaction, commonly
described in terms of a potential, its state of motion is irremediably
affected. The particle may escape eventually from the influence of the
potential, resuming its free journey until another obstacle is encountered.
This physical phenomenon is known as a scattering process, and if the
potential function is nonlinear, the particle might experience a chaotic
scattering seoane2013 . In this case, the particle will be subject to chaotic
dynamics that will cause its final state of motion to be extremely sensitive
to modifications in the initial conditions, which hinders its predictability
aguirre2009 . In most scattering processes, the particle experiences the
effect of chaotic dynamics just for a finite amount of time before escaping,
i.e., its motion undergoes a transient chaotic dynamics lai2010 ; tel2015 .
This phenomenon is typical in open 2D time-independent Hamiltonian systems,
among which the Hénon-Heiles system constitutes a paradigmatic example
aguirre2001 . Briefly, transient chaotic dynamics is governed by the presence
in phase space of the chaotic saddle, also called invariant nonattracting
chaotic set ott1993 , which rules how much time particles spend in the phase
space region called the scattering region. Chaotic scattering is a fundamental
field in nonlinear dynamics, which has many different applications in physics,
such as molecule scattering, advection of particles in fluid mechanics,
transition of materials, or even stars escaping from galaxies lin2013 ;
daitche2014 ; toledomarin2018 ; navarro2019 .

To study the probabilistic laws that govern the particle escape from the
influence of the potential, large ensembles of particles can be launched
inside the scattering region, where they interact with the potential well. In
this way, particles escape following an exponential decay law when hyperbolic
chaos governs the dynamics, whereas escapes take place according to algebraic
laws when the underlying dynamics is nonhyperbolic bauer1990 ; alt1996 ;
lau1991 ; kokshenev2000 ; bunimovich2005 . On the other hand, exponential
decay laws are sometimes associated with ergodic motions, so it is common to
find works in which the concepts of chaos and ergodicity are used in an
equivalent way zhao2007 . Boltzmann coined the terms _ergoden_ or _ergodische_
, whose Greek etymological origin means “energy path” moore2015 , to refer to
systems in which trajectories, if left to itself for long time enough, will
pass close to nearly all the dynamical states on the energy surface
lebowitz1973 , i.e., the dynamical states compatible with a constant energy.
It is widely known that in small systems with a few degrees of freedom chaos
is decisive to determine if the laws of statistical mechanics are satisfied,
and therefore it is the key to many proofs that demonstrate that systems are
ergodic berdichevsky1988 .

However, exponential decay laws based on the ergodic hypothesis, which we
shall name _ergodic decay laws_ hereafter for simplicity, are inadequate to
describe the rich phenomenology of open Hamiltonian systems, as confirmed by
the numerical results presented in this work. The reason is that assuming this
hypothesis in the open regime is an oversimplifying approximation for small
exits, since it implies disregarding the existence of KAM tori or the chaotic
saddle in the nonhyperbolic and hyperbolic dynamical regimes, respectively. As
is well known, KAM tori prevent particles to escape exponentially for long
periods of time due to the phenomenon of stickiness, by which particles
starting in a chaotic region can stick to the vicinity of the boundary of a
regular region before escaping hillermeier1992 ; lai1992 ; mackay1992 ;
contopoulos2010 ; bunimovich2012 . Therefore, the KAM stickiness clearly
disallows the assumption of equiprobability of the phase space states,
avoiding ergodic motions lebowitz1973 ; zheng1995 ; berdichevsky1991 .

On the other hand, whereas the presence of the chaotic saddle can make the
decay law exponential, it is also responsible for the existence of particles
that evolve and spend a finite amount of time near its vicinity before
escaping. Specifically, a significant number of particles approach the chaotic
saddle following its stable manifold, and move away from its vicinity
following its unstable manifold ott1993 . In this manner, despite the fact
that hyperbolic chaos rules the dynamics in these situations, the particle’s
motion is non-ergodic before escaping because the saddle vicinity occupies
only a subset of the energy shell. Then, once again, we can not expect
equiprobability of the phase space states, when following any escaping
trajectory compatible with some value of the energy. Briefly speaking, in
terms of measure theory, the measure that characterizes the escaping process
from the chaotic saddle along the so-called unstable manifold is very
different from the natural measure of a hyperbolic and closed Hamiltonian
system altmann2008  that is ergodic and, therefore, that can be described by
the microcanonical ensemble throughout the chaotic sea in its phase space
zheng1995 .

One of the purposes of the present work is to extend ergodic decay laws to
relativistic chaotic scattering, to study numerically their scope of
applicability regarding the system’s dynamical regimes, and to propose other
ways to widen these laws taking into account statistical methods that do not
rely on the ergodic hypothesis. Firstly, we provide a complete description of
the Hénon-Heiles model in Sec. II, including striking aspects of escapes, such
as the critical time. In Sec. III we thoroughly discuss the concept of
ergodicity in open Hamiltonian systems, and provide numerical evidence of the
limited applicability of ergodic decay laws at the end of the section.
Furthermore, we derive an analytic ergodic decay law in a relativistic version
of the Hénon-Heiles system in Sec. IV. We also derive this analytical decay
law considering that escapes are measured by a noninertial reference frame
comoving with a particle that never escapes from the scattering region.
Finally, in Sec. V, we conclude this work with a discussion of our findings
and suggesting possible ways to obtain a more accurate characterization of
decay laws utilizing statistical mechanics.

##  II Model description

The Hénon-Heiles system was introduced in 1964 to study the existence of a
third integral of motion in models of galaxies with axial symmetry henon1964 .
We utilize a dimensionless and conservative version of the Hénon-Heiles system
where the particle dynamics is governed by the Hamiltonian

$H(p_{x},p_{y},x,y)=\frac{p_{x}^{2}+p_{y}^{2}}{2}+V(x,y),$ (1)

where $p_{x}$ and $p_{y}$ are the momenta, and $V(x,y)$ is a nonlinear
potential function which depends on $x$ and $y$. In the Hénon-Heiles system,
the potential well is written as

$V(x,y)=\frac{1}{2}(x^{2}+y^{2})+x^{2}y-\frac{1}{3}y^{3}.$ (2)

The value of the mechanical energy is determined by the Hamiltonian, which is
a conserved quantity, $H(p_{x},p_{y},x,y)=E$. When the energy is above a
threshold value called the escape energy, $E_{e}=1/6$, most trajectories are
unbounded so that the system enters an open regime, where escapes are allowed
and the potential well exhibits three symmetric exits, as shown in Fig. 1. In
addition, we define the quantity $\Delta E\equiv E-E_{e}$ for convenience in
computing decay laws in next sections.

The corresponding equations of motion are given by

$\displaystyle\dot{x}=$ $\displaystyle\frac{\partial H}{\partial
p_{x}}=p_{x},$ $\displaystyle\dot{p_{x}}=-\frac{\partial H}{\partial
x}=-x-2xy,$ (3) $\displaystyle\dot{y}=$ $\displaystyle\frac{\partial
H}{\partial p_{y}}=p_{y},$ $\displaystyle\dot{p_{y}}=-\frac{\partial
H}{\partial y}=y^{2}-x^{2}-y.$

The fixed points of the system can be computed from these latter equations. In
this manner, we see that the system has three fixed points located at the
saddle points of the potential well, namely, $(x_{s},y_{s})=(0,1)$,
$(\sqrt{3}/2,-1/2)$ and $(-\sqrt{3}/2,-1/2)$, and another fixed point right in
the minimum of the well, $(x_{m},y_{m})=(0,0)$. In addition, there exist three
highly unstable periodic orbits known as Lyapunov orbits contopoulos1990 . If
a trajectory crosses through one of them, it escapes towards infinity and
never returns back to the scattering region. The Lyapunov orbits are placed
extremely close to the saddle points in the present version of the Hénon-
Heiles system, and therefore it is convenient to define the scattering region
as the allowed region of the potential well delimited by the lines that cross
the saddle points, as can be visualized in Fig. 1(b). Unless otherwise
specified, trajectories have been computed by means of an adaptive Runge-
Kutta-Fehlberg method burden2015  with a relative tolerance of $10^{-12}$ for
each solution of the system of equations and a maximum value of the
integration step of $5\cdot 10^{-3}$, small enough to conserve the mechanical
energy along the trajectory.

Figure 1: (a) Hénon-Heiles potential,
$V(x,y)=\frac{1}{2}(x^{2}+y^{2})+x^{2}y-\frac{1}{3}y^{3}$. (b) Closed and open
isopotential curves projected into the configuration space $(x,y)$. We
arbitrarily establish the exits of the potential well at the saddle points:
Exit 1 (red), Exit 2 (green) and Exit 3 (blue). In this way, the scattering
region is the allowed region of the potential well delimited by the saddle
points lines. The potential minimum is located right at the center of the well
(black dot). (c) Scattering region for $E=0.19$ (light gray) and a particle
trapped in a KAM structure describing a non-escaping orbit typical of the
nonhyperbolic regime. (d) Exit basins and non-escaping orbits in black
associated with $E=0.19$, as computed by means of the tangential shooting
method aguirre2001  with a 500$\times$500 resolution. (e) Similarly, for
$E=0.25$, an escaping orbit starting from the same initial condition that in
(c), but in the hyperbolic regime. (f) Exit basins for $E=0.25$.

The underlying dynamics of open systems is determined by the presence of KAM
tori in phase space. They are sticky sets of quasiperiodic orbits that may
delay escapes or even trap particles forever within the scattering region, as
shown in Fig. 1(c). Then, there exist two open dynamical regimes. On the one
hand, the nonhyperbolic regime where KAM tori coexist with the chaotic saddle
and the phase space exhibits regions where the dynamics is regular surrounded
by a chaotic sea sideris2006 . An example of the typical nonhyperbolic
topology of the exit basins is presented in Fig. 1(d). As already indicated,
the KAM tori stickiness may significantly influence the global properties of
the dynamics, such as the decay law, which can become algebraic at long times.

On the other hand, if the system energy is increased, the KAM tori that
dominate the dynamics disappear and thus the hyperbolic regime begins. The
destruction of such KAM structures has been studied recently, yielding
approximately an energy value of $E\approx 0.2309$ for the transition between
dynamical regimes nieto2020 . After this energy value, the chaotic saddle
dominates the dynamics, making the system chaotic (see, e.g., Figs. 1(e) and
1(f)). Nonetheless, some remnant KAM tori may survive for $E>0.2309$, after
the destruction of the last relevant KAM tori. The presence of such remnant
structures in the hyperbolic regime show that purely hyperbolic dynamics is
unlikely to exist duarte1994 . However, they have been extensively studied and
found to be extremely small and difficult to locate barrio2009 ; barrio2020 ,
so that their possible stickiness barely affects the decay law: an ensemble of
particles escapes conforming an exponential law at long times in the
hyperbolic regime. Finally, this situation is different in area-preserving
maps, where there exist apparently relevant sticky KAM structures for ranges
of parameter values arbitrarily large miguel2015 .

In order to illustrate how the volume in phase space occupied by the KAM tori
evolves as the system energy is varied, we analyze the well-known surface of
section $(y,p_{y})$ located at $x=0$, which is a convenient plane transverse
to the orbits. Similar surfaces of section have been extensively used in area-
preserving maps lai1992  and in continuous-time systems as well aguirre2001 .
Thus, particles are launched from the exits of the scattering region with
random initial shooting angles into its interior. These trajectories map the
entire surface of section except the area occupied by the KAM tori. Using a
partition of the surface of section with high resolution, it is easy to
compute what fraction of points have not been visited by the escaping
trajectories, which is the fraction occupied by the KAM tori (see Appendix at
the end of the manuscript for some method specifications). We show in Fig.
2(a) examples of which regions are mapped (darker gray) and those occupied by
the KAM tori (lighter gray) for three values of the energy, and in Fig. 2(b)
the fraction of KAM tori versus the system energy. The shape of the curve is
non-trivial due to the complexity of the KAM tori destruction barrio2020 .
Although the volume of the phase space occupied by tori KAM is residual for
$E>0.208$ (dashed red line), it has been shown that KAM structures can
influence up to $E\approx 0.2309$ still some global properties of the
dynamics, such as the unpredictability of the final state of the system
nieto2020 . For this reason, along this work, we set the hyperbolic regime in
the energy range $E>0.2309$.

Figure 2: Computing the fraction of KAM tori in the surface of section
$(y,p_{y})$ placed at $x=0$, where only escaping trajectories are launched
from the exits of the scattering region (more details in the text and in
Appendix). (a) For $E=0.19$, a trajectory launched from the Exit 3 with eight
crossings though the surface of section before escaping through Exit 3, again.
We also show three examples of surfaces of section for $E=0.17$, $0.19$ and
$0.25$, where the mapped regions (darker gray) and KAM tori regions (lighter
gray) are depicted. (b) Fraction of the surface of section occupied by KAM
tori $f_{k}$ versus the system energy, specifically, 50 equally spaced values
belonging to the interval $E\in[0.17,0.25]$. The dashed lines are placed at
$E=0.208$ (red) and $0.2309$ (black), respectively. We have set a
500$\times$500 grid resolution for the surface of section and $10^{6}$ initial
conditions have been launched for each energy value.

In addition to KAM tori, a relevant set of the phase space is the chaotic
saddle. Specifically, it is an invariant set with zero-Lebesgue measure
resulting from the intersection of the stable and unstable manifolds. The
stable (unstable) manifold consists of a self-similar fractal set of an
uncountable number of orbits that approach the saddle as $t\to\infty$
($t\to-\infty$) ott1993_2 . On the other hand, the stable manifold can be
understood as the boundary of the exit basins, which exhibits the Wada
property in the Hénon-Heiles system aguirre2001 . For this reason, the stable
manifold of the chaotic saddle (and therefore itself) can be characterized by
means of the fractal dimension of the exit basins boundary. As the energy is
increased, such boundary becomes smoother, the fractal dimension of the
chaotic saddle decreases and the system becomes more predictable nieto2020 .

An important concept, such as the critical time, has not been addressed in
this system yet. The critical time is defined as the escape time spent by a
particle launched from a saddle point of the potential towards the center of
the well bolotin2008 ; bolotin2010 . More specifically, it is the maximum time
that straight-line trajectories can spend inside the scattering region before
escaping. We name these trajectories as the critical trajectories and there
exist as many as saddle points the potential has (see Fig. 3(a)). We compute
this time in the Hénon-Heiles system utilizing the critical trajectory
associated with Exit 1,

$t_{c}(E)=2\int_{y_{s}}^{y_{f}}\frac{dy}{p}=2\int_{1}^{y_{f}}\frac{dy}{\sqrt{2E-y^{2}+\frac{2}{3}y^{3}}},$
(4)

where $y_{s}=1$ is the $y$-coordinate of the saddle point located at the Exit
1, and $y_{f}$ is the $y$-coordinate of the point of the configuration space
which satisfies the equality $E=V(0,y_{f})$. Finally,
$p=\sqrt{p_{x}^{2}+p_{y}^{2}}$ is the modulus of the particle momentum.

Interestingly, critical trajectories are not disturbed by the presence of KAM
tori in light of Fig. 3(b). In addition, already in Fig. 3(a), we also show
the critical trajectories even when the potential well is closed, for
$E=0.16$. In this example, the three trajectories remain bounded along
straight lines over 100 time units, bouncing back and forth against the
potential barriers. They can be simulated over any arbitrarily longer time and
would describe the same straight-line bounded behavior, because the so-called
critical trajectories are closely related to some of the stable normal modes
of the Hénon-Heiles system in its closed regime barrio2009 ; barrio2020 . More
examples of critical trajectories are shown when the system is open in Fig.
3(a) again. Therefore, as the critical trajectories are not deviated by the
likely influence of KAM tori or the chaotic saddle, we can safely state that
critical trajectories are the orbits in the open regime that spend the longest
time without suffering any effect of the sets that govern the dynamics. In
order to demonstrate the latter, we compute the exit basins shown in Fig. 3(c)
by shooting particles with a constant initial angle $\theta_{0}=3\pi/2$. We
notice that there exists a smooth region of initial conditions (depicted in
red) near the Exit 1 critical trajectory, that is surrounded by the typical
fractal regions produced by the influence of the sets that governs the
dynamics, where the final state of particles is subject to chaos. As this
smooth region becomes larger as the system energy increases, it implies that a
significant amount of particles starting from it might escape having only
suffered negligible effects from such mentioned sets on their trajectories and
their escape times.

Figure 3: (a) Critical trajectories (red, green and blue) in the configuration
space for $E=0.16$ (closed regime) and $E=0.17$, $0.19$ and $0.25$ (open
regime). We also simulate another trajectory (yellow) initialized near the
critical trajectory associated with the Exit 1. (b) Critical trajectories in
the space $(y,p_{y})$ when $x=0$. They allow to visualize the presence of KAM
tori (black) and their proximity to the critical trajectories. (c) Exit basins
as computed by shooting particles always with a constant initial angle,
$\theta_{0}=3\pi/2$. (d) Critical times versus the system energy as a
numerical result of Eq. (4) (black line) and as an average from the escape
times of three critical trajectories (red triangles). On the other hand,
$\bar{t}_{e}$ (red squares) represents the average escape times of $5\cdot
10^{5}$ particles launched randomly along the critical trajectory associated
with the Exit 1, and $\bar{t}_{e}$ (yellow) means the average escape times of
particles launched with $\theta_{0}=3\pi/2$ and randomly distributed within
the delimited region by the dashed lines depicted in (a) and (c). Finally, the
vertical dashed line splits the nonhyperbolic and hyperbolic regimes at an
approximate energy value $E=0.2309$.

The critical times versus the system energy as a numerical result of Eq. (4)
and as an average computed from the escape times of three critical
trajectories are in perfect agreement as shown in Fig. 3(d). As pointed out
above, the values of critical times are not disturbed by the presence of KAM
tori. Furthermore, we compute the average escape times of two different
ensembles of particles. On the one hand, we randomly launch particles along
the Exit 1 critical trajectory, i.e., with $x_{0}=0$ and $\theta_{0}=3\pi/2$
as initial conditions. As expected, the average escape times computed are
always less that the critical times, since simulated trajectories are just
shorter parts of the Exit 1 critical trajectory. On the other hand, we
simulate particles with $\theta_{0}=3\pi/2$ again, but randomly distributed
within the region delimited by the limits of the Exit 1 for $E=0.17$, i.e.,
$x_{0}=\mp\sqrt{2\Delta E/3}$, where $\Delta E=0.17-1/6$ (see the dashed black
lines in Fig. 3(a) or 3(c)). Simulations reveals that average escape times are
less than critical times for approximately $E>0.18$ because the majority of
these particles bounce against the potential barrier and escape through the
Exit 1, having evolved shorter trajectories that critical trajectories in
general. We recall that all these particles belong to the smooth region
mentioned above.

Finally, according to the results about critical times, we state that all
escaping orbits after the critical time are affected by KAM tori and the
saddle in the nonhyperbolic and hyperbolic regimes, respectively. Conversely,
as indicated above, there exist some escaping orbits that are negligibly
affected by these sets before the critical time. Importantly, this fact is
reflected in the decay curves: a sudden change of the tendency of the escapes
occurs when the critical time is achieved. For instance, escapes take place
more slowly after the critical time when the dynamics is nonhyperbolic because
all escaping orbits are affected by KAM tori stickiness, as shown in the next
section.

##  III Decay law in the Newtonian Hénon-Heiles system

Theoretically, a Hamiltonian system is ergodic when its energy surface
uniquely consists of a chaotic sea, due to the existence of a dense set of
highly unstable periodic orbits densely embedded in it ott1993_2 . However, a
problem may arise when testing if the motion is ergodic or not in open
systems, since it is necessary to simulate trajectories trapped in the
scattering region for a large amount of time. For example, in open systems,
simulated trajectories generally escape before accomplishing long enough times
to fill densely the energy surface and test the assumption of ergodicity. In
this regard, a modified Hénon-Heiles system with reflecting walls, where
elastic collisions occur, was proposed to prevent particles from escaping
zheng1995 . The authors provide numerical results claiming that this modified
Hénon-Heiles system behaves ergodically in the hyperbolic regime when the
walls are placed far from the origin of the potential well. Thus, one might
think that if the reflecting walls are placed infinitely far from the origin,
the Hénon-Heiles system would be recovered and, following the computational
results shown above, one would obtain that the open and hyperbolic Hénon-
Heiles system behaves ergodically. However, this limit case is unattainable in
computational practice, since the Hénon-Heiles system is open and does not
allow any escaped trajectory to return to the scattering region again. As we
will discuss below, the fact that trajectories escape forever causes the
escaping process to be non-ergodic.

Since the chaotic saddle is a zero-Lebesgue-measure set and the definition of
ergodicity only requires that almost all trajectories in the energy surface
fill it densely moore2015 , the presence of such an invariant set should not
hinder a priori the application of the ergodic hypothesis. However, saddle
sets separate the phase space into disjoint sets, what invalidates the ergodic
hypothesis. This is evinced by the fact that a trajectory launched directly
towards any of the exits can leave the scattering region without returning
(for instance, see any critical trajectory). Moreover, the phase space
fragmentation is manifested by the existence of large and differentiated exit
basins, as can be appreciated in Fig. 1(f). Therefore, the approximation of
ergodicity to derive decay laws can only be made for energies close to the
energy escape, i.e., when the exits are vanishingly small and the exit basins
exhibit a completely fractalized structure. Furthermore, and as we show
numerically right ahead, this hypothesis also demands to consider short enough
times. Finally, as is well known, we recall that the stickiness and the
possible entrapment of the particle developed by KAM tori prevent that a set
of trajectories with a Lebesgue measure different from zero satisfies the
ergodic hypothesis lebowitz1973 ; zheng1995 ; berdichevsky1991 .

Having stated the problems of assuming the ergodic hypothesis in open systems,
we elucidate below the range of applicability of an ergodic decay law to
describe exponential escapes. At first, the evolution of the number of
particles inside the scattering region for any hyperbolic system can be
expressed without loss of generality as

$N(t)=N_{0}e^{-\alpha(E)t},$ (5)

where $\alpha(E)$ is the decay rate and depends on the mechanical energy of
the system. A decay law, $\alpha_{e}$, has been derived analytically in the
case of the open Hénon-Heiles system zhao2007 , obtaining

$\alpha_{e}(\Delta E)=\frac{\sqrt{3}\Delta E}{S(\Delta E)},$ (6)

where the function $S(\Delta E)$ is the area of the scattering region as
defined in Sec. II, which can be computed by utilizing a Monte Carlo method.
The calculus of the ergodic decay rate, $\alpha_{e}$, is accomplished by means
of the microcanonical ensemble, which is the only ensemble able to describe
the statistical properties of ergodic systems gutzwiller1990 .

To test the ergodic decay law, we launch a large amount of particles
($N_{0}=5\cdot 10^{5}$) from different initial conditions randomly distributed
in the configuration space of the scattering region, i.e., with random values
for $x_{0}$, $y_{0}$ and $\theta_{0}$. For convenience, the evolution of the
number of particles inside the potential well is shown in natural logarithmic
scale, i.e., $\ln N(t)$ in Fig. 4. This representation enables the reader to
infer easily an exponential decay when a linear relation between $\ln N$ and
$t$ appears. We take into account several time intervals to fit the
exponential behaviors.

Figure 4: We consider short times ($t_{5\%}$, $t_{10\%}$ and $t_{15\%}$), the
critical time ($t_{c}$) and long times ($t_{2c}\equiv 2t_{c}$, $t_{99\%}$ and
$t_{99.9\%}$). (a) Decay curves (gray) when the underlying dynamics is
nonhyperbolic for $E=0.17$, $0.19$, $0.21$ (solid lines) and $0.23$ (dashed
line). (b) Exponential decay curves (gray) when the underlying dynamics is
hyperbolic, for $E=0.235$, $0.25$, $0.275$ (solid lines) and $0.29$ (dashed
line). (Inset plots) The first instants of the decay curves. (c) The evolution
of the fitting times versus the system energy. (d) Again the evolution of the
fitting times versus the system energy, but the evolution of the long times is
included.

In order to ascertain how the time intervals affects the decay, firstly, we
define the short-time quantities $t_{5\%}$, $t_{10\%}$ and $t_{15\%}$, where,
for instance, $t_{10\%}$ means the time at which the $10\%$ of the particles
starting inside the scattering region have escaped. These short times are
usually smaller than the critical time and thus are located at the beginning
of the time evolution. Secondly, we consider the critical time, $t_{c}$. As
indicated in the previous section, notice that the decay curve suffers a
slight change of tendency when this time is reached in both dynamical regimes,
although this change is even more subtle in the hyperbolic regime. Finally, we
define the long-time quantity $t_{2c}$, which is twice the critical time for
nonhyperbolic dynamics and, on the other hand, we also establish the long-time
quantities $t_{2c}$, $t_{99\%}$ and $t_{99.9\%}$ for hyperbolic dynamics. Note
that, in the nonhyperbolic regime, the exponential decay only appears during
the first instants of the time evolution, provided that KAM tori turn the
decay law algebraic at long times, as is clearly shown in Fig. 4(a). For this
reason, we only compute exponential fittings until the time $t_{2c}$, when KAM
tori are present. However, the presence of the chaotic saddle makes the
escapes exponential in the hyperbolic regime (see Fig. 4(b)), and thus we are
able to fit the exponential decay until longer times, such as $t_{99\%}$ or
even $t_{99.9\%}$. Interestingly, we observe that the fraction of trapped
particles for $E=0.23$ is greater than for $E=0.21$. This result is in
agreement with the fact that the phase space volume occupied by KAM tori
usually follows a non-trivial tendency as the energy changes. Furthermore, we
show in Figs. 4(c) and 4(d) how the time quantities defined above evolve as
the system energy increases.

The fitted decay rates $\alpha_{a,b}$ associated with the time intervals
$t\in[t_{a},t_{b}]$ are displayed in Fig. 5. For example, the exponential rate
$\alpha_{0,10\%}$ has been computed considering only a very small part of the
decay curve, $t\in[0,t_{10\%}]$. We observe that the shorter the times
considered in exponential fitting, the better the agreement is between the
ergodic decay rate, $\alpha_{e}$, and the fitted rate, such as
$\alpha_{0,5\%}$ (see Fig. 5(a)). This is explained because the ergodic decay
law does not take into account the KAM tori and the chaotic saddle, whose
traces in particle paths are more noticeable when the longer the time
particles have evolved inside the scattering region. Thus, as a consequence,
the ergodic decay law overestimates the decay rate $\alpha_{c,2c}$ when KAM
tori are present and escapes are delayed by them. However, in the hyperbolic
regime, the ergodic decay law underestimates all long-time decay rates, such
as $\alpha_{c,2c}$, $\alpha_{2c,99\%}$ and $\alpha_{99\%,99.9\%}$, because the
unstable manifold of the chaotic saddle affects the escapes at long times,
accelerating them (see Fig. 5(b)). In addition, we note that as the system
energy is increased, the chaotic saddle also enhances its effects on escapes
at short times. Hence, the ergodic decay rate underestimates escapes when the
energy system is high even at short times. This last argument can be applied
when analyzing why the fitted rate $\alpha_{0,c}$ does not agree with the
ergodic decay law neither, because the influence of the sets that governs the
dynamics is strong enough to prevent particles from being distributed
according to the microcanonical ensemble at the critical time as well.

Furthermore, we would like to highlight that the decay rates for long times
may be useful to quantify whether particles always escape in the same way over
the three different periods of the time evolution that we establish in the
hyperbolic regime, namely $\alpha_{c,2c}$, $\alpha_{2c,99\%}$ and
$\alpha_{99\%,99.9\%}$. For example, these decay rates can be affected by the
existence of extremely small KAM tori that may survive even after the
destruction of the dominant KAM tori barrio2009 ; barrio2020 . In our
simulations, we have not found trapped particles for energies greater than
$E>0.2309$, neither long-time algebraic decays for at least the $99.9\%$ of
the initialized particles in the scattering region. Therefore, such remnant
KAM structures are not decisive in the exponential escape regime. Nonetheless,
we observe in Fig. 5(b) that the decay rates of the two last periods
considered, $\alpha_{2c,99\%}$ and $\alpha_{99\%,99.9\%}$, exhibit a
fluctuating evolution when plotted against the system energy, whose cause
deserves further research. So far, we have found that hyperbolic decay rates
measured at any given sufficiently long time are underestimated by the ergodic
decay law, which implies a clear limitation of this kind of decay laws in
describing escapes in open systems.

Figure 5: Exponential fittings of the decay curves at (a) short times, (b)
critical and long times. We recall that the dashed line again splits the
nonhyperbolic and hyperbolic dynamics when $E=0.2309$. For the sake of
clarity, although the decay law $\alpha_{e}$ has been computed in terms of
$\Delta E$, we show the results of fittings in terms of $E$ for simplicity.

Now, we present the _decay basins_ for hyperbolic dynamics in Figs. 6(a–c) to
visualize which particles escape during each time interval and how they are
distributed inside the scattering region. Furthermore, they relate the escape
times and the exit basins built by utilizing the tangential shooting method
aguirre2001 . We observe that an important fraction of particles escapes
before the critical time. These particles start their escaping trajectories
far from the fractal boundaries of the exit basins, i.e., far from the
influence of the chaotic saddle. This confirms that the critical time is the
maximum escape time for a trajectory that is not affected by the sets that
governs the dynamics. Conversely, all the initial conditions associated with
escape times higher than the critical time are located close to the fractal
boundaries of the exit basins.

Figure 6: (a–c) Decay basins for $E=0.235$, $0.25$ and $0.275$, respectively.
The color code indicates the initial conditions from which escapes occur
before the short time $t_{15\%}$ (lighter red), the critical time $t_{c}$
(red), twice the critical time $t_{2c}$ (darker red), and finally those from
which particles have not already escaped by $t_{2c}$ (black). (Inset plot)
Exit basins associated with the shown decay basins. (d–f) Average escape
times, again, for $E=0.235$, $0.25$ and $0.275$. We launch $2.5\cdot 10^{3}$
particles for each initial condition $(x_{0},y_{0})$ with random shooting
angles belonging to the interval $\theta_{0}\in[0,2\pi)$. As explained in the
text, the average escape times can be interpreted as the regions where it is
more likely to find a particle in the configuration space. (Inset plot) Zoom-
in near the potential barrier between Exit 1 and Exit 3, where a low-
probability region is found in the three cases.

Nonetheless, we are interested in the escape time distribution by shooting
particles from initial conditions randomly. Hence, we show here the average
escape times in the configuration space in Figs. 6(d–f), where for each pixel
we launch a large amount of particles and compute the mean of their escape
times. In the three cases shown, it can be observed an apparently fractal
structure that is reminiscent of the typical structures that appear in
scattering functions when a Smale horseshoe-like set rules the dynamics, such
as the chaotic saddle. In fact, this structure is formed by the initial
conditions that spend long escape times, i.e., the initial conditions close to
the stable manifold of the chaotic saddle ott1993_2 . Therefore, the images of
average escape times provide a clear manifestation of the presence of the
chaotic saddle when particles are launched randomly.

Moreover, we study how the chaotic saddle explicitly affects particle
trajectories. Then, again, by launching particles with random values for
$x_{0}$, $y_{0}$ and $\theta_{0}$, and taking into account every position of
their orbits, we obtain the regions where a hypothetical particle is more
likely to be found before its eventual escape. Not surprisingly, the images
obtained are the same as those provided by calculating the average escape
times (see Figs. 6(d–f)), since trajectories starting close to the chaotic
saddle spend longer times within the potential well: they come close to the
chaotic saddle along its stable manifold, stay in its neighborhood, and leave
it along its unstable manifold. In light of the results, the particle is more
likely to be found near the chaotic saddle in its escaping process. In
addition, we also detect a region of low probability of finding the particle,
as can be visualized in the inset plots of Figs. 6(d–f).

Finally, in Fig. 7 we show the instantaneous distribution of particles at some
relevant times of the decay law. The regions that are emptied in the first
place are those close to the exits from which particles can escape if the
shooting angle points directly towards the exit. Subsequently, when the
critical time is reached, a symmetric structure can be inferred in the
particle distribution, which is closely related to the manifestation of the
chaotic saddle (see Figs. 6(d–f)). Interestingly, the fact that trajectories
spend a finite period of time in the chaotic saddle vicinity implies that
there are zones of the scattering region full of particles and others empty
during the time evolution. For instance, as indicated above, there is an
almost zero probability of staying in the space nearby the wall between two
exits, because the chaotic saddle does not occupy that zone (see the inset
plots of Figs. 7(a–b)). Hence, the presence of the chaotic saddle is
incompatible with an equiprobability of being in any dynamical state
compatible with energy, as assumed by the microcanonical ensemble.

Figure 7: Particle distribution at (a) $t_{15\%}$, (b) $t_{c}$ and (c)
$t_{99\%}$, respectively, by means of a histogram, for $E=0.25$. The histogram
color code consists of a fading between yellow (highly occupied area) and dark
blue (lowly occupied area). The white color means an area depleted of
particles. (Inset plots) Zoom-in of the particle distribution near the
potential barrier between Exit 1 and Exit 3. We depict the present particles
(black dots) in the allowed scattering region compatible with energy (gray).

##  IV Decay laws in the relativistic Hénon-Heiles system

We now employ a relativistic version of the Hénon-Heiles system to extent our
claims to the relativistic realm. Hence, we consider a particle interacting in
the limit of weak external fields whose dynamics is governed by the
conservative and dimensionless Hamiltonian lan2011 ; chanda2018 ; kovacs2011 ;
calura1997

$H(p_{x},p_{y},x,y)=c\sqrt{c^{2}+p_{x}^{2}+p_{y}^{2}}+V(x,y),$ (7)

where $c$ is the value of the speed of light and $V(x,y)$ is again the Hénon-
Heiles potential described above. Although the equations of motion are
different from the Newtonian ones, the four fixed points of the relativistic
system still remain at the three saddle points of the potential well and at
its minimum. For this reason, we have defined the relativistic scattering
region as in the Newtonian Hénon-Heiles system. Furthermore, the Lorentz
factor is defined as

$\gamma=\frac{1}{\sqrt{1-\frac{v^{2}}{c^{2}}}}=\frac{1}{\sqrt{1-\beta^{2}}},$
(8)

where $v$ is the speed of the particle and $\beta=v/c$ the ratio between the
speed of the particle and the speed of light. The Lorentz factor $\gamma$ and
the quotient $\beta$ are two ways to express how large is $v$ compared to $c$,
and they vary in the ranges $\gamma\in[1,+\infty)$ and $\beta\in[0,1)$,
respectively. We shall use $\beta$ as the relativistic parameter for
convenience in showing our numerical results.

We utilize a method introduced in bernal2017 ; bernal2018  to increase the
kinetic energy of the system to the relativistic regime (the value of
$\beta$). Firstly, the maximum speed that the particle can have along its
trajectory before escaping can be defined as $v_{m}\equiv\sqrt{2E}$, bearing
in mind that the maximum kinetic energy is found at the potential minimum,
$V(x_{m},y_{m})=0$, and $E$ is the total mechanical energy of the Newtonian
version of the Hénon-Heiles. Thus, the method is based on the fact that
$\beta$ is a quantity that equivalently depends on $v_{m}$ and $c$, hence
working on a dimensionless system it is possible to increase the value of
$\beta$ by fixing the value of $v_{m}$ and decreasing the numerical value of
the speed of light. In this manner, if the value of $\beta$ is almost null,
$v_{m}\ll c$, the speeds of the particle during its evolution within the
scattering region only represents a very low percentage of the speed of light.
In this case, we recover the Newtonian version of the Hénon-Heiles system. On
the contrary, if the value of $\beta$ is close to one, the speed of the
particle represents a high percentage of the speed of light. For the sake of
clarity, although the numerical value of $c$ is modified to select a value of
$\beta$, this numerical value of $c$ always represents the speed of light in
the system because its dynamics is governed by the Hamiltonian (7). For more
method’s details see fernandez2020 .

We have arbitrarily established $v_{m}\approx 0.679$ along this section, which
corresponds to the Newtonian energy $E\approx 0.2309$, to focus on hyperbolic
dynamics in this relativistic Hénon-Heiles system. Summarizing, the Newtonian
and hyperbolic system becomes relativistic and increasingly hyperbolic as
$\beta$ increases. In Fig. 8, we show the effect of $\beta$ on the dynamics
and the basin topology: the particle generally escapes more quickly and the
fractal boundaries of exit basins are smoother as $\beta$ increases. Finally,
the computations along this section have been carried out only until
$\beta=0.75$, since the relativistic Hamiltonian (7) has been derived in the
limit of weak external fields.

Figure 8: Scattering regions (pale gray) and escaping trajectories starting
from the same condition but for (a) $\beta=0.25$, (b) $0.5$ and (c) $0.75$,
where the underlying dynamics is hyperbolic. As $\beta=v_{m}/c$ always holds,
the computational value of the speed of light changes inevitably when the
value of $\beta$ changes. We have arbitrarily set $v_{m}\approx 0.679$, and
therefore the values of the speed of light in the computations are (a)
$c\approx 2.718$, (b) $1.359$ and (c) $0.906$, respectively. (Inset plots)
Exit basins associated.

We also concentrate our efforts on exponential decay laws along this section
and derive a decay rate based on ergodic motion inside the scattering region.
Firstly, instead of the momenta $p_{x}$ and $p_{y}$, Eq. (7) can be rewritten
in terms of the total momentum of the particle $p$ and the angle $\theta$
formed by the speed vector and the positive $x$-axis, obtaining thus

$H(p,x,y,\theta)=c\sqrt{c^{2}+p^{2}}+V(x,y).$ (9)

Importantly, we recall that the total momentum of the particle is related to
its speed in this dimensionless version of the system by $p\equiv\gamma v$.
Furthermore, we express the total relativistic energy as $\Delta E\equiv
E-E_{e}=\gamma c^{2}+V-1/6$ for convenience when computing the decay rate. In
addition, there exists a relation between $\Delta E$ and $\beta$, which is

$\Delta E(\beta)=\frac{c^{2}}{\sqrt{1-\beta^{2}}}-\frac{1}{6},$ (10)

since the potential energy is null, $V(x_{m},y_{m})=0$, when we define $\beta$
as the relativistic parameter of the system. The motion in the well is assumed
to be ergodic, and therefore the phase space distribution inside the
scattering region for a given energy $\Delta E$ is expressed by

$\rho(p,x,y,\theta)=\frac{\delta\left(E-H(p,x,y,\theta)\right)}{\int
dpdxdyd\theta\delta\left(E-H(p,x,y,\theta)\right)},$ (11)

where $\delta(E)$ represents the Dirac delta function. Here, we are interested
in the distribution of the variables in the configuration space. In this
manner, we integrate Eq. (11) along the interval $p\in[0,\infty)$, yielding

$\rho(x,y,\theta)=\frac{\gamma}{2\pi\int_{S(\Delta E)}dxdy\gamma},$ (12)

where the integral is bounded inside the scattering region, and the Lorentz
factor $\gamma$ depends on the spatial coordinates because the system is
conservative and every point of the well is associated with a specific kinetic
energy value. Notice that if the limit $c\to\infty$ (or $\gamma\to 1$,
equivalently) is considered, the Newtonian probability density
$\rho(x,y,\theta)=1/2\pi S(\Delta E)$ is recovered, as derived in zhao2007 .
For simplicity, we define $\Gamma(\Delta E)\equiv\int_{S(\Delta
E)}dxdy\gamma$, where the integral can be solved numerically by means of a
Monte Carlo method.

In resemblance to the Newtonian case, the number of particles remaining inside
the scattering region under hyperbolic dynamics at the time $t$ is
$N(t)=N_{0}e^{-\alpha_{e}t}$, where the $\alpha_{e}(\Delta E)$ is again the
decay rate, but relativistic. It can be computed by setting that the flux of
escaping particles through an exit of the well, for instance, the Exit 1, is

$\iint
dxd\theta\rho\textbf{v}\cdot\hat{\textbf{n}}=\int_{x_{a}}^{x_{b}}dx\int_{0}^{\pi}d\theta\rho(x,y,\theta)v(x,y)\sin\theta,$
(13)

where the quantity $\rho\textbf{v}$ means a current density vector and
$\hat{\textbf{n}}$ is the normal vector to the Exit 1, which points outwards
along the direction of the $y$-axis zhao2007 . We choose the Exit 1 for
convenience, since the integral then remains defined along the segment
delimited by $x\in[x_{a},x_{b}]$ and $y=1$, as shown in Fig. 9(a).
Specifically, the points

$x_{a,b}=\mp\sqrt{\frac{2}{3}(\Delta E-c^{2})}$ (14)

satisfy the condition $E=c^{2}+V(x_{a,b},1)$, where the particle kinetic
energy vanishes. Finally, bearing in mind the triangular symmetry of the well,
the flux of escaping particles is the same through all exits, and therefore
the decay rate of the system is three times the decay rate through the Exit 1,
yielding the integral

$\alpha_{e}(\Delta E)=\frac{3}{\pi c\Gamma(\Delta
E)}\int_{-\sqrt{\frac{2}{3}(\Delta E-c^{2})}}^{\sqrt{\frac{2}{3}(\Delta
E-c^{2})}}dx\sqrt{\left(\Delta E-\frac{3}{2}x^{2}\right)^{2}-c^{4}}.$ (15)

It is possible to solve this integral by means of elliptic functions,
obtaining the result

$\alpha_{e}(\Delta E)=\frac{4\sqrt{2}\Delta E\sqrt{\Delta
E+c^{2}}}{\sqrt{3}\pi c\Gamma(\Delta E)}\left\\{{\rm E}\left(\frac{\Delta
E-c^{2}}{\Delta E+c^{2}}\right)-\frac{c^{2}}{\Delta E}{\rm
K}\left(\frac{\Delta E-c^{2}}{\Delta E+c^{2}}\right)\right\\},$ (16)

where E$(m)$ and K$(m)$ are the elliptic integrals of first and second kind,
respectively, and $m$ means the squared elliptic modulus. In our case,
$m=(\Delta E-c^{2})/(\Delta E+c^{2})$.

The results are similar to the Newtonian case already discussed. Firstly, as
the value of $\beta$ is varied, the slope of the exponential decay curves
becomes greater, as shown in Fig. 9(b). On the other hand, the shorter the
fitting time of the decay curve is considered, the better the agreement
between the numerical fittings and the ergodic decay rate $\alpha_{e}$, as
displayed in Fig. 9(c). This is due to the fact that the chaotic saddle barely
affects the trajectories escaping during the very first instants of the time
evolution. Hence, as in the Newtonian case, the ergodic decay rate,
$\alpha_{e}$, only takes into account the initial distribution of particles,
in our case, randomly distributed inside the well. However, this initial
distribution of particles is modified by the saddle since the time evolution
begins. That is why only short-time fittings work. This last reasoning is also
valid to interpret why the critical-time fittings do not strictly agree with
the analytical decay rate, as seen in Fig. 9(d). Finally, when the
relativistic system becomes very energetic, the ergodic decay rate
underestimates all long-time fittings, $\alpha_{c,2c}$, $\alpha_{2c,99\%}$ and
$\alpha_{99\%,99.9\%}$, similarly to the Newtonian case.

Figure 9: (a) Scheme of the Exit 1 for $\beta=0.25$. (b) Hyperbolic particle
decay for $\beta=0.25$, $\beta=0.5$ (solid lines) and $\beta=0.75$ (dashed
line). (c) Exponential fittings of the decay curves for only short times in
the inertial relativistic approximation for $\beta\in[0.05,0.75]$. (d) Similar
fittings to (c) but for critical and long times. We clarify that although the
decay law $\alpha_{e}$ has been computed in terms of $\Delta E$, we show the
results of fittings in terms of $\beta$ for simplicity (see Eq. (10)).

Moreover, we would like to mention that the formalism of ergodic decay laws
discussed and developed along this work also fails when the regime is highly
energetic (Newtonian and relativistic schemes included). In Fig. 9(b), for
$\beta=0.75$, it is noticeable that the first change of tendency in the decay
curve happens before the indicated critical time (red triangle). Then, a
problem arises, because the critical time is providing us misleading or
unclear information regarding tendency changes. We argue that in very
energetic cases the scattering region that we define, delimited by the
potential saddle points, becomes far from the actual scattering region,
delimited by the Lyapunov orbits, and does not reflect the actual
phenomenology of escapes in the system. This produces the effect that the
decay curve shows changes in the tendency before (and after as well) the
critical time. In other words, when the system is highly energetic, the
critical time becomes irrelevant to study the decay laws.

Time is absolute in nonrelativistic systems, but the measure of time in
relativistic systems may depend on the choice of the reference frame. Escape
times have been measured by an inertial clock at rest in the potential well in
this work so far. However, it is possible to consider another noninertial
reference frame comoving with a particle that describes a trajectory for a
long time enough to measure the escape times of the other particles. The
hypothetical clock attached to this particle will suffer the well-known time
dilation phenomenon and measure the so-called proper time barton1999 , which
in few words corresponds to the travel time of the twin who leaves the Earth
in the twin paradox. Hence, given an infinitesimal time interval $dt$ as
measured by an inertial clock, the particle clock will measure the time

$d\tau=\frac{dt}{\gamma(t)},$ (17)

where $\gamma(t)$ is the Lorentz factor at time $t$. Both time intervals
satisfy that $d\tau\leq dt$, since the Lorentz factor is usually greater than
unity along the particle trajectory. Importantly, note that it is necessary to
know the complete Lorentz factor evolution of the particle clock to measure
the proper time of an escaping trajectory, which is the result of integrating
Eq. (17).

Nonetheless, as the particle dynamics is bounded under the same energetic
conditions, given a value of $\beta$, the Lorentz factor of all the
trajectories is similar on average at any time $t$ before escaping. Then,
there exists an approximate average value of the Lorentz factor and it can be
reasonably estimated as the arithmetic mean between the maximum and the
minimum values of the Lorentz factor inside the potential well, i.e.,

$\bar{\gamma}(\beta)=\frac{1+\sqrt{1-\beta^{2}}}{2\sqrt{1-\beta^{2}}},$ (18)

as introduced in fernandez2020 . Taking this into consideration, Eq. (17) can
be approximated as $d\bar{\tau}\equiv dt/\bar{\gamma}$, where $\bar{\gamma}$
is a constant. Hence, escape times as measured by an inertial clock will be
higher than escape times of the same trajectories as measured by the
noninertial particle clock. This latter affects quantitatively the exponential
decay laws. The decay rate in the inertial case, for instance, $\alpha_{e}$,
will be smaller than the decay rate in the noninertial case, yielding

$\bar{\alpha}^{\tau}_{e}\equiv\bar{\gamma}\alpha_{e},$ (19)

where the approximation of bounded dynamics have been taken into account.

Observing particle escapes from a noninertial clock makes them occur faster
due to the time dilation both at short and long times, according to the
results shown in Fig. 10. Moreover, the fitted decay rate
$\alpha^{\tau}_{0,5\%}$ for very short times agrees with the noninertial
ergodic decay rate $\bar{\alpha}^{\tau}_{e}$. Thus, the approximation of
assuming an average Lorentz factor inside the scattering region to
characterize escapes works, at least for short times, which corresponds to the
real range of applicability of ergodic decay laws, as proven in previous
sections of this work.

Figure 10: Exponential fittings for short times and long times to compare the
decay rates as measured by an inertial clock at rest in the potential well in
blue ($\alpha_{e}$, $\alpha_{0,5\%}$ and $\alpha_{c,99.9\%}$) and as observed
from a noninertial particle clock evolving chaotically inside the scattering
region in red ($\bar{\alpha}^{\tau}_{e}$, $\alpha^{\tau}_{0,5\%}$ and
$\alpha^{\tau}_{c,99.9\%}$).

##  V Conclusions

The indiscriminate use of the concepts of hyperbolicity, chaos and ergodicity
is sometimes found in the bibliography referring to dynamical systems with a
few degrees of freedom, because chaos is usually responsible for the
fulfillment of ergodicity in this kind of systems. In this manner, we have
contextualized in the present work each of these terms in the open regime of
the Hénon-Heiles system and provided solid evidence that, far from being
equivalent, chaos and the fact that the system is open can also be the precise
reason why its motion is non-ergodic. In particular, we refer to the existence
of a chaotic saddle, which provokes that particles explore all the regions of
the energy surface unevenly before escaping.

Therefore, we have studied numerically an exponential decay law based on
ergodic motions in the Newtonian Hénon-Heiles system. We conclude that this
ergodic decay law is able to quantify adequately escapes only for the very
first instants of the time evolution when particles are launched randomly
distributed. Thus, when the time evolution begins, particles remain no longer
equally distributed throughout the scattering region, and in this way no
longer conform to the microcanonical ensemble, on which the ergodic decay law
is based. Indeed, particles evolve according to other more complex probability
distributions given by the Liouville equation which dictates the evolution of
the probability density according to the sets that govern the dynamics, i.e.,
KAM tori or the chaotic saddle in the nonhyperbolic and hyperbolic regimes,
respectively.

In addition to the detailed study of the Newtonian case, we have extended the
formalism of the ergodic decay laws to a relativistic version of the Hénon-
Heiles system for the first time. The result, obtained by means of measuring
escapes from two inertial and noninertial reference frames, is that we have
encountered the limitations present in the Newtonian case, for the same
reasons. Nonetheless, the current work constitutes a step towards the solution
in the characterization of decay laws in open systems. Bearing in mind the
limitations of the microcanonical distribution, other distributions based on
the sets that governs the dynamics have to be considered in the calculus of
new decay laws. For instance, when KAM tori are present, the particle
distribution must be affected by the phenomenon of stickiness. Therefore, in
this case, there exist some regions near KAM tori where particles explore for
longer times that others. On the other hand, when the chaotic saddle governs
the dynamics, a fractalized particle distribution over the saddle can
definitely describe all the escaping orbits in the hyperbolic regime. In
general, we insist, the measure in the open regime evolves according to
Liouville equation and converges to a fractal probability density defined on
the chaotic saddle.

The critical time has to be considered as well, since it is part of the
phenomenology caused by the invariant sets. Importantly, our numerical study
has allowed us to interpret that the so-called critical time in open
Hamiltonian systems means the maximum escape time that a trajectory is able to
spend inside the scattering region without being affected at all by KAM tori
or the chaotic saddle. Then, the existence of escaping orbits unaffected by
them is reflected in the decay curves, which suffer a sudden change in their
tendency when the critical time is reached. This fact opens some room to
explore particle distributions that evolve in time.

To conclude, we mention some processes in nature with only a few relevant
degrees of freedom where our numerical work can be useful. For example, open
systems can be interpreted as scattering models of nuclear or chemical
reactions, where each exit represents a final rearrangement of atoms and
molecules zotos2015 . Interestingly, these kinds of reactions may occur
chaotically and, in light of the present results, non-ergodically. Another
relevant scattering phenomenon is the chaotic interaction that takes place
among the resulting components of molecules above their dissociation energy.
These processes constitute a sort of laboratory for studying the so-called
quantum-classical correspondences (QCC) present in coupled nonlinear
oscillators lin2013 , and perhaps our conclusions about ergodicity might be
helpful to understand the formation of different bounded and unbounded states.

On the other hand, relaxation of nonequilibrium systems such as glass-forming
melts or soft matter can be modeled by means of a modified Hénon-Heiles
system, whose saddle points are located at different values of the potential
energy toledomarin2018 . More specifically, processes of relaxation consist of
local reorientations of the masses that form the glass structure. In some
cases, these processes are governed by a chaotic saddle, and hence new studies
can be addressed to analyze how non-ergodic dynamics may also affect the
structure and properties of glasses. Finally, the interest in diffusion
phenomena in the so-called soft Lorentz gases has grown because these systems
are able to simulate the electronic transport in graphene-like structures
klages2019 . The unit cell of the triangular lattice of potentials in such
systems is reminiscent to the morphology of the Hénon-Heiles potential.
Numerical works have evinced that soft Lorentz gases exhibit normal diffusion
or superdiffusion due to a non-trivial interaction between trapped and
ballistic periodic orbits that form the KAM tori. Superdiffusion is related to
a dependence on initial conditions, and hence ergodicity is broken in such a
case. Thus, the present work can provide new approaches to elucidate this
phenomenology.

## ACKNOWLEDGMENTS

This work has been supported by the Spanish State Research Agency (AEI) and
the European Regional Development Fund (ERDF, EU) under Projects No.
FIS2016-76883-P and No. PID2019-105554GB-I00.

*

##  Appendix A Comments on computing the fraction of KAM tori in a surface of
section.

As indicated, an ensemble of particles are launched from the exits of the
scattering region with random initial shooting angles into its interior. Such
initial conditions, blue, red and green colored in Fig. 2(a), are all the
possible final conditions of trajectories just before escaping forever, but
with their velocity vector rotated by $\pi$ radians. For clarity, the limit
case of all possible final conditions is reached when the number of initial
conditions is infinite. Thus, due to the time symmetry of conservative
systems, we simulate all escaping trajectories but backwards in time.

Another basic method to compute the fraction of KAM tori consists of counting
how many trajectories do not escape of the total number of trajectories
initially launched nieto2020 . Such a method is usually computationally
expensive as it has to evolve trajectories belonging to KAM tori, which do not
escape by definition, in addition to trajectories with long transients due to
KAM stickiness. The latter trajectories, if the final computation time is less
than their escape times, can be interpreted as belonging to KAM tori, which
would yield an erroneous result of the fraction of KAM tori. Nonetheless, the
alternative method proposed here depends on the number of particles launched:
the greater the number of escaping particles, the greater the probability of
mapping the whole area of the surface of section occupied by the chaotic
regions. Specifically, the regions of the surface of section that are most
complicated to map are those away from the fractal boundaries of exit basins
and tori KAM. These regions, which grow as the energy of the system increases,
are formed by a multitude of trajectories that we must launch, which escape
without exhibiting long transients and map only a few points on the surface of
section.

For completeness, we add two more method’s aspects. On the one hand, note that
it is totally equivalent launching an infinite ensemble of particles by the
proposed method that simulating a single particle trajectory that elastically
bounces off the exits forever zheng1995 , since in both cases one would obtain
the same fraction of KAM tori. On the other hand, this computationally
affordable method can be applied to any system (including area-preserving maps
and continuous-time systems) that exhibits exits well-defined where particles
can be started from.

## References

  * (1) J. M. Seoane and M. A. F. Sanjuán, Rep. Prog. Phys. 76, 016001 (2013).
  * (2) J. Aguirre, R. L. Viana, and M. A. F. Sanjuán, Rev. Mod. Phys. 81, 333 (2009).
  * (3) Y.-C. Lai and T. Tél, _Transient Chaos: Complex Dynamics on Finite-Time Scales_ , Springer, New York (2010).
  * (4) T. Tél, Chaos 25, 097619 (2015).
  * (5) J. Aguirre, J. C. Vallejo, and M. A. F. Sanjuán, Phys. Rev. E 64, 066208 (2001).
  * (6) E. Ott, _Chaos in Dynamical Systems_ , Cambridge University Press, New York (1993).
  * (7) Y.-D. Lin, A. M. Barr, L. E. Reichl, and C. Jung, Phys. Rev. E 87, 012917 (2013).
  * (8) A. Daitche and T. Tél, New J. Phys. 16, 073008 (2014).
  * (9) J. Q. Toledo-Marín and G. G. Naumis, Phys. Rev. E 97, 042106 (2018).
  * (10) J. F. Navarro, Sci. Rep. 9, 13174 (2019).
  * (11) W. Bauer and G. F. Bertsch, Phys. Rev. Lett. 65, 2213 (1990).
  * (12) Y.-T. Lau, J. M. Finn, and E. Ott, Phys. Rev. Lett. 66, 978 (1991).
  * (13) H. Alt, H. D. Gräf, H. L. Harney, R. Hofferbert, H. Rehfeld, A. Richter, and P. Schardt, Phys. Rev. E 53 2217 (1996).
  * (14) V. B. Kokshenev and M. C. Nemes, Physica A 275, 70–77 (2000).
  * (15) L. A. Bunimovich and C. P. Dettmann, Phys. Rev. Lett. 94, 100201 (2005).
  * (16) H. J. Zhao and M. L. Du. Phys. Rev. E 76, 027201 (2007).
  * (17) C. C. Moore, PNAS 112, 7, 1907–1911 (2015).
  * (18) J. L. Lebowitz and O. Penrose, Physics Today 26, 2, 23 (1973).
  * (19) V. L. Berdichevsky, J. Appl. Math. Mech. 52, 738 (1988).
  * (20) C. F. Hillermeier, R. Blümel, and U. Smilansky, Phys. Rev. A 45, 3486 (1992).
  * (21) Y.-C. Lai, M. Ding, C. Grebogi, and R. Blümel, Phys. Rev. A 46, 4661 (1992).
  * (22) R. S. MacKay, Renormalisation in Area-Preserving Maps, Advanced Series in Nonlinear Dynamics 6, World Scientific (1992).
  * (23) G. Contopoulos and M. Harsoula, Internat. J. Bifur. Chaos 20, 2005–2043 (2010).
  * (24) L. A. Bunimovich and L. V. Vela-Arevalo, Chaos 22, 026103 (2012).
  * (25) Z. Zheng, G. Hu, and J. Zhang, Phys. Rev. E 52, 3440 (1995).
  * (26) V. L. Berdichevsky and M. V. Alberti, Phys. Rev. A 44, 858 (1991).
  * (27) E. G. Altmann and T. Tél, Phys. Rev. Lett. 100, 174101 (2008).
  * (28) M. Hénon and C. Heiles, Astron. J. 69, 73 (1964).
  * (29) G. Contopoulos, Astron. Astrophys. 231, 41 (1990).
  * (30) R. L. Burden, J. D. Faires, and A. M. Burden, _Numerical Analysis_ , Tenth Edition, Cengage Learning (2015).
  * (31) I. V. Sideris, Phys. Rev. E 73, 066217 (2006).
  * (32) A. R. Nieto, E. E. Zotos, J. M. Seoane, and M. A. F. Sanjuán, Nonlinear Dyn. 99, 3029–3039 (2020).
  * (33) P. Duarte, Ann. Inst. H. Poincaré Anal. Non Linéaire 11 (4) 359–409 (1994).
  * (34) R. Barrio, F. Blesa, and S. Serrano, New J. Phys. 11, 053004 (2009).
  * (35) R. Barrio and D. Wilczak, Nonlinear Dyn. 102, 403–416 (2020).
  * (36) N. Miguel, C. Simó, and A. Vieiro, Found. Comput. Math. 15, 89–123 (2015).
  * (37) E. Ott and T. Tél, Chaos 3, 417 (1993).
  * (38) Yu. L. Bolotin, V. A. Cherkaskiy, and G. I. Ivashkevych, Phys. Lett. A 372, 4080–4085 (2008).
  * (39) Yu. L. Bolotin, V. A. Cherkaskiy, G. I. Ivashkevych, and A. I. Kirdin, Ukr. J. Phys. 55, 841–850 (2010).
  * (40) M. C. Gutzwiller, _Chaos in Classical and Quantum Mechanics_ , Springer, New York (1990).
  * (41) B. L. Lan and F. Borondo, Phys. Rev. E 83, 036201 (2011).
  * (42) S. Chanda and P. Guha, Int. J. Geom. Methods Mod. Phys. 15, 1850062 (2018).
  * (43) T. Kovács, G. Bene, and T. Tél. Mon. Notices Royal Astron. Soc. 414, 2275 (2011).
  * (44) M. Calura, P. Fortini, and E. Montanari, Phys. Rev. D 56, 4782 (1997).
  * (45) J. D. Bernal, J. M. Seoane, and M. A. F. Sanjuán, Phys. Rev. E 95, 032205 (2017).
  * (46) J. D. Bernal, J. M. Seoane, and M. A. F. Sanjuán, Phys. Rev. E 97, 042214 (2018).
  * (47) D. S. Fernández, Á. G. López, J. M. Seoane, and M. A. F. Sanjuán, Phys. Rev. E 101, 062212 (2020).
  * (48) G. Barton, _Introduction to the Relativity Principle: Particles and Plane Waves_ , John Wiley & Sons (1999).
  * (49) E. E. Zotos, Nonlinear Dyn. 82, 357–398 (2015).
  * (50) R. Klages, S. S. G. Gallegos, J. Solanpää, M. Sarvilahti, and E. Räsänen, Phys. Rev. Lett. 122, 064102 (2019).