# Tracking Blobs in the Turbulent Edge Plasma of a Tokamak Fusion Device

Woonghee Han  MIT Plasma Science and Fusion Center, Cambridge, Massachusetts
02139, USA<EMAIL_ADDRESS>Randall A. Pietersen  MIT, Civil and
Environmental Engineering, Cambridge, Massachusetts 02139, USA  Rafael
Villamor-Lora  MIT, Civil and Environmental Engineering, Cambridge,
Massachusetts 02139, USA  Matthew Beveridge  MIT, Electrical Engineering and
Computer Science, Cambridge, Massachusetts 02142, USA  Nicola Offeddu  École
Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SPC), CH-1015
Lausanne, Switzerland  Theodore Golfinopoulos  MIT Plasma Science and Fusion
Center, Cambridge, Massachusetts 02139, USA  Christian Theiler  École
Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SPC), CH-1015
Lausanne, Switzerland  James L. Terry  MIT Plasma Science and Fusion Center,
Cambridge, Massachusetts 02139, USA  Earl S. Marmar  MIT Plasma Science and
Fusion Center, Cambridge, Massachusetts 02139, USA  Iddo Drori  MIT,
Electrical Engineering and Computer Science, Cambridge, Massachusetts 02142,
USA

###### Abstract

The analysis of turbulence in plasmas is fundamental in fusion research.
Despite extensive progress in theoretical modeling in the past 15 years, we
still lack a complete and consistent understanding of turbulence in magnetic
confinement devices, such as tokamaks. Experimental studies are challenging
due to the diverse processes that drive the high-speed dynamics of turbulent
phenomena. This work presents a novel application of motion tracking to
identify and track turbulent filaments in fusion plasmas, called blobs, in a
high-frequency video obtained from Gas Puff Imaging diagnostics. We compare
four baseline methods (RAFT, Mask R-CNN, GMA, and Flow Walk) trained on
synthetic data and then test on synthetic and real-world data obtained from
plasmas in the Tokamak à Configuration Variable (TCV). The blob regime
identified from an analysis of blob trajectories agrees with state-of-the-art
conditional averaging methods for each of the baseline methods employed,
giving confidence in the accuracy of these techniques. High entry barriers
traditionally limit tokamak plasma research to a small community of
researchers in the field. By making a dataset and benchmark publicly
available, we hope to open the field to a broad community in science and
engineering.

## Introduction

Due to the enormous quantity of energy released by the fusion reaction, the
virtually inexhaustible fuel supply on earth, and its carbon-free nature,
nuclear fusion is a highly desirable energy source with the potential to help
reduce the adverse effects of climate change. Fusion research is, therefore,
an ongoing worldwide effort [1, 2].

In order to approach the conditions necessary for sufficient fusion
reactivity, the fuel – a mixture of heavy isotopes of hydrogen – must be
raised to extremely high temperatures, above 100 million degrees Celsius [3,
1] – for comparison, the core of the sun is roughly 15 million degrees
Celsius. Under these conditions, the fuel, like all stars, is in the plasma
state and must be isolated from material surfaces. Several confinement schemes
have been explored over the past 70 years . Of these, the tokamak device, a
scheme first developed in the 1950s, is the best-performing fusion reactor
design concept to date . It uses powerful magnetic fields of several to over
10 Tesla to confine the hot plasma – for comparison, this is several times the
field strength of magnetic resonance imaging machines (MRIs). The Tokamak à
Configuration Variable (TCV) [4], sited in Lausanne, Switzerland and shown in
Figure 1, is an example of such a device and provides the data presented here.

Figure 1: The schematic of the Tokamak à Configuration Variable (TCV) (left)
and its interior (right). Credits: EPFL (left) and A. Herzog, EPFL (right).

The research addressed in this paper involves phenomena that occur around the
boundary of the magnetically confined plasma within TCV. The boundary is where
the magnetic field-line geometry transitions from being “closed" to “open
."The “closed" region is where the field lines do not intersect material
surfaces, forming closed flux surfaces. The “open" region is where the field
lines ultimately intersect material surfaces, resulting in a rapid loss of the
particles and energy that reach those field lines. Around this boundary
(called the Last Closed Flux Surface or LCFS) is a region of enhanced
turbulent transport across the field lines (“perpendicular" transport).
Therefore, the transport occurring here is of great interest because of its
role in plasma confinement and dealing with the exhaust heat and particle
loads on the intersected material surfaces.

Figure 2: Cross-section of a plasma in the Tokamak à Configuration Variable
(TCV) with the locations of Gas Puff Imaging (GPI) views near the Last Closed
Flux Surface (LCFS) (left). Snapshot of experimental GPI data capturing a blob
on the right-hand side moving radially outward (middle). Here empty (white)
spots correspond to dead GPI views. The brightness level is coded in the color
bar to the right, with low as blue and high as yellow. Snapshot of synthetic
data capturing a blob moving radially outward (right). A Gaussian ellipse
represents the blob with a major and minor axis marked by perpendicular black
lines.

Researchers use a technique called Gas-Puff Imaging (GPI) [5] to visualize the
phenomena occurring at and around the plasma boundary in both space and time.
A small amount of neutral gas is locally injected into the region of interest.
The visible-light emissions resulting from the interaction of the plasma with
this gas cloud are captured along sight lines that are tangential to the local
magnetic field. Analysis of the time sequences of the images produced by this
technique is the primary subject of this work. At TCV, the 2D GPI data are
collected at frame rates of $\sim$ 2 MHz, significantly faster than the
turbulence timescales of interest. The images represent what is occurring
perpendicular to the direction of the local magnetic field lines. This is
illustrated schematically in Figure 2 (left), showing cross-sections of plasma
flux surfaces and the TCV vacuum vessel. (The center-line of this toroidally
symmetric plasma geometry is to the left of the figure.) A $\sim$50$\times$40
mm cross-section of plasma, spanning the LCFS at the outer edge of the plasma,
is imaged with a 12$\times$10 pixel array. The images’ motion from left to
right increases the radial coordinate and moves away from the closed flux
surfaces. The most prominent features in the image sequences are 1) a bright
band of emission that is aligned with the LCFS and 2) the expulsion of
features from the LCFS region that are brighter than their surroundings and
the subsequent motion of these features. In the literature, these features are
commonly called “filaments” (because they are actually extended in the third
dimension, along the magnetic field lines) or “blobs” (because of their
appearance in cross-section in the images) [6]. These filaments are strong
perturbations (of order 1) in the local plasma conditions, with elevated
electron temperature and density [7, 8]. They have typical auto-correlation
times of roughly 10 $\mu$s and, therefore, can exist as distinct features for
multiple frames, but their shape and intensity can gradually change frame-to-
frame. In this work, we are mainly interested in the detection and tracking of
these “blob” features.

The primary problem we are addressing in this paper is how to process these
image sequences in ways that will ultimately allow estimation of the
contributions of this blob transport to the heat and particle fluxes leaving
the confined plasma. We are also addressing how to process the images to
facilitate continued and more detailed comparisons of the experimental
observations with theoretical models of blob dynamics. The theoretical models
predict scalings for the relation between blob size and radial velocity
depending on key plasma parameters, such as local temperature and
collisionality [9, 10].

This primary problem is essential since the expulsion and transport of these
objects represent particle and energy loss from the confined plasma. This has
enormous consequences for the interaction of the hot plasma with the plasma-
facing material surfaces that are in place to accept the exhaust heat and
particles. Parallel heat fluxes of hundreds of MW/m2 are obtained in the
narrow exhaust channels of present-day devices [11], and the parallel heat
fluxes will be even higher in the near-future net-energy-gain tokamaks like
SPARC [12] and ITER [13]. The radially-outward motion of the blobs can broaden
this exhaust channel, thereby reducing the peak heat- and particle fluxes, the
desired circumstance. However, the same motion can also increase an undesired
plasma interaction with the other plasma-facing components. These
considerations constitute the highest-level motivations for an accurate
assessment of the blobs’ occurrence, size, and motion. Combining knowledge of
the occurrence frequency, size, and radial velocity with temperature and
density within the blobs allows quantitative estimations of radial particle-
and energy fluxes. The additional problem that we are addressing, facilitating
closer experimental comparisons with theoretical models of blob transport, is
essential for experimental validation tests of these models. These tests of
the models (e.g., [14]) are made more valuable with accurate determinations of
the size, shape, and motion of blobs.

Traditional approaches to blob analysis include a family of conditional
averaging methods, custom-made workflows that track high signal regions, and
spatio-temporal cross-correlation techniques [14, 5, 15]. The conditional
averaging methods and the cross-correlation techniques have the limitations of
only providing averaged characteristics of blobs. The custom-made workflows
are non-standardized and are not benchmarked. The data sets that we are
working with typically include $\sim 10^{5}$ sequential images, having been
recorded over a time-subset shorter than the $\sim 1$ s plasma pulse on TCV.
30 such plasma discharges can be produced per each working day at TCV; hence
this data throughput is too high for human by-eye analysis. Machine learning
allows accurate and efficient analysis of turbulence characteristics utilizing
these extensive datasets.

This work presents a novel application of four well-known, standardized, and
benchmarked tracking methods to track blobs in GPI images. These methods are
trained to reproduce identification of blobs by humans as close as possible
since human subjectivity isolates blobs from the non-blobs in the GPI images.
We evaluate their performance by comparing them to the brightness-threshold
contours and human detection/tracking using a limited set of real-data images.
In Methods we describe the four models: 1) optical flow detection using
Recurrent All-pairs Field Transforms (RAFT) [16]; 2) mask detection using Mask
Region-based Convolutional Neural Network (Mask R-CNN) [17] in combination
with Bayesian optimization [18] for tuning the hyperparameters in the
training; 3) optical flow detection using Global Motion Aggregation (GMA)
[19]; and 4) optical flow detection using Flow Walk [20, 21]. In Methods we
also define the performance metrics and the workflow. The performance
comparisons are given in Results, as is a description of the results obtained
by human detection, which shows some ambiguity in what constitutes a blob in
real experimental data. Finally, the trained models are applied to an active
plasma-physics research topic, i.e., finding the regime of the blob dynamics,
as described in Myra _et al._ [9, 10], and the models’ results are compared
with the result from a Conditional Average Sampling (CAS) method that has been
used previously for this purpose [14]. The work is summarized in Conclusions.

## Methods

We train our models on synthetic data and test them on synthetic and real
data. Synthetic data affords several advantages for training a model. First,
our models use optical flow vectors (RAFT, GMA, Flow Walk) or masks (Mask
R-CNN) of the blob objects as ground-truth features that humans cannot
adequately label. The user has perfect knowledge of and control over these
features and noise in the synthetic dataset. By contrast, feature recognition
in available real data may be prohibitively costly to process, hard to
reproduce, poor quality, and not representative of general phenomenology.
Second, synthetic datasets of arbitrary size and complexity can be made more
generic than a small, real, vetted dataset, such that a model may be trained
to recognize features that do not appear in available real data but might
appear more broadly. There are a few edge cases in the real data to be covered
in training. Events such as splitting or merging blobs often occur but not
often in real data. We further iterate upon the generation of synthetic data
by adding, removing, or emphasizing characteristics like these events.
Moreover, synthetic data is more efficient in acquisition and storage than
real data. Real GPI data is costly to acquire; blob dynamics vary with the
plasma condition, and we need to run many experiments to generate a
significant variety of blob dynamics in the real world. Also, the algorithm to
create the synthetic data is more compact than a large real dataset and is
easily adapted to other situations where the effort expended in processing the
real dataset cannot be reused. Therefore, by using synthetic data, we can have
inexpensive, variable, and unlimited data generation [22, 23, 24].

An essential goal of this contribution is to solicit greater participation in
fusion research from the broader machine-learning community. Toward this end,
we make our synthetic training dataset and a real dataset available to
benchmark performance against other models. We hope this may inspire readers
not only to evaluate their model performance against this fusion-relevant task
but to seek out and engage more generally to help solve critical problems in
the field of fusion energy.

### Dataset

##### Real Data

Each experimental GPI video is a length $t$ series of grayscale images with
$12\times 10$ pixels, where $t$ is the number of frames in the video. A
snapshot of an experimental GPI video is shown in Figure 2 (middle), which
shows a blob (as a bright spot on the bottom-right), and the Last Closed Flux
Surface (LCFS), which is also approximately the position of a shear layer,
across which the vertical plasma background flow reverses direction. Before
inputting the experimental GPI data into the models, we standardize the
amplitude ranges between 0 and 1 and remaps the images to a finer spatial
grid. We subtract the mean brightness from the individual pixel brightnesses
and divide the result by the standard deviation. It is also clipped by
$\frac{1}{n}\sum_{t=1}^{n}C_{min}\left(t\right)$ and
$\frac{1}{n}\sum_{t=1}^{n}C_{max}\left(t\right)$, where $n$ is the total
number of frames, and $C_{min}\left(t\right)$ and $C_{max}\left(t\right)$ are
the minimum and maximum (standardized) brightnesses of the frame at time $t$,
respectively. In addition, although the GPI sensors output data with $12\times
10$ pixels per frame, the images become better interpreted by the models when
we upsample their resolution, in the present case, to $256\times 256$ pixels
per frame. We use radial basis function (RBF) interpolation with a cubic
function. This degree of upsampling allows adequate detail during analysis
without being too computationally demanding. We design the models to receive
inputs and produce outputs at this resolution.

##### Synthetic Data

We maintain the same upsampling convention as real data by generating
synthetic data for training at $256\times 256$ pixels per frame. We
approximate blobs in the synthetic data as ellipses, and a blob’s size and
speed at each moment are randomized so that a random process forms the
trajectory. Here, the ellipse boundary is chosen to be the full width at half
maximum (FWHM) of the brightness since the blob sizes are conventionally
estimated as the FWHM of the density perturbation [25] which is approximately
proportional to the brightness in the GPI data. Besides blobs, we simulate
background flows as slowly moving, elongated ellipses. Figure 2 (right) shows
a snapshot of such a synthetic video. We also simulate complex instances in
the real GPI data for these synthetic data, such as merging and splitting
blobs. The brightness of the synthetic data ranges between 0 and 1. We save
the mask of blobs and the velocities at pixels in blobs as an optical flow
image for every frame. We use these images as labels for training.

### Benchmark

Our solution combines optical flow/mask detection with Tracking-by-detection
to achieve accurate identification and tracking of the blobs. We provide and
compare four baselines for detecting blobs and measuring performance using
standard and specialized metrics.

#### Baselines

##### RAFT

Recurrent all-pairs field transforms (RAFT) [16] computes optical flow by
extracting features of pixels and building multiscale 4D correlation volumes
for all pairs of pixels. RAFT iteratively updates an optical flow field using
a recurrent unit that uses the correlation volumes.

##### Mask R-CNN with BO

We applied mask region-based convolutional neural network (mask R-CNN) [17],
which computes segmented masks of the data. Bayesian optimization is used to
find the hyperparameter values that minimize the loss of the mask R-CNN during
the training. We apply Bayesian optimization [18] with two levels of
hyperparameters: (1) learning rate, momentum, weight decay, and number of
epochs; and (2) data augmentation transformations $P_{horizontalFlip}$,
$P_{scale}$, $P_{translate}$, $P_{shear}$, $P_{rotate}$, as well as dropout
probability $P_{dropout}$ of each dropout layer in the model. The exploration
ranges for each hyperparameter are in the supplementary information. We first
optimize the group one hyperparameters and then keep those values fixed while
optimizing the second group of hyperparameters.

##### GMA

We also implemented global motion aggregation (GMA) [19] to estimate hidden
motions. It finds long-range dependencies between pixels in an image and
performs global aggregation of corresponding motion features using a
transformer model. GMA is tailored to perform well on occluded regions which
is not the case in our application, so this feature is not exploited.
Nonetheless, we implemented GMA for comparison.

##### Flow Walk

Flow Walk learns pixel trajectories with multiscale contrastive random walks
by computing the transition matrix between frames in a coarse-to-fine manner
[20, 21]. Flow Walk works well on detecting pixel-level changes of objects
with high spatial frequencies, which is not the feature that appears on our
data, as we are tracking smooth blobs. As with GMA, we implemented Flow Walk
for comparison.

#### Tracking-by-detection

After detecting blobs in each frame, temporal coherence between frames is
enforced based on a tracking-by-detection workflow (illustrated in the
supplementary information) using any of the baselines, which consists of four
steps:

  1. 1.

Object detection from the model. Given the input image sequence, the masks of
blobs are predicted for each frame by the model.

  2. 2.

Feature extraction. We select the objects to be tracked (i.e., blob objects)
by discarding predictions with scores below a threshold.

  3. 3.

Pairwise cost. We exploit the temporal coherence of the images by computing
the pairwise cost between the objects in the current and previous frames using
the VIoU cost metric.

  4. 4.

Bipartite matching. Using the cost matrix from the previous step, we assign
unique correspondence between objects with the constraint that no object
receives more than one assignment. If a new object appears in an isolated
frame (i.e., it has no correspondence either in the previous or the next
frame), it is ignored and discarded as a noise fluctuation. When a new object
appears in the current frame with no correspondence in the previous frame but
with a match in the next frame, we start a new track. By keeping track of the
active and finished tracks, we assign IDs to the blobs in the video and record
their trajectories.

#### Performance metrics

We use standard and specialized metrics to evaluate performance.

  * •

Endpoint error (EPE): The standard error measure for evaluating optical flow
is defined by $\frac{1}{n}\sum\|\vec{v}_{est}-\vec{v}_{gt}\|$, where $n$ is
the number of pixels in the image, $\vec{v}_{est}$ is the estimated flow field
and $\vec{v}_{gt}$ is the proxy ground truth flow field.

  * •

Volumetric IoU (VIoU): Object detection algorithms usually use intersection
over union (IoU) metrics to evaluate the prediction quality. For the
prediction/ground-truth pair (e.g., bounding boxes, mask), the IoU is computed
as the ratio between the intersection and union areas. While this approach is
quite robust when dealing with the mask detection of solid objects with well-
defined, sharp boundaries, IoU can mislead our application’s score when the
intersection area has a low brightness level, which is intuitively a poor
prediction. Note that the boundary definition of blob features within the
images is a function of the brightness threshold. Since this threshold may
vary from experiment to experiment, in our case, it is more appropriate to
describe the blob shape as a volume, with volumetric IoU (VIoU), using the
brightness level for the third axis (illustrated in the supplementary
information), and defined as:

$VIoU=\frac{\sum_{x,y\in A_{intersection}}B_{x,y}}{\sum_{x,y\in
A_{union}}B_{x,y}}$ (1)

where $B_{x,y}$ is the brightness at pixel located at $\left(x,y\right)$, and
$A_{intersection}$ and $A_{union}$ are, respectively, the set of pixels inside
the intersection and union of the predicted/ground-truth blob mask.

## Results

We present results in four parts: (1) training results for each model
described in Training, (2) testing scores of the trained models on both
synthetic and real GPI data shown in Testing and blob tracking, (3)
furthermore, the models are evaluated based on human-labeled blobs to
demonstrate their validity in Evaluation by human-labeled blobs, (4) the
tracking information from each model is used to identify the regime of the
blob dynamics in Identification of blob regimes to perform a task that
addresses an active plasma-physics research topic.

### Training

We train the tracking models using 30 synthetic videos, each with 200 frames.
We split the frames into training data (95 percent) and validation data (5
percent). We used a training pipeline similar to the original implementation
with default hyperparameters for RAFT and GMA. For Mask R-CNN, we find
hyperparameters by Bayesian optimization (see the supplementary information).
For Flow Walk, we implemented the training pipeline from RAFT. The training
scores based on the Performance metrics are shown in the supplementary
information, where RAFT is the best among the four models.

### Testing and blob tracking

We test the models on synthetic and real GPI data and summarize the results in
Table 1. For the real GPI data, we only compute VIoU since there are no labels
in optical flows. Among the four models, RAFT performs best on real data.

Table 1: Scores from testing for each model with corresponding score metrics, endpoint error (EPE) and volumetric IoU (VIoU). Lower EPE and higher VIoU are better. Among four models, RAFT performs best with the lowest EPE and highest VIoU. Mask R-CNN misses EPE because it is not an optical flow detection model but a mask detection model. EPE for real data is not presented because the real data does not have ground-truth optical flow velocities. |  | Testing score
---
on synthetic data
| Testing score
---
on real data
Models\Metric | EPE | VIoU | VIoU
RAFT | 2.701 | 0.783 | 0.781
Mask R-CNN | N/A | 0.638 | 0.700
GMA | 2.974 | 0.777 | 0.650
Flow Walk | 6.020 | 0.631 | 0.621
Figure 3: An example of blob tracking on experimental data. With the current
frame (left) and the next frame as input, RAFT predicts optical flows for each
pixel (middle). In the middle, the color of each pixel indicates the velocity
(speed and direction) of the optical flow predicted, corresponding to the
color palette shown at the bottom left. On the right, the red contour is drawn
for the boundary of the pixels having non-zero optical flow in the central
figure. The blue contour is the contour of $0.7\times B_{max}$ where $B_{max}$
is the maximum (standardized) brightness inside the red contour, which is used
for computing VIoU. The red contours chosen by several thresholds (VIoU
$>0.8$, $B_{max}>0.75$, lifetime $>15$ frames) are only tracked. The blob ID
(Blob 5) was assigned by Tracking-by-detection.

Figure 3 shows an example of blob tracking by RAFT on real experimental data.
RAFT takes two images as input: the current and the next frame. It then
predicts optical flows for each pixel, as shown in the middle of Figure 3. The
color of each pixel indicates the velocity of the predicted optical flow,
corresponding to the color palette shown at the bottom left. In order to
identify the mask of the blobs from the map of optical flows, which is not
segmented, the algorithm computes a set of contours of different values of the
magnitude of optical flows in each frame. We merge contours that overlap
within VIoU $>\alpha$. We find $\alpha=0.2$ suitable for our purposes. For
each contour, we compute the VIoU using the blue contour in the figure, which
is the contour of $\beta\times B_{max}$, where $B_{max}$ is the maximum
(standardized) brightness inside the red contour. Then the contours which give
a VIoU greater than $\gamma=0.8$ and $B_{max}$ greater than $\beta=0.7$ are
chosen to be tracked. The user sets the hyperparameters $\alpha$, $\beta$, and
$\gamma$ based on the GPI data, and we set them as above for the data used in
this work. In Figure 3 (right), the blob ID (Blob 5) is assigned by the
Tracking-by-detection. In order to enforce temporal coherence in blob
detection, only blobs with a lifetime longer than 15 frames (corresponding to
7.5 micro-seconds in real data) are chosen and drawn in the output video.
Unlike the optical flow prediction models, Mask R-CNN outputs segmented masks
of objects, and therefore the contour of an object is already obtained (i.e.,
we do not need to scan the values of contours). The rest of the tracking is
the same as above.

Figure 4: Examples of model prediction on the real data with human labels. The
labeler marks a dot (magenta) on the blob that he/she identifies in every
frame. As illustrated in the left figure, TP is true positive where the human-
labeled dot is contained inside the predicted blob boundary, and FN is false
negative where the human-labeled dot is not contained in any predicted blob
boundary. FP is a false positive (middle figure) where the predicted blob
boundary does not contain human labels. In the right figure, there are neither
predictions nor human labels in this frame, hence the case of true negative.

### Evaluation by human-labeled blobs

For the testing score on real data shown in Table 1, the label (blue contour
in Figure 3 (right)) is dependent on the prediction (red contour). In other
words, the label is drawn only for the predicted objects, and the model
performance is thus evaluated only for structures detected by our tracking
models (i.e., there are only true positives). In general, there is some
subjectivity in identifying blobs in real data. Therefore, we quantified how
close the machine-predicted blobs are to the human-labeled blobs selected by
different domain experts. We consider as “ground-truth" the cases that human
labelers have identified as blobs. We cover cases with false positives (the
model identified a blob where the human identified none), true negatives (did
not identify a blob where there was none), false negatives (did not identify a
blob where there was one), as well as true positives (identified a blob where
there was one), as defined in Figure 4. Each of the three domain experts
separately labeled the blobs in 3,000 frames by hand, and our blob-tracking
models are evaluated against these human-labeled experimental data based on F1
score, False Discovery Rate (FDR), and accuracy, as shown in Figure 5. These
are the average per-frame scores (i.e., the average across the frames), and we
did not use the score across all frames, which can be dominated by outlier
frames that may contain many blobs. Figure 6 displays the corresponding
confusion matrices. In this result, RAFT, Mask R-CNN, and GMA achieved high
accuracy (0.807, 0.813, and 0.740 on average, respectively), while Flow Walk
was less accurate (0.611 on average). Here, the accuracy of 0.611 in Flow Walk
is seemingly high, misleading because Flow Walk gave few predictions (low TP
and FP in Figure 6). This is because the data is skewed to true negatives as
many frames have no blobs, which is seen from the high true negatives of
confusion matrices in Figure 6. Thus, accuracy is not the best metric for the
data used. F1 score and FDR are more suitable for our purposes because they
are independent of true negatives. Indeed, other scores of Flow Walk are as
expected; the F1 score is low (0.036 on average) and the FDR is high (0.645 on
average). RAFT and Mask R-CNN show decently high F1 scores and low FDR. GMA
underperformed RAFT and Mask R-CNN in all metrics, but the scores are fairly
good. These observations are consistent across the comparisons with the three
labelers’ results. Overall, blobs predicted by three models (RAFT, Mask R-CNN,
and GMA) are quite similar to human-labeled blobs, exhibiting fairly good
scores in Figure 5.

Figure 5: F1 score, false discovery rate (FDR), and accuracy of the four
methods (RAFT, Mask R-CNN, GMA, and Flow Walk) on the real GPI data with 3,000
frames hand-labeled by three domain experts (A, B, and C). Metrics shown are
the average per-frame scores (i.e., the average across the frames). Figure 6:
Confusion matrices of the four methods (RAFT, Mask R-CNN, GMA, and Flow Walk),
for the real GPI data with 3,000 frames hand-labeled by three domain experts
(A, B, and C).

The trained models’ average running times per frame (without post-processing)
are milliseconds. Specifically, 104, 50, 109, and 23 milliseconds for RAFT,
Mask R-CNN, GMA, and Flow Walk, respectively, using 4 GPUs. For RAFT and GMA,
there are 12 iterations for each frame following the default setup, and we run
1 iteration for Mask R-CNN and Flow Walk. This indicates that $\sim 3$ hours
(RAFT, GMA) and $\sim 1$ hours (Mask R-CNN, Flow Walk) would be needed to
process $10^{5}$ images. This can become faster if we reduce the number of
iterations or downsample the input image at the expense of the model’s
performance. The post-processing takes seconds, precisely 12, 11, and 9
seconds per frame for RAFT, GMA, and Flow Walk, respectively, and $\sim 300$
hours for $10^{5}$ images. This is not short because of the search for the
contour value to be drawn from the optical flow for the blob contour
prediction. For Mask R-CNN, the per-frame post-processing time is 1 second,
much faster than the other three optical flow models. This is because Mask
R-CNN gives segmented masks, and there is no need to find a contour value.
Therefore, RAFT, GMA, and Flow Walk can be run for a small subset of the data
($\sim$1,000 frames), and Mask R-CNN can be run for the data with a larger
number of frames. Despite the long computation time, the post-processing
provides information of blobs for every frame, whereas CAS only provides
average results in a relatively shorter time.

### Identification of blob regimes

The trained models can now estimate various blob parameters, such as their
size, speed, and occurrence frequency in real GPI data. Here, we use such
information to identify the regime of the blob dynamics for two different
plasma conditions. The regime is identified based on location in the diagram
in Figure 7. Here, $\Theta$ (a normalized blob size) and $\Lambda$ (a
normalized plasma collisionality, which is proportional to the electron-ion
collision frequency) are defined as in [9, 10]

$\Theta=\hat{a}^{5/2}\quad\mathrm{and}\quad\Lambda=1.7\times
10^{-18}\frac{n_{e}L_{\|}}{T_{e}^{2}}$ (2)

where

$\hat{a}=\frac{a_{b}R^{1/5}}{L_{\|}^{2/5}\rho_{s}^{4/5}}\quad\mathrm{and}\quad\hat{v}=\frac{v_{R}}{c_{s}\left(2L_{\|}\frac{\rho_{s}^{2}}{R^{3}}\right)^{1/5}}$
(3)

$\hat{a}$ and $\hat{v}$ are the blob’s normalized radius and radial speed,
respectively. Here, $n_{e}$ and $T_{e}$ are the local electron density and
temperature, respectively, $L_{\|}$ is the parallel connection length, $a_{b}$
is the blob radius, $R$ is the major radius of the tokamak, $\rho_{s}$ is the
ion sound Larmor radius, $v_{R}$ is the radial speed of the blob, and $c_{s}$
is the sound speed. Our blob-tracking models allow us to estimate $a_{b}$ and
$v_{R}$. Other plasma parameters are measured using other diagnostics. There
are four regimes in the diagram, named “resistive ballooning" (RB), “resistive
X-point" (RX), “connected ideal-interchange" ($C_{i}$), and “connected sheath"
($C_{s}$). Theory predicts different scaling relationships between the
normalized radial speed ($\hat{v}$) and the radius ($\hat{a}$) of the blob
depending upon regime (see Figure 7). We have evaluated $\Theta$ and $\Lambda$
for two different plasma conditions (Plasma 1 and Plasma 2) using a
traditional method (CAS) and our four blob-tracking models (RAFT, Mask R-CNN,
GMA, and Flow Walk). This evaluation locates the two plasmas within the
theory-defined regimes, as shown in Figure 7. We identify plasma 1 and 2 as
being within $C_{s}$ and RX regimes, respectively, unanimously by four methods
(CAS, RAFT, Mask R-CNN, and GMA) except Flow Walk. The closeness of the data
points between the traditional method and the blob-tracking models
demonstrates the validity of the machine learning approach in blob-tracking
applied to an important research investigation. Note that the bars around the
centroids for the blob-tracking models are not error bars but rather the
spread in actual blob statistics. This is an advantage of the blob-tracking
methods since they allow statistics derived from individual blob measurements,
whereas CAS can only yield average results. Flow Walk’s large spread is due to
its poor prediction performance, as shown previously in Figure 5.

Figure 7: $\Theta$—$\Lambda$ diagram where $\Theta$ and $\Lambda$ are defined
in Equation 2. The diagram is divided into four blob regimes (RB, RX, $C_{i}$,
and $C_{s}$) depending on the relation between the normalized speed
($\hat{v}$) and radius ($\hat{a}$) of the blobs. Data from two different
plasmas (Plasma 1 and Plasma 2) are interrogated by four blob-tracking models
(RAFT, Mask R-CNN, GMA, and Flow Walk) and the traditional CAS method, with
the locations from each method indicated in the diagram. For these data, the
factor $\epsilon_{x}$ is 0.5. The bars correspond to the standard deviation of
$\log{\Theta}$ and $\log{\Lambda}$.

## Conclusions

We introduce a dataset and benchmark for tracking turbulent structures, called
blobs, in Gas Puff Imaging data from the boundary region of magnetically
confined tokamak plasmas created in TCV. We experimented with four baseline
models based on RAFT, Mask R-CNN, GMA, and Flow Walk, trained on synthetic
data, and tested on synthetic and real data. The synthetic data mimic the
realistic movement and evolution of blobs found in actual GPI videos of the
TCV boundary plasma. We optimize the training hyperparameters by Bayesian
optimization. The trained models show high VIoU testing scores on both
synthetic and real data and high accuracy scores on human-labeled real data.
As a demonstration of the validity of our tracking models, we achieve (nearly)
unanimous identification of the regime of the blob dynamics for two different
plasma conditions using blob statistics estimated by tracking models and the
traditional CAS method. At least two of these models (RAFT and Mask R-CNN) are
reliable tools for identifying and tracking blobs, which allow the estimation
of various measures of blob dynamics connected to the levels and effects of
turbulence on the edge of tokamak plasmas. This characterization of edge
turbulence furthers our understanding of transport processes in the plasma
boundary, contributing to essential knowledge for the practical generation of
fusion energy.

## Data availability

In the spirit of reproducible research, we make our data, models, and code
publicly available at https://github.com/harryh5427/GPI-blob-tracking.

## Acknowledgements

The support from the US Department of Energy, Fusion Energy Sciences, awards
DE-SC0014264 and DE-SC0020327, are gratefully acknowledged. This work was
supported in part by the Swiss National Science Foundation. Also, this work
has been carried out within the framework of the EUROfusion Consortium, funded
by the European Union via the Euratom Research and Training Programme (Grant
Agreement No 101052200 -EUROfusion). Views and opinions expressed are, however
those of the author(s) only and do not necessarily reflect those of the
European Union or the European Commission. Neither the European Union nor the
European Commission can be held responsible for them.

## Author contributions

W.H. constructed the initial concept of the study, curated the dataset, and
wrote the code. W.H. and I.D. defined the benchmark and reviewed the code. All
authors provided ideas and analyzed the results. W.H., R.A.P., R.V.L., N.O.,
T.G., C.T., J.T., and I.D. contributed to the writing and all authors reviewed
the manuscript. W.H., I.D., and R.V.L. created figures.

## References

  * [1] Bécoulet, A., Butler, E. & Whyte, D. G.  _Star Power: ITER and the International Quest for Fusion Energy_ (The MIT Press, 2021).
  * [2] Ball, P.  The chase for fusion energy.  _Nature_ (2021).  https://www.nature.com/immersive/d41586-021-03401-w/index.html.
  * [3] Tester, J. _et al._ _Sustainable Energy: Choosing Among Options_ (MIT Press, 2012).
  * [4] Reimerdes, H. _et al._ Overview of the TCV tokamak experimental programme.  _Nuclear Fusion_ 62, 042018, DOI: 10.1088/1741-4326/ac369b (2022).
  * [5] Zweben, S. J., Terry, J. L., Stotler, D. P. & Maqueda, R. J.  Invited review article: Gas puff imaging diagnostics of edge plasma turbulence in magnetic fusion devices.  _Review of Scientific Instruments_ 88, 041101 (2017).
  * [6] D’Ippolito, D. A., Myra, J. R. & Zweben, S. J.  Convective transport by intermittent blob-filaments: Comparison of theory and experiment.  _Physics of Plasmas_ 18, 060501, DOI: 10.1063/1.3594609 (2011).  https://doi.org/10.1063/1.3594609.
  * [7] Agostini, M. _et al._ Fast thermal helium beam diagnostic for measurements of edge electron profiles and fluctuations.  _Review of Scientific Instruments_ 86, 123513, DOI: 10.1063/1.4939003 (2015).  https://doi.org/10.1063/1.4939003.
  * [8] Kube, R. _et al._ Intermittent electron density and temperature fluctuations and associated fluxes in the alcator c-mod scrape-off layer.  _Plasma Physics and Controlled Fusion_ 60, 065002, DOI: 10.1088/1361-6587/aab726 (2018).
  * [9] Myra, J. R. _et al._ Blob birth and transport in the tokamak edge plasma: Analysis of imaging data.  _Physics of Plasmas_ 13, 092509, DOI: 10.1063/1.2355668 (2006).  https://doi.org/10.1063/1.2355668.
  * [10] Myra, J. R., Russell, D. A. & D’Ippolito, D. A.  Collisionality and magnetic geometry effects on tokamak edge turbulent transport. i. a two-region model with application to blobs.  _Physics of Plasmas_ 13, 112502, DOI: 10.1063/1.2364858 (2006).  https://doi.org/10.1063/1.2364858.
  * [11] Brunner, D., LaBombard, B., Kuang, A. & Terry, J.  High-resolution heat flux width measurements at reactor-level magnetic fields and observation of a unified width scaling across confinement regimes in the alcator c-mod tokamak.  _Nuclear Fusion_ 58, 094002, DOI: 10.1088/1741-4326/aad0d6 (2018).
  * [12] Kuang, A. Q. _et al._ Divertor heat flux challenge and mitigation in sparc.  _Journal of Plasma Physics_ 86, 865860505, DOI: 10.1017/S0022377820001117 (2020).
  * [13] Goldston, R.  Theoretical aspects and practical implications of the heuristic drift sol model.  _Journal of Nuclear Materials_ 463, 397–400, DOI: https://doi.org/10.1016/j.jnucmat.2014.10.080 (2015).  PLASMA-SURFACE INTERACTIONS 21.
  * [14] Offeddu, N. _et al._ Cross-field and parallel dynamics of sol filaments in tcv.  _Nuclear Fusion_ (2022).  In press.
  * [15] Häcker, R., Fuchert, G., Carralero, D. & Manz, P.  Estimation of the plasma blob occurrence rate.  _Physics of Plasmas_ 25, 012315, DOI: 10.1063/1.5008301 (2018).  https://doi.org/10.1063/1.5008301.
  * [16] Teed, Z. & Deng, J.  RAFT: Recurrent all-pairs field transforms for optical flow.  In _European conference on computer vision_ , 402–419 (Springer, 2020).
  * [17] He, K., Gkioxari, G., Dollár, P. & Girshick, R.  Mask R-CNN.  In _Proceedings of the IEEE international conference on computer vision_ , 2961–2969 (2017).
  * [18] Balandat, M. _et al._ Botorch: A framework for efficient monte-carlo bayesian optimization.  _Advances in Neural Information Processing Systems_ 2020-December (2020).  Funding Information: Andrew Gordon Wilson is supported by NSF I-DISRE 193471, NIH R01 DA048764-01A1, NSF IIS-1910266, and NSF 1922658 NRT-HDR: FUTURE Foundations, Translation, and Responsibility for Data Science. Publisher Copyright: © 2020 Neural information processing systems foundation. All rights reserved.; 34th Conference on Neural Information Processing Systems, NeurIPS 2020 ; Conference date: 06-12-2020 Through 12-12-2020.
  * [19] Jiang, S., Campbell, D., Lu, Y., Li, H. & Hartley, R.  Learning to estimate hidden motions with global motion aggregation.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 9772–9781 (2021).
  * [20] Bian, Z., Jabri, A., Efros, A. A. & Owens, A.  Learning pixel trajectories with multiscale contrastive random walks.  _CoRR_ abs/2201.08379 (2022).  2201.08379.
  * [21] Bian, Z., Jabri, A., Efros, A. A. & Owens, A.  Learning pixel trajectories with multiscale contrastive random walks.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ (2022).
  * [22] Nikolenko, S. I.  _Synthetic data for deep learning_ , vol. 174 (Springer, 2021).
  * [23] Andrews, C., Sirazitdinova, E., Hedges, D., Robinson, M. & Kulkarni, S.  Tutorial on creating and using synthetic data for computer vision applications (2022).
  * [24] Wood, E., Fidler, S., Urtasun, R. & Laserson, J.  Workshop on machine learning with synthetic data (2022).
  * [25] Fuchert, G., Carralero, D., Manz, P., Stroth, U. & and, E. W.  Towards a quantitative prediction of the blob detection rate.  _Plasma Physics and Controlled Fusion_ 58, 054006, DOI: 10.1088/0741-3335/58/5/054006 (2016).

## Appendix

We provide additional details on the methods and results. Figure 8 illustrates
Tracking-by-detection workflow, which is applied to the consecutive frames
containing blob contours predicted by the model and assigns blob IDs based on
closeness to the blobs in the previous frame. Figure 9 shows an example of
computing the VIoU between the contours of a prediction (red) and a label
(black). For the training of Mask R-CNN, we explore the hyperparameters in
Table 2 by Bayesian optimization and find the optimal values shown in Table 3.
The scores of the four models on training and validation synthetic data are
shown in Table 4 with corresponding score metrics (EPE and VIoU). The EPE
scores of Mask R-CNN are not applicable because it is not an optical flow
detection model but a mask detection model.

Figure 8: Tracking-by-detection consisting of four steps: (1) Object detection within each frame using a pre-trained model. Here the mask R-CNN produces a bounding box and a mask for each candidate object. (2) Extraction of features of interest (i.e., masks of blobs). (3) Computation of pairwise costs between objects in the current and previous frame. (4) Bipartite matching between objects assigning unique correspondence. Figure 9: Volumetric Intersection over Union (VIoU) computed for a prediction on an image containing a blob (top) using the volume under the profile of the brightness (bottom). Table 2: Exploration ranges of hyperparameters used in Bayesian optimization. BO level 1 hyperparameters: learning rate initial value, learning rate reduction period, learning rate decay factor, momentum, weight decay and number of training epochs. BO level 2 hyperparameters: probabilities for horizontal flip, scaling, translation, shearing, rotation, and dropout. BO Lv. 1 |  LR initial val |  LR reduction period |  LR decay | Momentum |  Weight decay | $n_{epochs}$
---|---|---|---|---|---|---
| Exploration
---
range
1e-6—5e-2 | 1—40 | 0.1—1.0 | 0.0—1.0 | 0.0—1e-3 | 2—40
BO Lv. 2 | $P_{horizontalFlip}$ | $P_{scale}$ | $P_{translate}$ | $P_{shear}$ | $P_{rotate}$ | $P_{dropout}$
| Exploration
---
range
0.0—1.0 | 0.0—1.0 | 0.0—1.0 | 0.0—1.0 | 0.0—1.0 | 0.0—0.5
Table 3: Optimal values of hyperparameters found from Bayesian optimization. BO Lv. 1 |  LR initial val |  LR reduction period |  LR decay | Momentum |  Weight decay | $n_{epochs}$
---|---|---|---|---|---|---
| Optimal
---
value
0.050 | 12 | 0.402 | 0.296 | 2.729e-5 | 28
BO Lv. 2 | $P_{horizontalFlip}$ | $P_{scale}$ | $P_{translate}$ | $P_{shear}$ | $P_{rotate}$ | $P_{dropout}$
| Optimal
---
value
0.341 | 0.560 | 0.261 | 0.564 | 0.368 | 1.283e-17
Table 4: Scores from training and validation for each model with corresponding score metrics, endpoint error (EPE) and volumetric IoU (VIoU). Mask R-CNN misses EPE because it is not an optical flow detection model but a mask detection model. |  | Training score
---
on synthetic data
| Validation score
---
on synthetic data
Models\Metric | EPE | VIoU | EPE | VIoU
RAFT | 0.161 | 0.869 | 0.338 | 0.879
GMA | 0.160 | 0.865 | 0.298 | 0.862
Mask R-CNN | N/A | 0.705 | N/A | 0.704
Flow Walk | 0.946 | 0.835 | 0.924 | 0.830


††thanks<EMAIL_ADDRESS>
# Light chaotic dynamics in the transformation from curved to flat surfaces

Chenni Xu1,2 Itzhack Dana1 Li-Gang Wang2 Patrick Sebbah1 1 Department of
Physics, The Jack and Pearl Resnick Institute for Advanced Technology, Bar-
Ilan University, Ramat-Gan 5290002, Israel  2 Zhejiang Province Key Laboratory
of Quantum Technology and Device, Department of Physics, Zhejiang University,
Hangzhou, 310027, Zhejiang, China

###### Abstract

Light propagation on a two-dimensional curved surface embedded in a three-
dimensional space has attracted increasing attention as an analog model of
four-dimensional curved spacetime in laboratory. Despite recent developments
in modern cosmology on the dynamics and evolution of the universe,
investigation of nonlinear dynamics of light in non-Euclidean geometry is
still scarce and remains challenging. Here, we study classical and wave
chaotic dynamics on a family of surfaces of revolution by considering its
equivalent conformally transformed flat billiard, with nonuniform distribution
of refractive index. This equivalence is established by showing how these two
systems have the same equations and the same dynamics. By exploring the
Poincaré surface of section, the Lyapunov exponent and the statistics of
eigenmodes and eigenfrequency spectrum in the transformed inhomogeneous table
billiard, we find that the degree of chaos is fully controlled by a single
geometric parameter of the curved surface. A simple interpretation of our
findings in transformed billiards, the ”fictitious force”, allows to extend
our prediction to other class of curved surfaces. This powerful analogy
between two a prior unrelated systems not only brings forward a novel approach
to control the degree of chaos, but also provides potentialities for further
studies and applications in various fields, such as billiards design, optical
fibers, or laser microcavities.

##  I Introduction

Einstein’s general theory of relativity (GR) unprecedentedly interprets
gravity in a geometrical framework, namely, the presence of a massive object
distorts the very fabric of the space and time in its vicinity. Despite of
GR’s triumph, gravitational effects are in principle too feeble for light to
be perceived in laboratory environment. One of the analog models of GR is two-
dimensional (2D) curved surfaces embedded in 3D space, whose theoretical
cornerstones are the 3+1 membrane paradigm Thorne1986  and embedding diagrams.
After taking a constant time and extracting the equatorial slice of a
spherically symmetrical space, the remnant metric can be visualized through a
2D curved surface. Ever since this innovative notion was brought up by Batz
and Peschel Batz2008 , electromagnetic (EM) dynamics and wave propagation on
2D curved surfaces have been prosperously developed both theoretically
Batz2010 ; Bekenstein2014 ; Lustig2017 ; XuPRA ; XuOE ; XuNJP ; Wang2018 ;
Jing2021  and experimentally Schultheiss2010 ; Schultheiss2015 ; Patsyk2018 ;
Bekenstein2017 ; Li2018 . On the other hand, GR notions could reciprocally
contribute to the engineering of nanophotonic devices Bekenstein2017 ;
Lebental2021  and transformation optical designs, ranging from EM invisible
cloaks Tyc2009 ; Zolla2009  to perfectly focusing lens LXu2019 ; LXu2019OE ;
Patrick2015 , based on manipulation of the curved geodesic paths.

Investigation of nonlinear dynamics in the context of GR may date back to
1960s, and has been extensively discussed in various aspects Radu2017 ;
Contopoulos2008 ; Aydiner2018 . Nonlinear dynamics, or chaos, widely exists in
diverse modern scientific disciplines. One characteristic of chaos is the
exponential sensitivity to initial conditions, with __ a well-publicized
metaphor being the “butterfly effect”. In the recent decades, instead of being
considered a nuisance, wave chaos has been taken advantage of to ameliorate
optical resonators and laser microcavities Liu2013 ; Cao2018 ; Cho1998 ;
Stone1997 ; Kim2004 ; Hentschel2008 ; Xiao2017 ; Rotenberg2015 ;
Fratalocchi2014 ; Fan2021 , such as enhancing energy storage Liu2013 ,
suppressing spatiotemporal lasing instabilities Cao2018 , realizing high power
directional emission Cho1998 ; Stone1997 ; Kim2004 , etc. To this end,
information about dynamical behaviors in 2D cavity, such as critical periodic
orbits and volume of chaotic area in phase space, is significant. One paradigm
model is a 2D table (or flat) billiard, which is typically a closed domain
wherein particles and light rays propagate freely, except for elastic
collisions on its boundaries. Various degrees of chaos in table billiards are
realized by either deforming the boundaries Sinai1963 ; Berry1981 ;
Bunimovich2001  or introducing external force Leonel2015 ; Berry1985 ;
Goodings1993 . In a recent work, a novel notion, transformed cavity, is
proposed with a nonuniform profile of refractive index present in a deformed
optical cavity Min2016 . Such transformed cavity opens up a new pathway of
chaos engineering.

Investigation of light chaotic dynamics in a non-Euclidean geometry remains
challenging from the complexity of the equations induced by the curved
geometry, as well as the ambiguity of concepts, such as the definition of
boundaries and indicators of chaos. To the best of our knowledge, the wave
chaos idea has been carried over in a few curved elastic systems
Sodergaard2010 ; Sodergaard2007 ; Rebinsky1994 , such as aluminum thin shells
Sodergaard2010 , yet has rarely been explored in optics, except for a very
recent study about ray chaotic behavior on a deformed toroidal surface
Chan2020 . Here we address this challenge by developing the analogy between a
curved surface and a 2D table billiard with nonuniform distribution of
refractive index, on the strength of transformation optics (TO). The
mathematical underpinning of TO is the form invariance of Maxwell equations
under general coordinate transformation, with the optical constitutive
parameters of the transformed media (i.e., engineering the permittivity
$\varepsilon$ and permeability $\mu$ tensor explicitly) Pendry2006 ;
Schurig2006 . As we will demonstrate, this approach of projection provides a
possibility to utilize physical systems on one surface to explore their
counterparts on the other, especially when the studies in one of the systems
are demanding to carry on.

In this work, we study classical and wave chaos on a special type of 2D curved
surfaces in 3D space, surfaces of revolution (SORs), by investigating light
propagation in its corresponding table billiard with azimuthally symmetric
nonuniform distribution of refractive index. We exemplarily choose a typical
family of SORs, the Tannery’s pears. Since such surfaces are integrable, we
consider half of the pear and its 2D nonuniform billiard where a circular hole
is pierced to introduce chaos. We first establish that these two systems share
the same dynamics by showing that a conformal coordinate transformation
preserves both geodesic equations and wave equations. Thanks to this
equivalence, we can assess the degree of chaos in one system (the curved
surface) by probing the nature of the trajectories in Poincaré surface of
section (SOS) and measuring the Lyapunov exponent in the other (the flat
billiard). We find that the degree of chaos is fully controlled by a single
geometrical parameter of the Tannery’s pear. This is also revealed in the
statistics of eigenmodes of the Tannery’s pears, as the wave equation is
invariant for both systems. At last, we show that the above results can be
generalized to arbitrary SORs, by coming up with a universal quantity in the
2D nonuniform billiards. This simple interpretation further demonstrates the
potentialities of our approach.

##  II Results

### Basic theory

As its name indicates, a SOR can be generated by revolving an arbitrary plane
curve (known as generatrix) around an axis of symmetry for a circle. The line
element of a typical SOR takes the form
$ds^{2}=g_{\alpha\beta}dx^{\alpha}dx^{\beta}=E(u)du^{2}+G(u)dv^{2}$, with $u$
along longitudinal direction and $v$, being the rotational angle, along
transverse direction, as shown in Fig. 1A. In curved spaces, light rays
propagate along geodesics, in analogy with straight lines in free space. On
the basis of geodesic equations, an arbitrary geodesic with given initial
conditions is depicted by (for mathematical details, see SI Appendix, section
1)

$dv=\eta\frac{\kappa\sqrt{E(u)}}{G(u)\sqrt{1-\frac{\kappa^{2}}{G(u)}}}du.$ (1)

Here $\eta=\text{sgn}\left[\left(du/ds\right)_{\text{initial}}\right]$ is
determined by the initial direction, and slant $\kappa$, defined as
$\left[G(u)dv/ds\right]_{\text{initial}}$, remains invariant as a consequence
of the conservation of angular momentum Wang2018 . When it comes to waves, the
time-harmonic scalar EM field follows the massless Klein-Gordon equation

$\frac{1}{\sqrt{E(u)G(u)}}\frac{\partial}{\partial
u}\left[\sqrt{\frac{G(u)}{E(u)}}\frac{\partial\Phi}{\partial
u}\right]+\frac{1}{G(u)}\frac{\partial^{2}\Phi}{\partial v^{2}}+k^{2}\Phi=0,$
(2)

with $k$ being the propagation constant note .

In this work, instead of directly investigating light propagation on the
surface, we employ an innovative method by projecting a SOR onto a plane with
distribution of refractive index, in light of conformal transformation optics
Leonhardt2006 ; Chen2015 . This notion is mathematically underlain by the
theorem that any 2D Riemannian manifold is conformally flat. Namely, the
metric of an arbitrary 2D curved surface $ds^{2}$ can be conformally related
to the metric of a flat plane $ds^{2}_{\text{f}}$ through $ds^{2}=\Lambda
ds^{2}_{\text{f}}$, where $\Lambda$ is a differentiable function inverno1992 .
Interestingly, the right-hand side, defined as $ds^{\prime 2}$, perfectly
describes an inhomogeneous planar dielectric medium with spatially varying
refractive index $n^{\prime}=\sqrt{\Lambda}$. For SORs, thanks to their
rotational symmetry, we can naturally suppose the projected plane is
azimuthally symmetric. Consequently, polar coordinate is chosen out of
convenience and the variation of refractive index rests exclusively on radial
component. Based on the premise that $ds^{2}=ds^{\prime 2}$, an equivalence
between these two systems can be established through the coordinate mapping

$r(u)=A\exp\left[\int^{u}\sqrt{\frac{E(u^{\prime})}{G(u^{\prime})}}du^{\prime}\right]\text{,}\
v=\varphi\text{,}$ (3)

with distribution of refractive index

$n^{\prime}(u)=\frac{\sqrt{G(u)}}{A}\exp\left[-\int^{u}\sqrt{\frac{E(u^{\prime})}{G(u^{\prime})}}du^{\prime}\right].$
(4)

Here $A$ is a positive integration constant, $r$ and $\varphi$ are radial and
azimuthal coordinates, respectively.

By virtue of this coordinate mapping, transformed version of Eq. (1) appears
as

$d\varphi=\eta^{\prime}\frac{\kappa^{\prime}}{n^{\prime}(r)r^{2}\sqrt{1-\frac{\kappa^{\prime
2}}{n^{\prime 2}(r)r^{2}}}}dr,$ (5)

with
$\eta^{\prime}=\text{sgn}\left[\left(dr/ds^{\prime}\right)_{\text{initial}}\right]$
and $\kappa^{\prime}=\left[n^{\prime
2}(r)r^{2}d\varphi/ds^{\prime}\right]_{\text{initial}}$. We prove in SI
Appendix 2 that Eq. (5) is the very solution of geodesic equations on the
projected plane. Thus the dynamics of light rays are preserved in both
systems. In terms of waves, Eq. (2) is transformed into

$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial\Phi}{\partial
r}\right)+\frac{1}{r^{2}}\frac{\partial^{2}\Phi}{\partial\varphi^{2}}+k^{2}n^{\prime
2}(r)\Phi=0,$ (6)

which is exactly the wave equation of electric field on a plane with
distribution of azimuthally symmetric refractive index $n^{\prime}(r)$. As a
result, the transformation preserves the solutions of wave equation in one
system to be the eigenmodes of another. In SI Appendix 2, we extend this
result to arbitrary curved surfaces, demonstrating the universality of the
transformability of both geodesic equations and wave equations. We have shown,
therefore, that the dynamical properties of both optical rays and waves are
absolutely identical on a homogeneous SOR and on its transformed inhomogeneous
plane, and we can safely investigate one system to infer about the other.

When a SOR is closed (e.g. sphere, spindle) or infinitely extends (e.g.
cylinder, cone), its projection covers the whole plane. A paradigmatic example
is the so-called “Tannery’s pear”, which is an object of interest for
mathematicians Besse1978  and is going to be the study case of this paper. The
family of Tannery’s pears can be parametrized by $E(u)=(c+\cos u)^{2}$,
$G(u)=\sin^{2}u$, with $c$ being a positive constant and in the meantime, the
single parameter to describe this family of surfaces. It is proved that such
surfaces have a constant period function $2c\pi$ Borzellino2006 . That is,
when the parameter $c$ is an integer, all the geodesics (except the ones along
longitudinal $u$ direction) are closed, and for a particle which departs from
an arbitrary position on the geodesic, moves along it and returns to the
starting point, the accumulated variation in coordinate $v$ is $2c\pi$, as
shown in Fig. 1B. In order to limit our study to a finite area on plane, a
mirror is put on the “equator” (i.e., the latitude with $u=\pi/2$, see the
black dashed lines in Figs. 1C and 1D) of the surface, which plays the role of
a perfectly reflecting boundary, and only the lower half of surface is taken.
On the projected plane, the equator corresponds to a circle with unit radius,
while the lower half of surface corresponds to the area inside the circle,
which together naturally form a circular billiard, as shown in the right
column of Figs. 1C and 1D. These two systems are absolutely integrable, since
the number of constants of motion (energy and orbital angular momentum) equals
to the dimensionality of the system Gutzwiller1990 , owing to the rotational
symmetry of SORs. To introduce chaos in this system, an off-centered disk-
shaped area is eliminated, leaving a circular hole with a specular boundary in
the billiard. The choice of a disk shape is out of consideration of
simpleness, yet more complex shapes could be explored. The right column of
Figs. 1C and 1D shows two typical holes whose centers are located at $0.2$ (in
unit of radius) from the centers of the billiards, and whose radius are $0.3$
and $0.1$, denoted by dark green and orange dotted lines, respectively. By
taking the coordinate transformation Eq. (3), the correspondence of these two
holes on the surface are shown in the left column. It is seen that the
refractive index approaches to infinity at the center of the projected
billiard, corresponding to the apex of the surface in the bottom, which is
basically a singularity for light rays. We focus in the rest of the paper on
the cases where the hole includes the billiard center, resulting in the
removal of the bottom apex. Its correspondence on the surface is a boundary
unparallel to latitudes, rather than a hole.

In what follows, in order to distinguish them from their 2D projected
billiards, SORs are referred to as “3D surfaces”, based on the fact that they
are embedded in 3D space in spite of being 2D per se.

Figure 1: (A) Sketch of a typical SOR with the orthogonal curvilinear
coordinates on it. The inset illustrates its generatrix. (B) Tannery’s pear
with $c=2$ (left) and its projected billiard (right). A typical closed
trajectory is depicted by yellow solid line. (C),(D) Truncated Tannery’s pears
which have a circular hole with (C) radius $0.3$ (dark green dotted line) and
(D) radius $0.1$ (orange dotted line) in their projected billiards. The
equators are denoted by black dashed lines in the projected billiards. Two
typical trajectories are plotted by green and red solid lines, respectively.

### Poincaré surface of section

Figure 2: Poincaré surfaces of section of Tannery’s pears with different
parameter $c$. (A1) Sketch of Birkhoff coordinates. (B1),(C1),(D1) Sketches of
Tannery’s pears with $c=1$, $2$ and $5$, respectively. The transparent parts
correspond to areas outside the projected billiard and inside the hole, which
are excluded in the study. (A2),(B2),(C2),(D2) Poincaré surface of section of
Tannery’s pears with $c=0$, $1$, $2$ and $5$, respectively. Each color
represents an arbitrary trajectory. The diamonds and pentagrams respectively
indicate typical stable and unstable periodic orbits, whose representations in
Poincaré surface of section are a succession of points. These periodic orbits
in real space are exhibited in (A3), (B3), (C3), (D3), respectively,
corresponding to the points with the same label and color.

Ray trajectories become rapidly intractable in real space after few bounces,
it is therefore more advisable to inspect ray motions in phase space which are
composed of two spatial dimensions and their conjugate momenta. For billiard
systems, an effective approach to extract information from 4D phase space is
to record the states of the trajectory exclusively when it collides on the
outer boundary. This 2D section of phase space is known as Poincaré surface of
section (SOS), and is conveniently depicted in the framework of Birkhoff
coordinate, as sketched in Fig. 2A1. In practice, the trajectories between two
bounces are traced by Eq. (1). When the trajectories collide on the hole/lower
unparallel boundary, the rule of specular reflection is applied in the 2D
projected billiard, after projecting the trajectories on 3D surface back to
the inhomogeneous plane. This step saves one from complicated or even
impractible calculation on 3D surface, since the analytical expression of the
lower unparallel boundary is inaccessible.

Figure 3: (A) Variation of Lyapunov exponents versus the parameter $c$ of
Tannery’s pears and three different means of truncation. The distribution of
refractive index are illustrated in the inset with corresponding labels, and
the black cross denotes the center of the billiard. The red square data
markers in line (ii) correspond to the cases illustrated in Fig. 2. (B) The
quantity $\lambda_{N}$ at $N$th bounce, calculated by Eq. (7). Five pairs of
extremely close trajectories are randomly chosen at the case denoted by the
pink square data marker in (A), whose initial conditions are shown in the
legend with $\theta_{2}=\theta_{1}+10^{-4}$,
$p_{\theta_{2}}=p_{\theta_{1}}+10^{-4}$. All the five pairs of trajectories
converge to the same asymptotic value, which is defined as Lyapunov exponent.

When $\left|\sin\chi\right|$ is large enough, trajectories propagate near the
outer boundary of 2D projected billiard/equator of 3D surface and are
therefore free from colliding on the hole/lower unparallel boundary. Their
tangential momenta $\left|\sin\chi\right|$ remain invariant, so that in
Poincaré SOS these trajectories are represented by horizontal straight lines.
In contrast, trajectories with small $\left|\sin\chi\right|$ collide on the
hole/lower unparallel boundary, resulting in their representations ergodically
and irregularly distributed in the central part of Poincaré SOS, which
therefore forms a chaotic sea. Among these trajectories, some return to their
initial conditions after several bounces, and retrace themselves repeatedly.
Such periodic orbits (POs) are represented as a succession of points in SOS,
with the number of points equal to the number of bounces on the outer
boundary/equator in each period. For some POs, the trajectories whose initial
conditions are slightly deviated from the critical conditions of POs are able
to survive in the vicinity of the latter. Such POs, considered as stable POs,
form regular islands among the chaotic sea in Poincaré SOS together with their
quasi-periodic orbits. As a contrast, the other unstable POs are submerged in
the chaotic sea and can be scarcely recognized. Some typical stable and
unstable POs are illustrated by diamonds and pentagrams respectively in the
middle row of Fig. 2, and the corresponding trajectories in the billiards are
depicted in the last row with same labels and colors. An interesting feature
of PO with label $1$ is explicated in SI index, section 3.

A generic presentation about the Poincaré SOSs of four typical Tannery’s pears
with $c=0$, $1$, $2$, $5$ is given in Fig. 2. Different shapes of the surfaces
are clearly exhibited in Figs. 2B1, 2C1 and 2D1. Specially, when $c=0$, the
metric of the surface degenerates to $ds^{2}=\cos^{2}udu^{2}+\sin^{2}udv^{2}$,
which is equivalent to a polar coordinate if one takes $\sin u$ as the radial
component and $v$ as the azimuthal component. Thus, surface with $c=0$
corresponds to the flat billiard and serves as reference. The Poincaré SOSs
reveal that the phase spaces of all the four surfaces are mixed, with
coexistence of both regular and chaotic trajectories. One can observe that the
area of chaotic sea remarkably enlarges with the parameter $c$, implying an
increasingly chaotic phase space structure of Tannery’s pear with greater
parameter $c$. Another phenomenon one may meanwhile notice is an obvious
reduce in the amount of stable POs and the area of islands, especially in
Figs. 2B2, 2C2 and 2D2.__ These two hints suggest an increasing proportion of
trajectories transferring from regular to chaotic, and consequently indicate a
more significant degree of chaos.__ Note that __ this observation is true both
on 3D surface and in 2D billiard with refractive index landscape, because of
their equivalence.__

### Lyapunov exponent

As was mentioned, a hallmark of chaos is the sensitivity to initial
conditions, i.e., the exponential divergence of two extremely nearby
trajectories in phase space. The speed of this exponential divergence is
characterized by the (maximum) Lyapunov exponent $\lambda$. Here, we adopt the
method of “billiard map” in terms of the collisions of trajectories on the
outer boundary, instead of the “billiard flow” with continuous time
Datseris2019 . In this way, Lyapunov exponent is defined as

$\lambda=\lim_{N\rightarrow\infty}\frac{1}{N}\sum\limits_{i=1}^{N}\ln\frac{\left|\boldsymbol{\delta}_{i}\right|}{\left|\boldsymbol{\delta}_{0}\right|},$
(7)

where $\left|\boldsymbol{\delta}_{i}\right|$ is the distance of the two chosen
trajectories in phase space at the $i$th bounce. In practice, two sets of
extremely close initial conditions are randomly chosen from the chaotic area
in Poincaré SOS, and distance $\left|\boldsymbol{\delta}\right|$ is collected
everytime two trajectories collide on the outer boundary. Theoretically, if
the number of bounces $N$ is large enough, Eq. (7) is supposed to approach to
a constant which is independent of the choice of initial conditions. Technical
details about the calculation of Lyapunov exponents are specified in SI index,
section 4.

By this method, Lyapunov exponents of a series of Tannery’s pears with varying
parameter $c$ and under three different means of truncation are calculated and
illustrated in Fig. 3A. Each data point in Fig. 3A is the result of averaging
over 4-6 randomly-chosen pairs of trajectories, where each pair of
trajectories experiences several tens of thousands of bounces on the outer
boundary/equator, and an approximate convergence has been reached (Example of
five typical pairs of trajectories is shown in Fig. 3B). The figure manifests
an overall increasing tendency of Lyapunov exponents with increasing parameter
$c$ of Tannery’s pears, indicating that in general, when parameter $c$
increases, the speed of divergence between two close trajectories, or the
instability of the trajectories increases, which validates the increase of
degree of chaos both on the 3D surface and in the 2D projected billiard.
Furthermore, the invariance of this increasing tendency in spite of different
positions and radii of the pierced hole in the projected billiard attests the
universality of this feature, as shown in Fig. 3A. Compared with the
qualitative interpretation obtained from Poincaré SOS with a coarse sampling
of parameter $c$, this signature brings a more quantitative insight, which
enables a finer exploration and hence reveals more details about the
dependence of chaoticity on parameter $c$. One example is that an exception to
this increasing trend is observed near $c=1$, elucidating that the variation
of chaoticity with parameter $c$ is not strictly monotonous.

### Statistics of eigenmodes and eigenfrequency spectrum

Figure 4: (A)-(H) Spatial statistics of eigenmodes. (A),(B) The field
distributions (amplitude) of a typical ergodic eigenmode (frequency
$\omega=2.9639$ GHz) of truncated Tannery’s pear with $c=1$ in projected
billiard (A) and on 3D surface (B). Its distribution of amplitude and
intensity are plotted in (C) and (D), respectively, with red solid line being
Gaussian and Porter-Thomas distribution as reference. The field distributions
(amplitude) of a typical scarring eigenmode (frequency $\omega=2.8863$ GHz)
and their statistics are shown in (E)-(H). The white solid lines in (E) and
(F) indicate its corresponding unstable periodic orbit. (J),(K) Statistics of
eigenfrequency spectrum. (J) Cumulative eigenfrequency density $N(E_{i})$ of
Tannery’s pear with $c=2$ versus eigenfrequency $\left\\{E_{i}\right\\}$. (K)
Distribution of nearest neighbor spacing of Tannery’s pears with $c=0$
(purple), $1$ (blue), $2$ (red) and $5$ (green). The Poisson (black dashed)
and Wigner-Dyson (black short dashed) distribution are also plotted as
reference. (Inset) Zoom-in of the area marked by red dotted rectangle.

When it comes to quantum or wave realm, notions of phase space and
trajectories are not properly defined due to the uncertainty principle.
Therefore, signatures in classical chaos are no longer valid, but are replaced
by quantum mechanical criteria based on energy spectrum, energy eigenvectors,
entanglement, temporal evolution of expectation values, etc. In this section,
we are going to focus on the first two fingerprints.

We have demonstrated that the wave equation on a SOR is identical with that in
its projected billiard, bridged by the coordinate transformation. Namely, the
performed calculation as well as the concepts and phenomena in one system can
be automatically extended to the other. Here, we obtain eigenmodes in both
systems with the help of COMSOL Multiphysics 5.2. In practice, the simulation
is conducted in the 2D inhomogeneous billiard, because of the difficulty in
constructing 3D surface (especially when the expression of the lower
unparallel boundary is unknown), and the eigenfunction on 3D surface is
subsequently obtained by projection from 2D billiard via the coordinate
transformation. Three typical eigenmodes in both systems are illustrated in SI
index, Fig. S4.

One credible imprint of ray chaos resides in the spatial statistics of
eigenmodes Legrand2013 . To be specific, the wave function of a typical
ergodic eigenmode, whose classical correspondence has a stochastic motion,
distributes uniformly over the whole available area of the phase space which
is ergodically visited by its classical trajectories. Such modes are thus
conjectured as a random superposition of plane waves with different
amplitudes, phases and directions but same wavenumber, manifesting an
analogous pattern of laser speckles. As per central limit theorem, such field
is Gaussian random, implying the amplitude of eigenmode follows the Gaussian
statistical distribution

$P\left(\Phi\right)=\frac{1}{\sqrt{2\pi\left\langle\Phi\right\rangle}}\exp\left(-\frac{\Phi^{2}}{2\left\langle\Phi\right\rangle}\right),$
(8)

while the probability distribution of intensity $I=\Phi^{2}$ is subjected to

$P\left(I\right)=\frac{1}{\sqrt{2\pi I/\left\langle
I\right\rangle}}\exp\left(-\frac{I}{2\left\langle I\right\rangle}\right),$ (9)

known as Porter-Thomas (PT) distribution. Figs. 4A and 4B illustrate a typical
ergodic eigenmode of truncated Tannery’s pear with $c=1$, both on 3D surface
and 2D projected plane. The good agreement of its $P\left(\Phi\right)$ and
$P\left(I\right)$ with Gaussian and PT distribution respectively validates the
ergodicity of the eigenmode.

Besides ergodic modes with speckle statistics, there exists a special class of
modes with enhanced amplitude in the vicinity of single unstable periodic
orbits in the corresponding classical system. This ubiquitous yet remarkable
phenomenon is well known as Quantum Scarring Heller1984 . Unlike the
enhancement of eigenstate intensity near stable periodic orbits which are well
understood by semiclassical theory of integrable systems, scar was initially
an unexpected phenomenon and later found to be a significant correction to
predictions from random matrix theory and Gutzwiller trace formula. One of
scarred modes of truncated Tannery’s pear with $c=1$ on projected plane and 3D
surface is exhibited in Figs. 4E and 4F, superimposed by a white solid line
indicating its corresponding unstable periodic orbit. A prominent deviation of
$P\left(\Phi\right)$ and $P\left(I\right)$ from ergodicity can be clearly
observed in Figs. 4G and 4H respectively.

Another classic and widely recognized indicator for the randomness of a
quantum system is the distribution of nearest neighbor spacing (NNS) of
eigenfrequency spectrum Zimmermann1986 . In a classically integrable system,
successive eigenfrequency distributes randomly, and the NNS $s$ follows the
Poisson law

$P_{\text{Poisson}}(s)=\exp(-s),$ (10)

with its peak located at zero. While in classically chaotic systems, the
presence of level repulsion leads NNS to best fit Wigner-Dyson (WD)
distribution

$P_{\text{WD}}(s)=\frac{\pi s}{2}\exp\left(-\frac{\pi}{4}s^{2}\right).$ (11)

Contrast with Poisson distribution where level spacing can be zero with
highest probability, a salient consequence of level repulsion is the
vanishment of infinitesimal NNS, i.e., $P(s)\rightarrow 0$ when $s\rightarrow
0$. This behavior can be well explained by the Bohigas-Giannoni-Schmit
conjecture Schmit1992  that spectra of time-reversal-invariant systems reveal
the same fluctuation properties as a Gaussian orthogonal ensemble in random
matrix theory. When the system has mixed dynamics, its distribution of NNS is
intermediate between the two limiting cases. In virtue of COMSOL Multiphysics
5.2, we got access to the eigenfrequency spectra and consequently performed an
unfolding procedure (For more details, see SI index, section 5). The resulting
distributions of NNS of truncated Tannery’s pears with $c=0,1,2,5$ are plotted
in Fig. 4K. In order to eliminate the effect of artificially-chosen interval
on the distribution curve, we adopt the cumulative distribution instead of
widely-used histogrammic representation. One can clearly observe that
distributions of NNS of all the four cases interpolate between Poisson and WD
distributions, indicating that these surfaces are mixing systems. Comparing
the distributions of the four cases, especially in the range of small spacing
$s$ which is zoomed in in the inset, the distribution curve gradually deviates
from Poisson distribution and inclines to WD distribution with the increase of
parameter $c$, manifesting a tendency of increasing chaoticity.

Figure 5: “Fictitious force” in projected billiards of Tannery’s pears with
different $c$. (A) Comparison of trajectories with same initial conditions,
$\theta=0$ and $p_{\theta}=0.8$, in projected billiards of Tannery’s pears
with $c=0$ (green), $1$ (red), $2$ (blue) and $5$ (orange). The star marks the
position of the center of billiard. (b) Comparison of the four surfaces on
$\left|\boldsymbol{\mu}(r)\right|$, the quantity defined to characterize the
degree of trajectories’ deflection.

##  III Discussion

The above-mentioned three signatures explicitly demonstrate an increasing
tendency of chaoticity with the geometrical parameter $c$. In this section we
propose a qualitative interpretation on this phenomenon. One should note that
it is the unparallel lower boundary on the 3D surface/hole in the projected
billiard, or namely, the breaking of symmetry that gives birth to the chaotic
dynamics in the first place Bunimovich2018 . On this basis, the presence of
curved surface further alters the chaoticity. There are two parallel
approaches to conduct the analysis, either on 3D surfaces or in 2D projected
billards. When considering from the perspective of 3D surfaces, on different
SORs trajectories follow their respective geodesics, the undisturbed natural
paths. Yet the lower unparallel boundaries, which are the reciprocal
projection of the hole in the projected billiards, have different shapes,
leading to different degrees of chaos. However, the various properties of
surfaces, such as the expression of generatrix, curvature, etc, could perplex
the discussion. Whereas, after being projected onto a plane, all the
information relevant to chaotic dynamics is completely and merely embodied in
the distribution of reflective index (since the location and radius of the
hole are fixed in all projected billiards with varying $c$), which simplifies
the problem to a great extent. Fig. 5A compares trajectories with identical
initial conditions in the projected billiards of Tannery’s pears with
parameter $c=0,1,2,5$. One can observe visually that with the presence of
nonuniform refractive index, the trajectories deflect towards the center of
the billiard, as if they were subjected to a “fictitious force”. More
importantly, in the projected billiard with greater parameter $c$, the
trajectories experience more conspicuous deflections, enhancing the chance of
collision on the hole. Put in another way, a greater initial tangential
momentum $p_{\theta}$ is required to realize the whispering gallery orbits,
resulting in a larger area of chaotic sea and consequently a greater
chaoticity of the system. We surmise that the degree of deflection is
proportional to the rate of variation of refractive index, and hypothetically
quantify the “fictitious force” to be

$\boldsymbol{\mu}=\frac{r}{n(r)}\boldsymbol{\nabla}n(r),$ (12)

with $r$ being the radial coordinate, and $n(r)$ being refractive index
profile. Details about construction of Eq. (12) are revealed in SI index,
section 6. Owing to the azimuthal symmetry of SORs, $\boldsymbol{\mu}$ is
along the radial direction. Further calculation shows that the magnitude of
$\boldsymbol{\mu}$ is determined uniquely by the metric of the surface (see SI
index, section 7). Besides, one may observe that the term in the absolute
value sign is always negative, indicating that in the projected billiards of
any SORs, the “fictitious force” always points towards the center of billiard.
Fig. 5B illustrates the increasing trend of $\left|\boldsymbol{\mu}(r)\right|$
with parameter $c$, which coincides with our assumption. Our interpretation is
supported by Fig. 3A, which shows the dependence of the Lyapunov exponent with
parameter $c$ for 3 different positions and diameters of the hole. The
difference between these three curves is a direct consequence of both the
pierced hole and presence of refractive index. Indeed, when the hole is moved
away from the center of billiard, or its size reduces, the billiard is farther
from being symmetric, leading to a higher degree of chaos. In the meantime,
areas with higher refractive index are free from being eliminated and are
revealed, further amplifying this difference.

We further inspect the general applicability of quantity
$\left|\boldsymbol{\mu}(r)\right|$ in SI index, section 7, by exploring some
other typical SORs and comparing their $\left|\boldsymbol{\mu}(r)\right|$s
with the volume of chaotic area in their Poincaré SOSs. The results confirm
the universality of $\left|\boldsymbol{\mu}(r)\right|$. When the studied
objects are generalized from Tannery’s pears to general SORs, different SORs
are parameterized by different parameters, and even for SORs which are from a
same family, the relation between their parameters and the degree of chaos
might not be as simple as it is in Tannery’s pears: that is where a universal
quantity could play a significant role. Another remarkable advantage of
defining $\left|\boldsymbol{\mu}(r)\right|$ is that given a SOR, one is able
to approximately estimate its degree of chaos by simply calculating its
$\left|\boldsymbol{\mu}(r)\right|$, which can be directly obtained by its
metric, instead of investigating more involved signatures of chaos, as we did
earlier.

##  IV Conclusion

In conclusion, we have investigated photonic chaotic dynamics (ray and wave
behaviors) on both 3D uniform SORs and in 2D table billiards with nonuniform
refractive index. By applying a conformal coordinate transformation, these two
systems are proved to be essentially equivalent for both light rays and waves.
This complete equivalence enables these two different systems to share
interesting phenomena and serve one another to solve in a simpler manner a
complex problem, since a geometry can appear easy and analytical in one
system, but untraceable in the other. The present proposal, projecting a 3D
SOR into a nonuniform table billiard, serves as an innovative pathway to help
both solve calculation and simulation problems and provide with a natural
interpretation of the role of curvature in 3D chaotic systems, which probably
could be done on the 3D SOR as well but are much simpler and more
straightforward to investigate in the latter. Compared with 3D surfaces, a
nonuniform 2D table billiard turns out to be a more promising candidate for
experimental design and realization, on account of the difficulty in
restricting and controlling lights on the former. Reciprocally, the landscape
of refractive index in nonuniform table billiards can be complex, while
transferring it to its corresponding 3D SOR could greatly reduce the
parameters. Our investigation also raises an interesting proposal on the
design of refractive index in a nonuniform billiard or cavity to achieve
expected features, taking advantage of special trajectories or geometrical
properties of its corresponding 3D surface. For instance, it is still an open
question for how to transfer an arbitrarily deformed chaotic billiard into
integrable, by correlating its refractive index profile with a SOR. More
underlying physics and potential applications about this system equivalence
remain to be explored.

Although we have taken a typical family of SORs, Tannery’s pears, as an
example, to study the variation of the degrees of chaos with respect to its
control parameter, the model can certainly be generalized to arbitrary SORs.
The parameter-dependent property of chaotic dynamics on 3D SORs enlightens an
efficient and neat mechanism to control and utilize chaos, and consequently
opens up many perspectives, for example, exploration and control of wave chaos
in multimode fiber amplifiers to enhance pump absorption efficiency
Legrand2003 ; Mortensen2007  with nonuniform transverse, or the extensive
study of a new type of microcavity Cho1998 ; Stone1997 ; Kim2004 ;
Hentschel2008 . Furthermore, exploration of chaotic dynamics within the
context of non-Euclidean geometry may open a new perspective in studies of
chaos in cosmology Bruni2019 ; Motter2003 ; Shellard1998 .

## Acknowledgments

The authors thank Tsampikos Kottos from Wesleyan University for fruitful
discussion on calculation of Lyapunov exponent, and Kun Tang from Bar-Ilan
University for help in COMSOL simulation. C. X. acknowledges the 2019 Israeli
'Sandwich Doctorate Program' for international students funded by the Council
for Higher Education at Bar-Ilan University. P. S. is thankful to the CNRS
support under grant PICS-ALAMO. This research is also supported by The Israel
Science Foundation (Grants No. 1871/15, 2074/15 and 2630/20), the United
States-Israel Binational Science Foundation NSF/BSF (Grant No. 2015694),
Zhejiang Provincial Natural Science Foundation of China under Grant No.
LD18A040001, the National Natural Science Foundation of China (NSFC) (grant
No. 11674284 and 11974309), and National Key Research and Development Program
of China (No. 2017YFA0304202).

## References

  * (1) K. S. Thorne, R. H. Price, D. A. MacDonald, Black Holes: The Membrane Paradigm (Yale University Press, 1986).
  * (2) S. Batz, U. Peschel, Linear and nonlinear optics in curved space. Phys. Rev. A 78, 043821 (2008).
  * (3) S. Batz, U. Peschel, Solitons in curved space of constant curvature. Phys. Rev. A 81, 053806 (2010).
  * (4) R. Bekenstein, J. Nemirovsky, I. Kaminer, M. Segev, Shape-preserving accelerating electromagnetic wave packets in curved space. Phys. Rev. X 4, 011038 (2014).
  * (5) E. Lustig et al., Curved-space topological phases in photonic lattices. Phys. Rev. A 96, 041804(R) (2017).
  * (6) C. Xu, A. Abbas, L. -G. Wang, S. -Y. Zhu, M. S. Zubairy, Wolf effect of partially coherent light fields in two-dimensional curved space. Phys. Rev. A 97, 063827 (2018).
  * (7) C. Xu, A. Abbas, L. -G. Wang, Generalization of Wolf effect of light on arbitrary two-dimensional surface of revolution. Opt. Express 26, 33263-33277 (2018).
  * (8) C. Xu, L. -G. Wang, Gouy and spatial-curvature-induced phase shifts of light in two-dimensional curved space. New J. Phys. 21, 113013 (2019).
  * (9) D. Wang, C. Liu, H. Liu, J. Han, S. Zhang, Wave dynamics on toroidal surface. Opt. Express 26, 17820-17829 (2018).
  * (10) H. W. Wu et al., Observation of protected localized states induced by curved space in acoustic topological insulators. arXiv [Preprint] (2021). arXiv: 2101.03294 (9 January 2021).
  * (11) V. H. Schultheiss et al., Optics in curved space. Phys. Rev. Lett. 105, 143901 (2010).
  * (12) V. H. Schultheiss, S. Batz, U. Peschel, Hanbury Brown and Twiss measurements in curved space. Nat. Photon. 10, 106-110 (2016).
  * (13) A. Patsyk, M. A. Bandres, R. Bekenstein, M. Segev, Observation of accelerating wave packets in curved space. Phys. Rev. X 8, 011001 (2018).
  * (14) J. Zhu, Y. Liu, Z. Liang, T. Chen, J. Li, Elastic waves in curved space: mimicking a wormhole. Phys. Rev. Lett. 121, 234301 (2018).
  * (15) R. Bekenstein et al., Control of light by curved space in nanophotonic structures. Nature Photon. 11, 664-670 (2017).
  * (16) Y. Song et al., Möbius strip microlasers: a testbed for non-Euclidean photonics. arXiv [Preprint] (2020). arXiv: 2011.12088 (24 November 2020).
  * (17) U. Leonhardt, T. Tyc, Broadband invisibility by non-Euclidean cloaking. Science 323, 110-112 (2009).
  * (18) A. Nicolet, F. Zolla, Cloaking with curved spaces. Science 323, 46-47 (2009).
  * (19) L. Xu et al., Light rays and waves on geodesic lenses. Photonics Res. 7, 1266-1272 (2019).
  * (20) L. Xu, T. Tyc, H. Chen, Conformal Optical devices based on geodesic lenses. Opt. Express 27, 28722-28733 (2019).
  * (21) G. Lefebvre et al., Experiments on Maxwell’s fish-eye dynamics in elastic plates. Appl. Phys. Lett. 106, 024101 (2015).
  * (22) P. V. P. Cunha, C. A. R. Herdeiro, E. Radu, Fundamental photon orbits: black hole shadows and spacetime instabilities. Phys. Rev. D 96, 024039 (2017).
  * (23) A. E. Motter, P. S. Letelier, FRW cosmologies between chaos and integrability. Phys. Rev. D 65, 068502 (2008).
  * (24) E. Aydiner, Chaotic universe model. Sci. Rep. 8, 721 (2018).
  * (25) C. Liu et al., Enhanced Energy storage in chaotic optical resonators. Nat. Photon. 7, 473-478 (2013).
  * (26) S. Bittner et al., Suppressing spatiotemporal lasing instabilities with wave-chaotic microcavities. Science 361, 1225-1231 (2018).
  * (27) C. Gmachl et al., High-power directional emission from microlasers with chaotic resonators. Science 280, 1556-1564 (1998).
  * (28) J. U. Nöckel, A. D. Stone, Ray and wave chaos in asymmetric resonant optical cavities. Nature 385, 45-47 (1997).
  * (29) M. S. Kurdoglyan, S. -Y. Lee, S. Rim, C. -M. Kim, Unidirectional lasing from a microcavity with a rounded isosceles triangle shape. Opt. Lett. 29, 2758-2760 (2004).
  * (30) J. Wiersig, M. Hentschel, Combining directional light output and ultralow loss in deformed microdisks. Phys. Rev. Lett. 100, 033901 (2008).
  * (31) X. Jiang et al., Chaos-assisted broadband momentum transformation in optical microresonators. Science 358, 344-347 (2017).
  * (32) C. Liu et al., Triggering extreme events at the nanoscale in photonic seas. Nat. Phys. 11, 358-363 (2015).
  * (33) C. Liu, A. Di Falco, A. Fratalocchi, Dicke phase transition with multiple superradiant states in quantum chaotic resonators. Phys. Rev. X 4, 021048 (2014).
  * (34) L. Fan, X. Yan, H. Wang, L. V. Wang, Real-time observation and control of optical chaos. Sci. Adv. 7, eabc8448 (2021).
  * (35) Ya. G. Sinai, On the foundations of the ergodic hypothesis for a dynamical system of statistical mechanics. Sov. Math Dokl. 4, 1818-1822 (1963).
  * (36) M. V. Berry, Regularity and chaos in classical mechanics, illustrated by three deformations of a circular 'billiard'. Eur. J. Phys. 2, 91-102 (1981).
  * (37) L. A. Bunimovich, Mushrooms and other billiards with divided phase space. Chaos 11, 802-808 (2001).
  * (38) D. R. da Costa, C. P. Dettmann, E. D. Leonel, Circular, elliptic and oval billiards in a gravitational field. Commun. Nonlinear Sci. Numer. Simulat. 22, 731-746 (2015).
  * (39) M. Robnik, M. V. Berry, Classical billiards in magnetic fields. J. Phys. A: Math. Gen. 18, 1361-1378 (1985).
  * (40) T. Szeredi, D. A. Goodings, Classical and quantum chaos of the wedge billiard. I. classical mechanics. Phys. Rev. E 48, 3518-3528 (1993).
  * (41) Y. Kim et al., Designing whispering gallery modes via transformation optics. Nat. Photonics 10, 647-653 (2016).
  * (42) M. Avlund, C. Ellegaard, M. Oxborrow, T. Guhr, N. Søndergaard, Observation of periodic orbits on curved two-dimensional geometries. Phys. Rev. Lett. 104, 164101 (2010).
  * (43) G. Tanner, N. Sødergaard, Wave chaos in acoustics and elasticity. J. Phys. A: Math. Theor. 40, R443-R509 (2007).
  * (44) A. N. Norri, D. A. Rebinsky, Membrane and flexural waves on thin shells. J. Vib. Acoust. 116, 457-467 (1994).
  * (45) D. Wang, C. Liu, S. Zhang, C. T. Chan, Chaotic photon spheres in non-Euclidean billiard. Nanophotonics 9, 3367-3372 (2020).
  * (46) J. B. Pendry, D. Schurig, D. R. Smith, Controlling electromagnetic fields. Science 312 (5781), 1780-1782 (2006).
  * (47) D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, D. R. Smith, Metamaterial electromagnetic cloak at microwave frequencies. Science 314 (5801), 977-980 (2006).
  * (48) The precise evolution of vector light fields is described in Section II of Ref. Batz2008 , starting from 3D Maxwell equations written in general coordinates. Specifically, for EM wave polarizing strictly along transverse direction, its wave equation also takes the form of Eq. (2), except for an extra curvature-induced potential, whose effect is discussed in Ref. XuNJP .
  * (49) U. Leonhardt, Optical conformal mapping. Science 312 (5781), 1777-1780 (2006).
  * (50) L. Xu, H. Chen, Conformal transformation optics. Nature Photon. 9, 15-23 (2015).
  * (51) R. D’Inverno, Introducing Einstein’s Relativity (New York, Oxford University Press, 1992).
  * (52) A. L. Besse, Manifolds all of whose geodesics are closed (Berlin-New York, Springer-Verlag, 1978).
  * (53) J. E. Borzellino, C. R. Jordan-Squire, G. C. Petrics, D. M. Sullivan, On the existence of infinitely many closed geodesics on orbifolds of revolution. arXiv [Preprint] (2006). arXiv: math/0602595 (27 Feburary 2006).
  * (54) M. C. Gutzwiller, Chaos in Classical and Quantum Mechanics (New York, Springer-Verlag, 1990).
  * (55) G. Datseris, L. Hupe, R. Fleischmann, Estmating Lyapunov exponents in billiards. Chaos 29, 093115 (2019).
  * (56) O. Legrand, F. Mortessagne, ”Wave chaos for the Helmholtz equation” in New Directions in Linear Acoustics and Vibration: Quantum Chaos, Random Matrix Theory and Complexity, M. Wright, R. Weaver (Cambridge, GB. Cambridge University Press, 2010).
  * (57) E. J. Heller, Bound-state eigenfunctions of classically chaotic Hamiltonian systems: scars of periodic orbits. Phys. Rev. Lett. 53, 1515-1518 (1984).
  * (58) Th. Zimmermann, H. -D. Meyer, H. Köppel, L. S. Cederbaum, Manifestation of classical chaos in the statistics of quantum energy levels. Phys. Rev. A 33, 4334-4341 (1986).
  * (59) E. B. Bogomolny, B. Georgeot, M. J. Giannoni, C. Schmit, Chaotic billiards generated by arithmetic groups. Phys. Rev. Lett. 69, 1477-1480 (1992).
  * (60) L. A. Bunimovich, Mechanisms of chaos in billiards: dispersing, defocusing and nothing else. Nonlinearity 31, R78-R92 (2018).
  * (61) L. Philippe et al., Experimental study of pump power absorption along rare-earth-doped double clad optical fibres. Opt. Commun. 218, 249-254 (2003).
  * (62) N. A. Mortensen, Air-clad fibers: pump absorption assisted by chaotic wave dynamics? Opt. Express 15, 8988-8996 (2007).
  * (63) C. Ganguly, M. Bruni, Quasi-isotropic cycles and nonsingular bounces in a mixmaster cosmology. Phys. Rev. Lett. 123, 201301 (2019).
  * (64) A. E. Motter, Relativistic chaos is coordinate invariant. Phys. Rev. Lett. 91, 231101 (2003).
  * (65) N. J. Cornish, E. P. S. Shellard, Chaos in quantum cosmology. Phys. Rev. Lett. 81, 3571-3574 (1998).