11institutetext: Warsaw University of Technology, Warsaw, Poland
11email<EMAIL_ADDRESS>
22institutetext: Systems Research Institute Polish Academy of Sciences
22email<EMAIL_ADDRESS>
# Applying machine learning to predict behavior of bus transport in Warsaw,
Poland

Łukasz Pałys 11 Maria Ganzha 11 0000-0001-7714-4844 Marcin Paprzycki 22
0000-0002-8069-2152

###### Abstract

Nowadays, it is possible to collect precise data describing movements of
public transport. Specifically, for each bus (or tram) geoposition data can be
regularly collected. This includes data for all buses in Warsaw, Poland.
Moreover, this data can be downloaded and analysed. In this context, one of
the simplest questions is: can a model be build to represent behavior of
busses, and predict their delays. This work provides initial results of our
attempt to answer this question.

##  1 Introduction

One of key aspects of development of Smart Cities is the growing number of
installed sensors. To make city really smart, streams of data generated by
city-sensors are made open to analysis. Moreover, sometimes, privately owned
sensors open data to contribute to the well-being of city inhabitants.

Following this trend, Warsaw, Poland, opened its data resources. Here,
available data includes, among others, culture, education and urban transport.
For instance, data about location of each city bus is regularly collected. The
question that arises is: how this data can be used to deliver service to
people using public transport? To answer this question we have applied machine
learning to predict delays in public transport. Note that such solution (1)
could improve effectiveness of movement by means of public transport.
Currently, it often happens that public transport vehicles (in particular,
buses) arrive inconsistently with their timetable. A delay prediction system
(2) would save time spent waiting for a late bus(es). It would also allow for
a (3) better planning of journeys requiring a change, taking into account the
risk of a late arrival at the bus change location. Another potential advantage
of the proposed solution is (4) possibility of predicting “typical delays”
and, on their basis, modify the timetables. Another way of using the results
could be to (5) modify traffic organization, e.g. by introducing bus-passes,
thus helping busses to smoothly pass through problematic locations. Moreover,
(6) from the results it may be possible to obtain information, which factors
influence the delays to the greatest extent. It may be also possible (7) to
establish which factors can be successfully captured in machine learning
models, and which factors would require special treatment and may be difficult
to incorporate. As a result, in the future, (8) it may be possible to
undertake actions aiming at removing, or reducing, some causes of delays.

Reported work is based on actual data provided by the city of Warsaw. This
data contains bus schedules and their geolocation. Note that tram information
is also available. However, we have decided to start our explorations from bus
lines. Since our results are promising, we plan to extend our investigations
further.

##  2 Related work

Predicting urban bus travel time has been discussed in the literature. Applied
methods can be divided into several groups. First, the statistical methods,
such as kNN [8, 2], regression [2], prediction based on sequence patterns [10]
or time series [11], and application of Kalman filter [6]. Second, methods
based on observations, which include the historical means approach [1, 5].
Next, methods based on machine learning, which include: back propagation
neural network [12], radial basis function network [13], multilayer perceptron
[5, 1]. The last group consists of hybrid methods that combine multiple
algorithms into a single model [4, 14]. Let us now discuss the selected
methods in more detail.

The Pattern Sequence-based Forecasting technique was used in [10], to predict
bus travel times in Chennai, India. Discrete Wavelet Transform, and two
clustering algorithms were used to extract travel time trends. Results showed
that predictions based on the k-means algorithm were more accurate.
Predictions have been tested on sections with low, medium and high travel time
variability. For sections with low and medium travel time variance, the
results were satisfactory. However, for sections with large time differences,
they were not.

The time series approach was used in [11] to forecast bus travel time in Lviv,
Ukraine. The authors used STL (Seasonal and Trend decomposition using Loess)
time series decomposition to identify patterns. Travel time data, from (5)
working days, from September-October 2018 was collected for a single bus line,
in both directions. The training set consisted of 9 weeks, while test data was
from week 10. Systematic variability of traffic intensity was noticed. For
travel towards the city center, the highest congestion occurred in the
morning, due to commuting to work. For travel from the city, the greatest
crowd was observed in the evening. The predictive model obtained a mean
absolute percentage error (MAPE) of less than 8%. The average error, of
predicting travel time towards city center was approximately 3 minutes, and 2
minutes in the opposite direction.

In [13], authors used Radial Basis Function Neural Network (RBFNN) to predict
bus travel times in Dalian, China. The model was trained on historical data,
with additional correction of results. The input data of the model was: time
of arrival at stop “i”, distance between stops “i” and “i + 1”, number of
passengers leaving/entering at stop “i”, delay, average speed between stops
“i” and “i + 1”. Correction was based on estimating travel time between
consecutive stops, using Kalman Filter. Obtained results were compared with
use of BPNN and RBFN networks, without data correction. Overall, the RBFNN
model, with online tuning, was the best, with the MAPE at 7.59%. In
comparison, the MAPE values of the non-tuned BPNN and RBFNN models were 17.41%
and 15.98%.

Bus travel time prediction, described in [14], is based on deep learning and
data fusion. Here, travel time was divided into driving time and stopping
time. Bus GPS data was applied. The models used fourteen features, including:
stop ID, day of the week, time, current speed (based on GPS), stopping time at
a stop, travel time between stops. Data used for testing came from Guangzhou
and Shenzhen (China) and represented a single line in each city. To avoid one-
hot encoding111https://en.wikipedia.org/wiki/One-hot, they used embedding was
used for for: time of day, day of the week, stop ID. Moreover, since the size
of the “discrete features” was much smaller than that of the continuous ones,
data was unbalanced. Hence, the Wide and Deep model was used. The Wide
component was responsible for the embedded discrete data, while the Deep
component analyzed the continuous data using the Deep Neural Network (DNN),
and results of both components were concatenated into a feature list. Next,
using the Recurrent Neural Network (RNN), the feature weights were determined.
Travel time was predicted using DNN. Prediction results were compared with
historical averages. Proposed model outperformed other approaches with MAPE at
8.43% (Gangzhou, line 261). Moreover, MAPE of predictions for the peak hours
was 4-7% lower than they were 9.51% (for morning hours) and 6.43% (for evening
hours). In [1], dynamic travel time prediction model was developed, to predict
approximate arrival times at subsequent stops. GPS data, from November, 2008
to May, 2009 from Macea (Brazil) for single bus line with 35 stops was used.
Here, Historical Average, Kalman Filtering and Artificial Neural Network were
tried. Authors discovered considerable variation in average travel time at
different times of the day. For the northern direction, in rush hours
(17:00-18:00), travel time was 30% longer than the average time. For the south
direction, the mean time during peak hours (7:00-7:30) was about 20% longer.
This fact was taken into account when developing the model. However, Ducan’s
multiple range test showed that there were no significant differences in
travel time distributions in both directions.

The 3-layer perceptron achieved the best result, with overall MAPE at 18.3%.
Importantly, the article shows that bus travel times can be reasonably
predicted using only arrival and departure times at stops.

In [2], the review of methods of predicting bus arrival times, found in [3],
was extended. Authors used: k-Nearest Neighbor, Artificial Neural Network, and
Super Vector Regression to predict bus travel times in Trondheim (Norway).
Used data included: the route of the vehicle, planned and actual time of
arrival, and stopping time. To check the influence of individual time
intervals, 423 different data sets were created. To evaluate model
performance, mean absolute error (MAE) was used. Depending on data set, (best)
MAE varied between 61 and 86 seconds. Moreover, each method delivered a
“winner” (best result for a data set). Separately, additional attributes,
e.g.: weather, football matches, tickets, were tried. However, this
information did not bring significant improvements. The ticket data improved
the results in some training periods of the ANN and KNN models.The weather
impact on foreseen methods also have been researched and authors discovered
that using travel time instead of delay time will provide more accurate and
stable results. The result of the study was the confirmation of this thesis.

Comparing the described classifiers, the authors found that the best model for
real-time prediction is kNN. Together with SVR, it achieved the best results,
it’s training time was by far the shortest and forecasts were the most stable.
The analysis of the data sets allowed the authors to conclude that the best
way to predict travel time on a given day is to train on the previous three
days of the same type. This means that if the goal is to predict the coming
Friday, the model should be trained on the last three fives.

In [4], an ANN was used with historical GPS data and an automatic toll
collection system data. Moreover, impact of intersections with traffic lights
was taken into account. Data came from Jinan (China), bus line 63, and was
collected for 30 days from November 1, 2010. In the model, travel time and
time distances between successive buses were used. The distribution of time
intervals between buses on weekdays and at weekends was analyzed. It was found
that these values fluctuate significantly during weekdays, while during
weekends the distribution was more stable. Moreover, the average travel time
and travel time variances for the following days of the week were determined.
Initially an ANN with single hidden layer was used. However, due to travel
time fluctuations, finally, a hybrid ANN (HANN) model was used. It consists of
separate subnets, trained in real-time on different subsets of data, to
predict of travel time within specific period, depending on the days (working
days, weekends) and peak hours. Comparison of prediction results using the
relative prediction error (RPE) index demonstrated that the ANN and HANN
models are more reliable than the Kalman Filter. Moreover, the HANN model is
better suited for short-distance prediction.

A comparison of the methods of predicting trams travel time in Warsaw, is
presented in [5]. Two data sources were used: (1) tram timetables, including
locations of stops and routes, and planned times of arrivals; and (2) the GPS-
based positions of trams. Authors tested: delay propagation, historical
average (HA), and an ANN model. It was found that, overall, the HA was the
best method. However, ANN was best in predicting delays in rush hours.

Overall, the main findings from the analysed literature can be summarized as
follows. (a) Reported results have been focused mostly on a single city. Only
in one case two different cities have been approached. In all cases a single
bus line was studied. Hence, there is no “benchmark” data that performance of
proposed approaches can be assessed against. (b) The main methods used to
predict bus travel time (delays) were: (i) statistical methods – k-nearest
neighbours (kNN), regression model, Kalman filter, (ii) historical observation
methods (HA), (iii) machine learning methods – back propagation neural
networks, radial basis function networks, multilayer perceptrons, and (iv)
hybrid methods combining the above algorithms into one model (HANN). Out of
these methods, the best results have been reported for HANN, HA, RBFN and
multilayer perceptron. However, it was also found that use of HANN requires
substantially larger set of training data. (c) Quality of predictions has been
measured using: (i) mean absolute error (MAE; in seconds); (ii) mean
percentage absolute error (MAPE); (iii) standard deviation (STD; in seconds).
(d) The simplest approaches used GPS data alone. Other popular data elements
were: information about sold tickets and about bus speed. Additionally,
effects of non-travel events (e.g. games or weather) have been
(unsuccessfully) tried. (e) Best accuracy was reported in [2], where MAE was
40s. However, predictions were for 1 to 16 stops only. A good accuracy for a
longer bus line was with MAE of the order of 70. However, for the same line,
the kNN model performed better, with MAE at 61s. On the other hand, for long
Warsaw tram lines (more than 40 stops), the MAE was at 148s. However, in this
case the HA reduced MAE to 123s.

Based on these findings, the approach and the main contributions of this work
can be summarized as follows. (A) The MLP, RBFN and HA models were selected.
(B) While still using data from a single city, bus lines, with different
characteristics (e.g. normal, express, cross-town, long rout, short rout,
etc.) have been studied. (C) Performance of different approaches to predicting
travel time at subsequent stops has been evaluated. (D) Influence of
additional features (data from timetables, number of intersections with
traffic lights, information whether the distance between stops includes a bus
lane, weather) on the prediction accuracy has been evaluated. (E) Hybrid
model, composed of models predicting travel time for various distances, has
been proposed.

##  3 Available data, its preprocessing, and experimental setup

Let us now describe the data set, and applied preprocesseing. As a part of the
project Open data in Warsaw222https://api.um.warszawa.pl/, data that includes
the exact location of public transport vehicles, reported in real time, is
available. From there, 30 days worth of data about busses was harvested. On
weekdays, approximately 210MB of data were downloaded. On Saturdays and
Sundays, the data was approximately 160 MB. Harvesting lasted from March 8,
2021 till May 10, collecting about 10 GB of data. Specifically, in order to
download the data, a script was written, which queried the API every 30
seconds, to retrieve locations of all buses in use. In response, JSON with a
list of busses and their exact locations was received. A single CSV record
contained: (1) line number, (2) vehicle ID, (3) brigade number, (4) time, (5)
geographic coordinates of the location: Based on analysis, presented in
Section 2, it has been decided to use inputs similar to these reported in the
literature. In this way, while direct comparison of results is impossible, an
indirect comparison will be attempted. Hence, preprocessing was applied,
resulting in records containing: (1) line number, (2) departure time from the
last stop, (3) current percentage of distance traveled between adjacent stops,
(4) time of the last GPS signal, (5) current time, (6) driving direction.
Here, use of time of arrival at the stop was based on [5].

Additionally, file containing timetable of buses, their routes, including list
and GPS coordinates of stops, and departure times from individual stops, was
downloaded from the ZTM website333https://www.wtp.waw.pl. A special script has
been prepared to extract all necessary data from this file. Pre-processing of
this file created two tables, containing: (1) geographical coordinates of all
stops, and (2) lists of all stops for a given line and direction. After
preprocessing, file with: line number, time, vehicle and brigade
number,driving direction, number of the next stop to visit on the route,
information whether the vehicle is at the stop, the percentage of distance
traveled between consecutive stops, was created (see, Figure 1).

Figure 1: Preprocessed data sample (A/B indicate different directions of the
same route

Separately, file with additional features has been prepared, containing:
section congestion, timetable, and weather information. The congestion
captures number of busses traveling at a given moment the route between
neighboring stops. While in [2] it was claimed that weather information did
not significantly improve the prediction quality, it was decided to test this
for Warsaw busses.

From available data, 29 bus lines, that sent an average of more than 22,000
GPS signals a day, were selected. Based on careful analysis of characteristics
of their routes, within the context of the city itself, selected bus lines
have been manually grouped into eight groups with similar characteristics:

  1. 1.

Long routes within periphery North-South: bypassing the City Center, running
on the western side of the Vistula River; lines: 136, 154, 167, 187, 189.

  2. 2.

Long routes within periphery West-East: bypassing the City Center, crossing
the Vistula River; lines 112, 186, 523.

  3. 3.

Centre-periphery: routes with one end in the City Center, running to the
peripheries; not crossing the river; lines 131, 503, 504, 517, 518.

  4. 4.

Long routes through Center (with ends on peripheries); lines 116, 180, 190.

  5. 5.

Express: a fairly straight long lines with small number of stops (typically
around 1/3 stops of “normal” lines); lines 158, 521, 182, 509.

  6. 6.

Centre-Praga: short routes starting in the City Center; crossing the Vistula
River; lines 111, 117, 102.

  7. 7.

Short lines within peripheries: located in Western part of Warsaw; lines 172,
191, 105.

  8. 8.

Short lines within the Center: not crossing the river; lines 128, 107, 106.

The base MLP and RBFN models described in [5] and [13] used vehicle location
data in time. We experimented with these approaches by adding:

  * •

current traffic volume between stops,

  * •

weather,

  * •

number of buses on a section between stops,

  * •

information whether the distance between stops includes a bus lane,

  * •

number of intersections with traffic lights,

  * •

schedule data.

Multiple architectures, accepting bus location and different subsets of
available additional data have been tested. For each model, the architecture
that gave best results for a given set of experiments, is reported. They are
compared to the architectures that use only bus location data.

Travel time distributions differs between days of the week (primarily between
working days and weekends). Therefore, following [2], models were trained on
data from the same day, from three consecutive weeks, and tested on data from
the last week. Finally, following work reported in Section2, in the
experiments, the accuracy was measured using MAE and STD.

##  4 Experimental results and their analysis

Let us now summarize experimental results. They are reported representing the
“process of discovery” that took place in our work.

Total travel time prediction Here, two methods have been tried: “recursive”
and “long distance”. In both: number of the next stop, travel times between
last three stops, current percentage of the distance traveled between stops,
number of the destination stop, current time, route direction, were used.
Applied models returned total travel time to the destination.

In the recursive method, the model is trained on a data set containing
information about travel time to the nearest stop (i.e. focus on prediction of
travel time between neighboring stops). Estimating travel time from stop $n$
to stop $n+k$ consists of predicting $k$-steps using the trained model.
Specifically, the result from the previous stop is included in the input data
for the next prediction.

The long distance method is based on prediction of travel time to a specific
stop. Here, the model is trained on data containing information about total
pertinent travel time. The assumption was that this approach would be worse
for short distances, but may be better at predicting travel time of long(er)
trips.

The comparison was made for the bus line 523. The training set consisted of
data from Thursdays: March 11, 18, 25. The testing data was from Thursday,
April 1. The number of records in the training set was over 1 million, and in
the test set over 300,000. Because the goal of this experiment was to compare
the methods, not the models, both approaches used the same neural network
architecture, i.e. MLP with two hidden layers consisting of 6 and 24 neurons,
with ReLU activation function. results are presented in Tables 1 and 2.

Table 1: Prediction of total travel time using the recursive method (bus line
523)

distance time | 7:00 - 10:00 | 10:00 - 14:00 | 14:00 - 19:00 | 19:00 - 23:00
---|---|---|---|---
MAE | STD | MAE | STD | MAE | STD | MAE | STD
1 - 3 | 58.24s | 62.76s | 60.01s | 79.94s | 59.51s | 75.74s | 47.02s | 54.01s
4 - 8 | 125.52s | 129.44s | 130.89s | 141.44s | 131.74s | 124.96s | 92.10s | 88.65s
9 - 15 | 211.36s | 165.19s | 214.10s | 179.41s | 240.66s | 170.14s | 156.74s | 121.10s
16 - 20 | 294.83s | 188.44s | 310.42s | 219.57s | 344.49s | 201.74s | 222.66s | 137.57s
21 - 27 | 360.19s | 200.56s | 398.48s | 266.61s | 407.73s | 217.13s | 277.52s | 154.40s

Table 2: Prediction of total travel time using the long distance method (bus
line 523)

distance time | 7:00 - 10:00 | 10:00 - 14:00 | 14:00 - 19:00 | 19:00 - 23:00
---|---|---|---|---
MAE | STD | MAE | STD | MAE | STD | MAE | STD
1 - 3 | 56.72s | 66.06s | 59.97s | 69.67s | 59.07s | 64.44s | 54.82s | 51.11s
4 - 8 | 81.07s | 97.31s | 85.59s | 106.85s | 78.52s | 87.18s | 66.44s | 68.93s
9 - 15 | 107.01s | 121.46s | 101.93s | 129.62s | 105.89s | 113.19s | 85.65s | 87.78s
16 - 20 | 130.41s | 144.66s | 140.38s | 170.57s | 130.76s | 138.80s | 113.58s | 106.45s
21 - 27 | 153.70s | 161.04s | 170.72s | 205.03s | 155.09s | 152.63s | 172.30s | 135.29s

As can be seen, from Tables 1 and 2, for all distances longer than 4 stops,
the long distance method was more accurate than the recursive method.
Moreover, for distances longer than 8 stops, in most time intervals, the
results were more than twice as accurate. For 1 to 3 stops, results of both
methods were comparable.

Architecture comparison The goal of next experiment was to compare
effectiveness of different RBFN and MLP architects. Here, bus lines were
considered in eight groups with different characteristics (see, Section 3).
The goal was to to establish which method is best (MAE), and which is most
stable (STD). Overall, five MLP architectures and four RBFN architectures were
applied.

Training days included the same working days of the week from three
consecutive weeks (e.g. Wednesdays March, 10, 17 and 24). The test data was
from the following week (March 29 to April 2). In this way, five models were
created for a single bus line for each architecture. Only data from working
days was used, due to the difference in weekend schedules. Moreover, weekends
have different traffic patterns, and we wanted to exclude this factor. In
total, 29 bus lines and 9 different architectures were tried, which gives a
total of 1305 models (29 lines * 9 architectures * 5 days).

The RBFN architectures were implemented in Python using the Keras library, and
the RBF [15] code. The hidden layer used Gaussian radial basis function:
$exp(-\beta r^{2})$. The tested architectures have $M=10,15,25,35$ neurons in
the hidden layer (denoted as RBFN M). The results are presented in Table 3.

Table 3: RBFN-based prediction results

group model | RBFN 10 | RBFN 15 | RBFN 25 | RBFN 35
---|---|---|---|---
MAE | STD | MAE | STD | MAE | STD | MAE | STD
Center–Praga | 83.58s | 104.18s | 93.28s | 114.87s | 103.77s | 134.16s | 121.55s | 143.93s
Center–periphery | 114.46s | 145.65s | 115.23s | 154.1s | 103.76s | 139.74s | 104.87s | 135.96s
Long with in periphery North-South | 122.82s | 162.4s | 108.89s | 146.75s | 102.91s | 140.21s | 111.14s | 149.3s
Long with in periphery East–West | 118.28s | 159.77s | 123.06s | 162.73s | 113.49s | 152.23s | 109.46s | 149.66s
Long throught the centre | 98.13s | 127.49s | 100.1s | 127.63s | 95.91s | 122.0s | 100.29s | 117.7s
Express | 102.01s | 133.38s | 100.02s | 125.26s | 100.84s | 130.55s | 105.04s | 134.87s
Short with in center | 94.72s | 118.68s | 86.0s | 109.44s | 87.33s | 111.41s | 103.81s | 132.3s
Short with in periphery | 89.07s | 116.0s | 89.17s | 115.54s | 89.32s | 115.87s | 102.78s | 132.06s
All groups | 102.88s | 133.44s | 101.96s | 132.04s | 99.66s | 130.77s | 107.36s | 136.97s

It can be seen that for short routes (Center–Praga, short within Center, short
within periphery) and Express routes, the most accurate results were obtained
by the “smaller” RBFN architectures (RBFN 10 and RBFN 15). For the remaining
routes, RBFN 25 and RBFN 35 performed better. This is likely because, for long
paths, more diverse data “needs more neurons” in the hidden layer.

Overall, the RBFN 25 seems to be the “optimal architecture”. Only, for
Express, Short within center, and Short within periphery, its slightly worse
than RBFN 10 and RBFN 15. The difference in MAE was between 0.25s and 1.33s,
and in STD between 0.33s and 5.29s, which is not significant.

For MLP (implemented in Python with Keras) five architectures were tried:
networks with 2 or 3 hidden layers and: [6, 12], [12, 32], [6, 8, 12], [8, 8]
neurons. Here, ReLU and tanh were compared as activation functions. The
results are presented in Table 4.

Table 4: MLP-based prediction results

grupa model | MLP 6, 12 ReLU | MLP 12, 32 ReLU | MLP 12, 32 ReLU tanh | MLP 6, 8, 12 ReLU | MLP 8, 8 tanh
---|---|---|---|---|---
MAE | STD | MAE | STD | MAE | STD | MAE | STD | MAE | STD
Center–Praga | 112.33s | 143.9s | 96.43s | 119.94s | 94.08s | 118.30s | 114.5s | 147.78s | 115.93s | 146.75s
Center–periphery | 101.54s | 133.53s | 106.9s | 135.26s | 110.86s | 146.74s | 110.69s | 145.8s | 143.31s | 180.73s
Long within periphery North-South | 139.06s | 178.95s | 119.79s | 159.51s | 114.39s | 164.76s | 120.05s | 164.01s | 136.25s | 173.61s
Long within periphery East–West | 122.47s | 167.79s | 113.67s | 151.78s | 133.6s | 181.72s | 124.23s | 162.63s | 154.59s | 196.67s
Long through the centre | 109.35s | 146.14s | 97.71s | 128.49s | 118.17s | 153.16s | 105.89s | 139.59s | 109.16s | 143.21s
Express | 127.52s | 168.45s | 92.46s | 123.74s | 118.71s | 157.53s | 116.61s | 154.02s | 132.97s | 174.49s
Short within center | 97.02s | 130.73s | 95.42s | 121.13s | 108.07s | 139.06s | 106.31s | 132.19s | 101.55s | 135.38s
Short within periphery | 107.13s | 133.08s | 90.24s | 122.65s | 96.44s | 124.07s | 99.89s | 129.8s | 104.84s | 139.98s
All groups | 114.55s | 150.32s | 101.57s | 132.81s | 111.78s | 148.16s | 112.27s | 146.97s | 124.82s | 161.35s

Overall, the most effective and stable architecture for all groups was two
hidden layers consisting of 12 and 32 neurons, with ReLU as the activation
function. Note, however, that the results are somewhat more complex. The
network with the tanh activation function in the last layer, provided slightly
better results for the Center–Praga lines. Moreover, for the Long on the
north-south periphery lines it had value of MAE lower by 5.40s, but higher STD
value by 5.25s. On the other hand, the [6, 12] ReLU network, for the
Center–periphery lines, had MAE lower by 5.36s, and STD lower by 1.73s.
Overall, the [12, 32] ReLU network can be seen as the “winner”. Therefore it
was used in subsequent experiments. Finally, the MLP using only tanh achieved
worst results.

The next step was to compare RBFN 25 and MLP [12,32] ReLU. Based on the travel
time ratio (average travel time in a given hour/shortest travel time during
the day), two rush hour periods have been established: (1) the morning peak –
7:00 am to 10:00 am, (2) the afternoon peak – 3:00 pm to 7:00 pm. During peak
hours, the predictive models encounter problems to correctly determine travel
time. Therefore, a detailed comparison of the results of RBFN 25 and MLP [12,
32] ReLU was performed for the “worst-case conditions”.

Based on the analysis of obtained results it has been concluded that: (a) on
short distances in the morning peak, both models performed quite similarly.
Only in the case (b) Long at the the north-south periphery lines, the
effectiveness of the RBFN was better than the MLP. On the other hand, for (c)
Express routes, a significant advantage of the MLP was observed. For (d)
distances ranging from 6 to 15 stops, models are comparable, with minor
differences for Center – Praga and Short in the periphery, lines, where the
MLP has advantage. In the case of (e) long distances, from 16 to 25 stops, for
the Express lines, the RBFN model significantly exceeded the MLP.
Interestingly, (f) for the same group, the MLP performed significantly better
over short distances. For the (g) remaining long-distance groups, models
achieved similar performance, with a slight advantage of the RBFN for the four
groups and a slight advantage of the MLP for the other three groups. Similar
tendencies have been observed (h) for predictions during the afternoon peak.
Interestingly, (i) performance of models in the afternoon peak was better than
in the morning peak. This is probably due to the fact that there are more
buses on the route in the afternoon than in the morning hours, so the training
set contains more data. Hence, models could better capture specificity of
travel time during afternoon hours.

Overall, forecasting travel time, based only on the location data, the RBFN
and MLP provided similarly accurate forecasts, with a slight advantage of
RBFN. Note, however, that for long distances from 16 to 25 stops, in the
morning and afternoon rush hours, RBFN achieved significantly better
prediction accuracy for 3 groups, slightly better accuracy for 8 groups, and
for two groups, the accuracy of both models was almost equal. The MLP, for
this range, was much better only for the Center–Praga lines in the afternoon,
and slightly better for the Center–periphery and the Center–Praga lines in the
morning peak.

Separately, consider accuracy of the predictions of the best models, for each
bus group, in the most difficult conditions, i.e. over long distances during
rush hours. For 6 out of 8 groups, the accuracy of approximately 75% of
predictions had an error of less than 200s. Only for the Express and Long
within periphery West-East lines, the accuracy of about 75% of the predictions
had an error of less than 260s. Material presented in what follows aims at
improving the results obtained using only location information.

Next experiment used Historical Averages method (HA, [1, 5]), which applies
average travel times from previous days to estimate the current travel time.
Here, prediction depends on days of the week and rush hours. Hence, averages
are calculated separately for working days and weekends. Based on the GPS
data, the departure time from the last visited stop can be determined. Then,
based on historical data, the estimated time of arrival at the next stops is
determined. For each bus line, the data from the training sets (i.e. all
working days from March,8 to 26) was divided into 20-minute groups, i.e.
records sent every day between 10:00-10:20 belonged to one group. The average
travel times between current locations, and all stops to the end of the route
were calculated. For example, if for time period 10:00-10:20 data included 10
buses located on the 6-th stop of the line, the average times of travel this
stop to all stops, till the end of their lines, were calculated for these 10
buses. In this way, the HA algorithm approximated the estimated travel time
between any stops on the route for any given hour. The average travel times,
determined using this algorithm, have been stored as the training sets. Then,
for each bus line, for all test sets (i.e. data for working days from March,29
to April,2), travel time predictions were calculated using the HA algorithm.
In addition, analogous calculations were carried out while the average travel
times were calculated for the same day of the week only. For example, travel
time predictions for March 31 (Wednesday) were based on data data from March
10, 17, 24. Table 5 presents the MAE values of travel time predictions using
HA method.

Table 5: Prediction MAE for timetable, HA for same days, HA for working days.

group | time table | HA(the same day ) | HA(all days )
---|---|---|---
Long within periphery North-South | 79.22s | 81.53s | 71.94s
Long within periphery West-East | 80.85s | 75.02s | 62.38s
Center–periphery | 64.76s | 62.02s | 61.62s
Long through the centre | 71.58s | 68.75s | 67.64s
Express | 60.38s | 64.65s | 52.70s
Center–Praga | 64.56s | 52.89s | 50.22s
Short within periphery | 55.87s | 50.92s | 35.65s
Short within centre | 60.41s | 53.02s | 49.02s

As can be seen, HA based on data from all working days provides better
accuracy then when using data from the same day of the week. Moreover, observe
significant improvement in the accuracy of prediction using the HA algorithm,
compared use of travel time found in the timetables. Based on these
observations, the idea was born to replace the travel time predicted by
distribution data with travel time predicted by the HA algorithm in the NN-
based travel time prediction. This led to creation of a hybrid model. Here,
averages for all working days were used. The determined estimated travel times
replaced the estimated travel time from the timetable data feature, and the
last stop delay feature in the predictions with MLP and RBFN. The graphs 2, 3,
4

Comparison of experimental prediction accuracy of RBFN and MLP models, using
HA and not, brought the following conclusions. At distances from 1 to 5 stops,
models using distribution data were more accurate than predictions using HA,
for 6 of 8 groups: Long within periphery North–South (RBFN and MLP), Long
within periphery East–West (RBFN and MLP), Long within the center (RBFN and
MLP), Express (MLP), Short within periphery (MLP), Short within the center
(RBFN). For the remaining groups, the prediction errors were comparable. The
only exception was the MLP in the Center-periphery group, where the results
were more accurate when using HA data. On the other hand, for distances from 6
to 15 stops, in the chart 3 better results when using HA, were established for
all groups. The greatest improvement, when using HA, was noted for predictions
for long distances from 15 to 26 stops. In the 4 chart, you can see a
significant reduction in the prediction error of models using HA compared to
models using distribution data. Here, for Center–periphery, Center–Praga,
Short within periphery, Short within the center lines, for MLP, MAE decreased
almost twice. For the RBFN, drops of MAE were also significant, reaching up to
30s, for Long within periphery North-South lines.

Figure 2: Comparison of the prediction accuracy of the RBFN and MLP models
with the additional features of the predicted scheduled travel time and delay,
and of the models with the additional feature of the predicted travel time
determined by the HA algorithm. Results include predictions for distances from
1 to 5 stops during rush hour.

Figure 3: Comparison of the prediction accuracy of the RBFN and MLP models
with the additional features of the predicted scheduled travel time and delay,
and of the models with the additional feature of the predicted travel time
determined by the HA algorithm. Results include predictions for distances from
6 to 15 stops during rush hour.

Figure 4: Comparison of the prediction accuracy of the RBFN and MLP models
with the additional features of the predicted scheduled travel time and delay,
and of the models with the additional feature of the predicted travel time
determined by the HA algorithm. Results include predictions for distances from
16 to 25 stops during rush hour.

Based on these results it has been concluded that for 1 to 5 stops, RBFN or
MLP models, using the estimated travel time feature, calculated on the basis
of data from timetables, work better. However, models using HA-based data,
deliver slightly better performance for 6 to 15 stops. In the case of
predictions for 16 to 25 stops, the use of HA significantly reduces the
prediction error.

This suggests a “hybrid model” for prediction of bus travel time. Hence, in
Table 6, information, which models provide best prediction accuracy, depending
on the group of bus lines, and the travel distances is summarized.

Table 6: Models with best travel time prediction by distance and group of bus
lines.

group distance | up to 5 stops | 6 to 15 stops | more than 16 stops
---|---|---|---
Long within periphery | MLP + table time | MLP + HA | MLP + HA
Long within periphery East-West | MLP + time table | MLP + HA | MLP + HA
Center–periphery | MLP + table time | MLP + HA | MLP + HA
Long through the center | MLP + time table | MLP + HA | MLP + HA
Express | MLP + time table | MLP + HA | MLP + HA
Center–Praga | RBFN + time table | MLP + HA | MLP + HA
Short within periphery | MLP + time table | RBFN + HA | MLP + HA
Short within center | RBFN + time table | MLP + HA | MLP + HA

Based on results from Table 6, a hybrid model was developed, consisting of:
RBFN 25 or MLP [12, 32] with the features: predicted travel time based on
schedule data and delay at the last stop, and RBFN 25 or MLP [12, 32] with
estimated travel time using HA. When making predictions, depending to which
group given bus line belongs, and the distance for which the prediction is
performed, the hybrid approach uses models specified in Table 6.

The created hybrid model was compared with: MLP and RBFN models using a basic
set of features, the HA algorithm, and predictions based on data from
timetables. Figure 5 represents comparison of MAE values (in seconds) of the
predictions made by the created hybrid model with the previously mentioned
methods, depending on the distance (number of stops) for which the prediction
was made. Due to space limitation, only two, distinct, cases are reported. Bus
lines belonging to the same group had different lengths of routes. Therefore,
black triangles mark distances for which the prediction of travel time ends
for one of bus lines belonging to the group described by the graph. For this
reason, the graphs report rapid changes in MAE for adjacent distances, as in
the case of Long within periphery group, where lines have total of 30, 35, 36
and 38 stops.

Figure 5: Comparison of the MAE value of prediction methods at different
distances

Results represented in figures, and the remaining experimental results, show
that for short-distance predictions, for all groups, the HA algorithm and
distribution times delivered much better accuracy than the hybrid model.
However, for the groups Short within center, Short within periphery, Long
within periphery West East, Long within periphery North-South, Center-Prague,
for long distances (from 15 stops for Short within center, from 17 stops for
Short within periphery, from 27 stops for Long within periphery West East,
from 16 stops for Long within periphery North-South group, and from 16 stops
for Center-Prague group), the hybrid model was most accurate.

More specifically, for the Long within periphery North-South group, the hybrid
model provides the best accuracy for as much as 60% of distances. In the group
Long within periphery West-East, for travel time predictions from 27 stops,
where the predictions were made only for line 102 (only line this long), the
MAE of the hybrid model was 95-110s, while the MAE of the second best method
(HA), was 130-150s. In the case of the Center-periphery group, for distances
from 13 stops, the accuracy of the hybrid model was at the level of the time
table method. Here, the HA algorithm had the lowest MAE for all distances. In
the Long through center group, for prediction from 23 stops, the MAE of the
hybrid model was similar to the HA algorithm. For the Express route group,
apart from the travel time predictions for distances greater than 21 stops,
the MAE of the hybrid model was much higher than the MAE of the HA algorithm.
It is also worth noting that the accuracy of the basic RBFN and MLP was
significantly inferior to all other methods that have been experimented with.

##  5 Comparison with the models from the SotA

Let us now report the best MAE results found in the literature, and “compare”
it with the results reported here. The top 5 lowest MAE values are as follows.
(1) ANN from [2]; MAE at 40s; for predictions on distances of up to 16 stops.
(2) Hybrid model reported here; MAE at 67.17s; for predictions for the Short
through the center group, for time interval 19:00 to 23:00. (3) ANN from [1];
MAE at 70s. (4) BP from [12]; MAE at 125s; for prediction for morning rush
hours. (5) MLP from [5]; MAE at 138.40s; for prediction of travel time of
trams for noon hours.

However, it has to be stressed that these comparisons represent only a meta
level view of possible accuracy of travel time prediction. In each article,
results were obtained for bus (or tram) lines from different cities, where
each city has very different road and traffic structure and public transport
characteristics. For example, in more populated cities, or those with less
developed road infrastructure, there may be more delays due to heavy traffic,
which may influence the accuracy of model prediction. Moreover, Trondheim
(city population is around 200 thousands, while metropolitan population is
around 280000) is a much smaller city than others, which may influence
prediction accuracy. Further, Warsaw is split by a river, with limited number
of bridges, through which all bus lines have to pass to get “to the other
side”. Besides, each of analyzed bus/tram lines had different lengths.
Finally, no other work used data for multiple bus lines individually or in
groups. Regardless, it can be argued that the results reported in this
contribution are very competitive and worthy further explorations.

##  6 Concluding remarks

The aim of this work was to extend knowledge concerning applicability of
machine learning to prediction of travel times of public transport. The main
contributions of this work are related to the systematic study of
interrelations between characteristics of groups of bus lines and applied ML
algorithms. While it is extremely difficult to fairly compare results obtained
for cities located in India, Northern Europe and China, it is clear that the
obtained results are very competitive. Moreover, even though there was no
space to provide detailed results, performed experiments showed that use of
auxiliary features, e.g. number of busses moving between stops, or number of
crossings with lights does not lead to important improvement of accuracy of
prediction.

Pre-processed data is available at https://tinyurl.com/mvt4hu3f.

## References

  * [1] Wei Fan, Zegeye Gurmu, _Dynamic Travel Time Prediction Models for Buses Using Only GPS Data_ , International Journal of Transportation Science and Technology, 2015, 353-366.
  * [2] Erlend Dahl, Aleksander Sjåfjell, Simen Skogen, _On Implementations of Bus Travel Time Prediction Utilizing Methods in Artificial Intelligence_ , Norwegian University of Science and Technology, 2014
  * [3] Aleksander Sjåfjell, Erlend Dahl, Simen Skogen, _Intelligent Transportation Systems and Artificial Intelligence – a State of the Art Review_ , Norwegian University of Science and Technology, 2013
  * [4] Yongjie Lin, Xianfeng Yang, Nan Zou, Lei Jia, _Real-Time Bus Arrival Time Prediction: Case Study for Jinan, China_ , Journal of Transportation Engineering, 2013
  * [5] Adam Żychowski, Konstanty Junosza-Szaniawski, Aleksander Kosicki, _Travel time prediction for trams in Warsaw_ , CORES, 2017
  * [6] Yang, J.S., _Travel time prediction using the gps test vehicle and kalman filteringtechniques. In: Proceedings of the 2005_ , American Control Conference, 2005, IEEE(2005) 2128–2133
  * [7] Chun-Hsin Wu, Jan-Ming Ho, D. T. Lee, _Travel-Time Prediction With Support Vector Regression_ , IEEE Transactions on Intelligent Transportation Systems 5(4): 276-281, 2005
  * [8] Myung Jiwon, Dong-Kyu Kim, Seung-Young Kho, Chang-Ho Park, _Travel Time Prediction Using k Nearest Neighbor Method with Combined Data from Vehicle Detector System and Automatic Toll Collection System_ , Transportation Research Record Journal of the Transportation Research Board 2256:51-59, 2012
  * [9] Haykin, S., _Neural Networks: A Comprehensive Foundation. Prentice Hall PTR_ , Upper Saddle River, NJ, USA (1998)
  * [10] Hima Elsa Shaji Arun K. Tangirala, Lelitha Vanajakshi, _Prediction Of Trends In Bus Travel Time Using Spatial Patterns_ , Transportation Research Procedia 48 998–1007, 2020
  * [11] Antonio Comi, Mykola Zhuk, Volodymyr Kovalyshyn, Volodymyr Hilevych, _Investigating bus travel time and predictive models: a time series-based approach_ , Transportation Research Procedia 45 692–699, 2020
  * [12] Jia Fei, Ying Lu, Yunchuan Guo, Huibing Zhang, _Predicting Bus Arrival Time Using BP Neural Network and Dynamic Transfer_ , Procedia Computer Science 174 95-100, 2019
  * [13] Lei Wang, Zhongyi Zuo, Junhao Fu, _Bus Arrival Time Prediction Using RBF Neural Networks Adjusted by Online Data_ , Procedia - Social and Behavioral Sciences 138 67-75, 2014
  * [14] Lei Wang, Zhongyi Zuo, Junhao Fu, _Bus Dynamic Travel Time Prediction: Using a Deep Feature Extraction Framework Based on RNN and DNN_ , Electronics, 2020
  * [15] Vidnerová, Petra, _RBF-Keras: an RBF Layer for Keras Library._ , 2019, Available at https://github.com/PetraVidnerova/rbf\\_keras


# Table-Top Tunable Chiral Photonic Emitter

Lu Wang  State Key Laboratory of Magnetic Resonance and Atomic and Molecular
Physics, Wuhan Institute of Physics and Mathematics, Innovation Academy for
Precision Measurement Science and Technology, Chinese Academy of Sciences,
Wuhan 430071, China<EMAIL_ADDRESS>Marcelo Fabián Ciappina
Department of Physics, Guangdong Technion - Israel Institute of Technology,
Shantou 515063, China  Technion - Israel Institute of Technology, Haifa 32000,
Israel  Guangdong Provincial Key Laboratory of Materials and Technologies for
Energy Conversion, Guangdong Technion - Israel Institute of Technology,
Shantou 515063, China  Thomas Brabec  Department of Physics, University of
Ottawa, Ottawa, Ontario K1N 6N5, Canada  Xiaojun Liu  State Key Laboratory of
Magnetic Resonance and Atomic and Molecular Physics, Wuhan Institute of
Physics and Mathematics, Innovation Academy for Precision Measurement Science
and Technology, Chinese Academy of Sciences, Wuhan 430071, China

###### Abstract

The increasing interest in chiral light stems from its spiral trajectory along
the propagation direction, facilitating the interaction between different
polarization states of light and matter. Despite tremendous achievements in
chiral light-related research, the generation and control of chiral pulse have
presented enduring challenges, especially at the terahertz and ultraviolet
spectral ranges, due to the lack of suitable optical elements for effective
pulse manipulation. Conventionally, chiral light can be obtained from
intricate optical systems, by an external magnetic field, or by metamaterials,
which necessitate sophisticated optical configurations. Here, we propose a
versatile tunable chiral emitter, composed of only two planar Weyl semimetals
slabs, addressing the challenges in both spectral ranges. Our results open the
way to a compact tunable chiral emitter platform in both terahertz and ultra-
violet frequency ranges. This advancement holds the potential to serve as the
cornerstone for integrated chiral photonics.

Controlling the interplay of light and its interaction with matter has long
been a fundamental objective in optics. In particular, chiral light has
attracted substantial research interests, as the light polarization spiraling
along the propagation direction provides a variety of opportunities for
control and spectroscopy of matter. To name a few, chiral light is applied in
subcycle optical gates such as angular-dependent dynamics studies Bloch  _et
al._ (2021); Chen  _et al._ (2017), chiral molecule spectroscopy Tanaka  _et
al._ (2018); Richardson and Riehl (1977), mechanical separation of chiral
objects Canaguier-Durand  _et al._ (2013), and chiral quantum optics Lodahl
_et al._ (2017).

Despite the tremendous achievements in chiral light related research, the
generation and control of ultraviolet and near-infrared, chiral pulses have
continued to present challenges for the following reasons. Firstly, although
chiral objects are widespread in nature, such as DNA and certain proteins Mun
_et al._ (2020), the chiral light-matter interaction is typically very weak
due to the large mismatch between their atomic feature sizes and optical
wavelengths Forbes (2022); Mun  _et al._ (2020). Although one can use external
magnetic fields to enhance the material chirality, the typically required
magnetic field strengths are of a few Tesla Fan  _et al._ (2021); Stanciu  _et
al._ (2007), leading to cumbersome and costly magnetic structures that are
unsuited to large-scale manufacturing. On the other hand, metasurfaces and
plasmonics Hendry  _et al._ (2012); Zhang  _et al._ (2022); Mun  _et al._
(2020) can be utilized for manipulating chiral light. However, the necessary
complex optical configurations require nano-scale fabrication with multiple
lithographic processes Santos  _et al._ (2015); Chen and Pepin (2001).

Generation and control of chiral light in the ultra-violet frequency range
(800 THz-30 PHz) is even more difficult, due to the lack of suitable optical
elements Borrego-Varillas  _et al._ (2018); Brahms and Travers (2023). This
impedes progress in many areas, from high-resolution spectroscopy Perkampus
(2013) and photo-lithography del Barrio and Sánchez-Somolinos (2019), to
controlling atomic/ molecular electronic transitions Ranitovic  _et al._
(2014) and inactivation of microorganisms Li  _et al._ (2009).

Control and manipulation of terahertz (0.1 THz-30 THz) chiral pulses, on the
opposite end of the optical spectrum, have proven equally challenging Fan  _et
al._ (2021). Terahertz pulses can serve as an active gate to enhance the
spatio-temporal resolution of spectroscopy Cocker  _et al._ (2013), realize
relativistic electron acceleration Hibberd  _et al._ (2020), and be the key to
6G communication Rikkinen  _et al._ (2020). Nevertheless, the terahertz
frequency can not be accessed easily, because it is in between the standard
electronic and photonic manipulation methods Vitiello and Tredicucci (2021);
Pang  _et al._ (2022). The difficulties in handling Terahertz pulses are
mainly due to two reasons. One is the large diffraction owing to their long
wavelength nature. The second aspect involves the absorption within solids, as
this frequency range is in resonance with lattice phonons Kozina  _et al._
(2019). Consequently, minimizing interactions with solid-state optical
elements is imperative to circumvent these challenges. The above difficulties
have been referred to as the "terahertz gap" Vitiello and Tredicucci (2021);
Pang  _et al._ (2022).

Figure 1: Elliptically polarized harmonics induced by linearly polarized
driving field along $y$ dimension. Panel a illustrates the proposed planar
chiral emitter i.e. two Weyl semimetal (WSM) slabs separated by air. Panels
b,c (d,e) are the reflected (transmitted) spectra of the harmonics as a
function of the air separation distance $L$ for $y-$ and $z-$ polarization,
respectively. The color bar represents $\log{[E(\omega)]}$. Panels f-i and j-m
show the 3D temporal evolution of given harmonic orders (both transmitted and
reflected) with $L=150$ nm and $L=450$ nm, respectively. These two air gap
distances are indicated by horizontal dashed lines in b-e. The results are
calculated with driving pulse parameters as, pulse duration $\tau=10$ fs,
center wavelength $\lambda_{0}=1\,\mu$m, and peak field strength
$E_{0}=2\times 10^{8}\,\text{V}/\text{m}$ corresponding to intensity $\sim
10^{14}\text{W}/\text{m}^{2}$. Both the WSM slabs are of 100 nm thickness. The
Weyl nodes separation is 2b along $k_{x}$ dimension where $b=0.06\pi/a_{0}$,
$a_{0}=3.4\r{A}$ is the lattice parameter along $x$ Lv  _et al._ (2015). The
Fermi energy is $100$ meV. The harmonics are at the ultra-violet spectral
range.

In this work, we introduce a versatile tunable chiral emitter addressing the
challenges in both the ultra-violet and terahertz spectral ranges at the same
time. The proposed setup employs only two planar Weyl semimetal (WSM) slabs
and a linearly polarized field driving high harmonic generation (HHG), as
illustrated in Fig. 1a. Our approach utilizes the chirality of the high-
harmonic generation process in the WSM planar structures, which presents the
basis for tunable chiral ultraviolet emission from a table-top source.
Elimination of complex optical elements simplifies the experimental
realization of the proposed setup. The theoretical model used in our analysis
solves the two-band equations for semimetals and Maxwell’s equations self
consistently. The model applies to WSMs, Dirac semimetals, and 2D semimetals
such as graphene.

WSMs are a newly discovered class of quantum materials with unique optical
properties and an extremely strong nonlinearity Patankar  _et al._ (2018); Wu
_et al._ (2017). In the momentum space, the points where valence and
conduction bands cross are Weyl nodes, around which electrons behave as
massless Weyl fermions Yan and Felser (2017). The Weyl nodes appear in pairs
with opposite chirality, making them the perfect source and detector for
chiral light Ma  _et al._ (2017); Yoshikawa  _et al._ (2022). An inherent
feature of the chiral properties in WSMs is the intrinsic coupling of
different polarization states. In other words, by choosing the right
combination of the driving field’s polarization direction with respect to the
Weyl nodes separation direction, it is possible to achieve elliptically
polarized emission using a purely linearly polarized driving field, in the
absence of any external magnetic field Wang  _et al._ (2023); Wang (2023);
Avetissian  _et al._ (2022); Gao  _et al._ (2020).

When exposed to a strong laser field, the nonlinear currents arise from the
movement of electrons within the WSM through two primary pathways: the motion
within the individual bands (intra-band) and the polarization accumulation
between bands (inter-band) Ghimire and Reis (2019); Wu  _et al._ (2015); Vampa
_et al._ (2014). The nonlinear current emits photons with energies that are
integer multiples of the driving field. This process is known as high-harmonic
generation in solids Ghimire and Reis (2019); Wu  _et al._ (2015); Vampa  _et
al._ (2014). In momentum space, the Weyl points act like magnetic monopoles.
The non-vanishing Berry curvature flows from the source $\chi=1$ to the sink
$\chi=-1$. These magnetic monopoles create magnetic field equivalent effects.
Under the influence of the magnetic field, the electrons move toward a
direction perpendicular to the applied field and velocity, resulting in
harmonic emissions with more than one polarization. Details can be found in
supplementary material (SM) sec. III sup . The generated harmonic pulses
exhibit excellent coherence properties (Salieres  _et al._ , 1997). Combining
the elliptical emission from the WSM and the two-slab structure, the proposed
scheme realizes easy tunability of the chiral pulses via interference along
the propagation direction. This is similar to the principles of a Fabry–Pérot
cavity. Moreover, since the harmonics are sensitively related to the
properties of the driving field, the chiral/3D polarization states of the
high-harmonic radiation can be fully controlled through the intense pulse
laser.

Results.
Modifying chiral harmonic properties by tuning the distance $L$. In Figs.
1b-e, the harmonic spectra, both reflected and transmitted, are depicted as a
function of the separation distance between the two slabs, denoted as $L$. The
theoretical method is briefly discussed below and in more detail in the SM
Sec. III sup . The calculations are based on standard experimental available
parameters. The driving field strength is $E_{0}=2\times
10^{8}\text{V}/\text{m}$ Ghimire  _et al._ (2011), the WSM slab thickness is
assumed to be (100 nm)Chi  _et al._ (2018); Bedoya-Pinto  _et al._ (2020), and
standard WSM material parameters are used Lv  _et al._ (2015).

Since the model we chose to describe WSM breaks only the time-reversal
symmetry and preserves the inversion symmetry, only odd harmonics are present
(Methods sections Eqs.2-4). The two-slab structure forms a Fabry–Pérot
interference cavity. By alternating the distances between the two WSM slabs,
the boundary conditions change for the electric field (i.e. superposition of
the harmonics and the driving field) interacting with the structure, resulting
in modified harmonics. Visible spectral amplitude modification can be observed
(note that the color bar is in log scale). The harmonic spectra are expected
to exhibit periodicity concerning $L$, with a period of $\sim\lambda_{0}/2$,
where $\lambda_{0}=1\,\textmu$m is the center wavelength of the driving field.
The spectra observed in reflection (Figs. 1b,c) and transmission (Figs. 1d,e)
manifest distinctive characteristics. For the reflected wave, the dominant
emission originates directly from the immediate reflection, i.e., the wave
reflected directly from the surface. Consequently, the spectra in reflection
demonstrate weaker nonlinear effects, specifically lower-order harmonics, in
contrast to the transmitted case, where more pronounced harmonic orders are
evident due to the interaction of the driving pulse through two WSM slabs. It
can be seen that our compact chiral light source offers the potential for
generating diverse polarization states across various harmonic orders (see
Fig.1f-m). These exotic polarization states can trigger nuanced angle-resolved
electron dynamics that will enhance our understanding of light-matter
interaction with a sub-cycle temporal resolution Jiang  _et al._ (2022).

Polarization and ellipticity evolution in time.

Figure 2: The polarisation and ellipticity evolution for transmitted harmonics
with the air gap distance $L=450\,$nm. In the upper panel, the solid and
dashed curves depict the amplitude of the electric field for the specified
harmonic in the $y$ and $z$ dimensions. The blue and red curves represent the
ratio of the amplitude $E_{0z}/E_{0y}$, and the phase difference
$\phi_{z}-\phi_{y}$, respectively. The 3D temporal evolution of b and c are
presented in Figs. 1k,m.

To further characterize the harmonic pulses and their temporal evolution, we
examine the polarization state of a given harmonic as a function of time. When
the driving pulse is polarized along the $y$ direction, the emitted harmonics
are exclusively along the $y$ and $z$ directions (SM Sec. IV sup ). The
ellipticity of a pulse is determined by both the relative amplitude and phase
relations of the two orthogonal polarizations. For instance, the electric
field of a particular harmonic order, polarized along a given direction can be
expressed as $E_{i}(t)=\\{E_{0i}(t)\exp[i\phi_{i}(t)]+c.c\\}/2$, where
$i\in{y,z}$ and $E_{0i},\phi_{i}\in\Re$. For a linearly polarized pulse, it
necessitates $\phi_{z}-\phi_{y}=n\pi$, where $n=0,\pm 1,\pm 2,...$. Similarly,
a circularly polarized light demands $E_{0z}/E_{0y}=1$ and
$\phi_{z}-\phi_{y}=\pi/2+n\pi$, where $n=0,\pm 1,\pm 2,...$. In Figs. 2a-c we
depict the 1st, 3rd, and 5th harmonics, respectively, in the transmitted
direction with an air gap distance $L=450\,$nm. Figures 1k,m and 2b,c
represent the same electric fields. In the upper panel of Fig. 2, the solid
(dashed) black curve represents the electric field envelope of $E_{0y}$
($E_{0z}$). To quantify the pulse evolution dynamics, we define two
parameters: the ratio of the amplitude $E_{0z}/E_{0y}$, and the phase
difference $\phi_{z}-\phi_{y}$. These two parameters are indicated by blue and
red curves on the lower panels, respectively. The temporal evolution of the
harmonic fields displays distinctive and dynamic changes throughout the entire
pulse duration. In other words, for a given harmonic order, the ellipticity
varies with respect to time within the pulse. Since the multi-photon
ionization is highly sensitive to the driving pulse polarisation and carrier-
envelope phase, these unusual and exotic polarization changes of the generated
harmonics across a single pulse can be harnessed for ultrafast electron
dynamics control Eickhoff  _et al._ (2021); Song  _et al._ (2023); Li and He
(2021). In addition, these time-varying polarisation states of the fields may
bring new symmetry to the free electron wave packets via the photon electron
ionization, resulting in a powerful tool for coherent control of the quantum
process Kerbstadt  _et al._ (2019).

Influence of the Fermi Energy on Carrier Envelope Phase (CEP) at the Terahertz
regime.

Figure 3: The color plot a, b represent the current normalized current
$J_{y}(t)$, $J_{z}(t)$, respectively. The results are calculated by
$\lambda_{0}=10\,\textmu$m corresponds to frequency $f=10\,$THz and $\hbar
2\pi f\approx 120\,$meV, pulse duration $\tau=40\,$ fs, electric field is
polarized along $\hat{y}$ with peak field strength $E_{0}=2\times 10^{7}V/m$.
The transmitted electric fields after two Weyl semimetal slabs (100 nm each)
separated by 2µm in the air are presented in c,d where the Fermi energies are
$0$ and $500$ meV, respectively.

The proposed scheme in Fig. 1a is also valid for the terahertz frequency
range. Thus, we will not repeat the same calculations. Instead, we analyze how
the waveform of the terahertz pulse varies as a function of the Fermi energy
$\text{E}_{\text{f}}$. Since terahertz pulses are frequently applied to
electron acceleration and communication, the capacity to tailor their waveform
is very important. With increasing Fermi energy, the inter-band transition is
suppressed, so that the intra-band transition starts to dominate. Since the
inter- and intra-band transitions are not always in phase, the final induced
current has a phase shift (see Figs. 3a,b). The peak phase differences of the
current are defined as $\Delta=(t_{y}-t_{z})\omega_{0}$ where $t_{i}$ is the
time corresponding to the peak of the electric field polarized along the $i$
direction with $i\in\\{y,z\\}$. The value of $\Delta$ as a function of
$\text{E}_{\text{f}}$ is plotted as the dotted curve in Fig. 3a with the
horizontal axis on top. The transmitted electric fields with
$\text{E}_{\text{f}}=0$ and $\text{E}_{\text{f}}=500$ meV are presented in
Fig. 3c and d, respectively. The two Weyl semimetal slabs (100 nm each) are
separated by 2µm in the air.

Discussion and Conclusion
Several WSMs have been experimentally discovered, such as
$\text{Y}_{2}\text{Ir}_{2}\text{O}_{7}$, $\text{HgCr}_{2}\text{Se}_{4}$, TaAs,
TaP and $\text{Co}_{3}\text{Sn}_{2}\text{S}_{2}$ Kumar  _et al._ (2017); Wan
_et al._ (2011); Xu  _et al._ (2011); Yan and Felser (2017); Okamura  _et al._
(2020), which could be used for an experimental implementation of our scheme.
Furthermore, previous experiments have successfully grown WSMs from
nanostructures Bachmann  _et al._ (2017) to crystals, grown by the chemical
vapor transfer method Ma  _et al._ (2017). They are accessible within a
considerable size range, spanning a few nanometers Bedoya-Pinto  _et al._
(2020) to millimeters Chi  _et al._ (2018).

The compact chiral light source we propose opens the door to various
applications. For example, it creates a new dimension for information coding
in optical communication Forbes (2022). Besides, since chiral molecules
respond differently to chiral light, chiral light sources can be used to sort
molecules and achieve all-optical chiral discrimination Genet (2022); Zhang
_et al._ (2017). Furthermore, this easily tunable setup can be integrated with
chemical experiments and serve as a sensitive trigger for photochemical
reactions Rikken and Raupach (2000).

Our methodology also has potential applications in the realm of band structure
spectroscopy of WSMs. For example, the chirality of the emitted harmonics is
determined by the Weyl nodes separation $b$, given the fact that the emission
is only linearly polarized if $b=0$ (see Method section and SM Sec. IV sup ).
As a result, emissions observed along the $\hat{z}$ direction become a pivotal
metric for deducing the Weyl nodes’ separation, $b$. Besides, since the
nonlinearity of the emitted harmonics is proportional to the Fermi velocities
Wang  _et al._ (2022) as shown in Eqs. (2-4), the Fermi velocities of WSMs can
be retrieved by obtaining the amplitude of a given harmonic.

Method. To achieve the proposed objective mentioned above, we choose the most
straightforward model to describe the WSM. In this approach, the time-reversal
symmetry is broken whereas the inversion symmetry is preserved. We focus on a
two-level system with one pair of Weyl nodes of opposite chirality separated
by a distance 2$b$ along the $\hat{k}_{x}$ dimension in the momentum space.
Specifically, the Hamiltonian is written as Avetissian  _et al._ (2022); Hasan
_et al._ (2017); Armitage  _et al._ (2018); Ishikawa (2010); Lim  _et al._
(2020):

$H_{\chi}=\chi\sum_{i\in\\{x,y,z\\}}v_{i}\sigma_{i}\pi_{i},\quad\chi=\pm 1,$
(1)

where $x,\,y,\,z$ are Cartesian coordinates, $v_{x}=v_{y}=v_{z}=1\times
10^{6}\,\text{m}/\text{s}$ are the Fermi velocities Armitage  _et al._ (2018),
$\sigma_{i}$ with $i\in\\{x,\,y,\,z\\}$ define the Pauli matrices, and $\chi$
represents the chirality of the node. The kinetic momentum of the electron
under the influence of an external field (for a detailed derivation see SM
Sec. III sup ) is described by $\pi_{i}=\hbar(k_{i}-\chi b_{i})-qA_{i}$, where
$\hbar k_{i}$ is the electron momentum in the absence of an external field,
$q=-e$, with $e=\lvert e\rvert$ the elementary charge, and $\bm{b}=(b,0,0)$ is
the node off-set with respect to the origin in the momentum space. The vector
potential is defined as
$A_{i}=-\int_{-\infty}^{t}E_{i}(t^{\prime})dt^{\prime}$, where $E_{i}$ is the
electric field polarized along $i$. The total current along the $i$ direction
takes into account the contributions from both nodes, i.e.
$J_{i}(t)=J_{i}(\chi=+1,t)+J_{i}(\chi=-1,t),i\in\\{x,y,z\\}$ Avetissian  _et
al._ (2022).

The general form of the current can be found in SM Sec. III Eqs. S34-S36 sup .
These expressions are highly comprehensive. Setting $b=0$ results in the
expression for Dirac semimetals Lim  _et al._ (2020). Based on this, by
further eliminating the third dimension $z$, one obtains the expression for
graphene Ishikawa (2010).

In our configuration, the Weyl nodes separation is along $\hat{x}$ and the
driving field is linearly polarized along $\hat{y}$. Applying the frozen band
approximation and accounting for k-space symmetry (refer to SM Sec. IV sup ),
we simplify the analytical expressions for the current as the following:

$\displaystyle j_{x}(t)=0$ (2) $\displaystyle
j_{y}(t)=-2\chi\mathcal{M}_{0}\frac{\Omega_{0}}{\omega_{0}}\text{Im}\left[\sum_{\begin{subarray}{c}n,m=-\infty\\\
n-m\in\text{
odd}\end{subarray}}^{\infty}\mathcal{C}_{2}(N,n){C}_{3}(N,m)^{*}\large\right]$
(3) $\displaystyle
j_{z}(t)=\frac{\sin(\theta)}{\cos(\phi)}\mathcal{M}_{0}\sum_{\begin{subarray}{c}n,m=-\infty\\\
n-m\in\text{
odd}\end{subarray}}^{\infty}\mathcal{C}_{1}(N,n)\mathcal{C}_{1}(N,m)^{*}$
$\displaystyle-2vq\frac{\sin(\theta)^{2}}{\cos(\theta)}\text{Re}\left[\sum_{\begin{subarray}{c}n,m=-\infty\\\
n-m\in\text{
odd}\end{subarray}}^{\infty}\mathcal{C}_{1}(N,n)\mathcal{C}_{3}(N,m)^{*}\right].$
(4)

The related variables are defined as

$\displaystyle\mathcal{C}_{1}(N,n)=\frac{n\text{J}_{n}(N)\exp(i\omega_{0}n)}{n+\Omega_{0}/\omega_{0}}\exp(-\frac{t^{2}}{\tau^{2}}),$
$\displaystyle\mathcal{C}_{2}(N,n)=\frac{\text{J}_{n}(N)\exp(i\omega_{0}n)}{n+\Omega_{0}/\omega_{0}}\exp(-\frac{t^{2}}{\tau^{2}}),$
$\displaystyle\mathcal{C}_{3}(N,n)={\text{J}_{n}(N)\exp(i\omega_{0}n)},$
$\displaystyle\mathcal{M}_{0}=qv\frac{\sin(\theta)[\cos(\phi)^{2}+\sin(\phi)^{2}\cos(\theta)^{2}]}{\cos(\phi)\cos(\theta)},$
$\displaystyle N=\frac{-\chi
q\exp(-{t}^{2}/\tau^{2})\cos{\phi}\cos{\theta}E_{0}}{({2\omega_{0}\hbar|k|})}$

where $\text{J}_{n}$ is the Bessel function of the first kind,
$\Omega_{0}=4v|q|E_{0}/\hbar\pi\omega_{0}$, $\theta=\arccos(k_{z}/|k|)$,
$|k|=\sqrt{(k_{x}-\chi b)^{2}+k_{y}^{2}+k_{z}^{2}}$, and
$\phi=\arctan[k_{y}/(k_{x}-\chi b)]$. It can be seen that only odd harmonics
survive (see SM Sec. IV for details). In particular, Eq.(4) is proportional to
the Weyl nodes separation $b$. If the two Weyl nodes of the opposite chirality
merge to one place in k space (i.e. $b=0$), $j_{z}(t)=0$ [see SM Sec. IV
Eq.(S84)].

The planar structure of the WSM slabs is modeled via a customized Finite-
Difference Time-Domain (FDTD) method. The emission current (we have used the
general form as shown in Eqs. S34-S36) is integrated into the Maxwell
equations. The numerical details are presented in SM Secs. VI and VII sup .
The total current requires integration over the 3D $k$-space. The most
standard mesh point selection would be an even-distance mesh, where up to 500
mesh points along each dimension are typically required Avetissian  _et al._
(2022). We choose the parameters for one of the most standard Weyl semimetals
i.e. TaAs. The lattice parameters at 300 K are $a_{0}=b_{0}=3.4\r{A}$ and
$c_{0}=11.6\r{A}$ Lv  _et al._ (2015). We choose to integrate over the entire
Brillouin space with the integration range as $[-\pi/r,\pi/r]$, where
$r\in\\{a_{0},b_{0},c_{0}\\}$. For very large $k$ values, there is almost no
transition Ishikawa (2010). This can be understood intuitively. With a very
large $k_{z}$ for example, in the $k_{x}-k_{y}$ plane, the Weyl semimetal
behaves as a material with a very large band gap ($\sim 2\hbar v_{z}\lvert
k_{z}\rvert$). Thus, the inter-band transitions are harder. Since the dominant
inter-band transitions are around the Weyl nodes, in the entire Brillouin
zone, many mesh points can be ignored, i.e., the current is almost zero in
regions far from the Weyl nodes (SM Sec. V sup ). Consequently, we use non-
equal-distance mesh sizes, which are denser around the cone (SM Sec. V sup
Fig.S7). With this approach, we obtained comparable results as Avetissian  _et
al._ (2022) (details see Sec. VI Fig. S9) using a mesh size of $60\times
60\times 60$, which is $\sim 10^{3}$ times smaller than the conventional
method Avetissian  _et al._ (2022); Vampa  _et al._ (2014). For our work, one
single job of the high-harmonic generation via the 2 Weyl semimetal slabs
(mesh size 60$\times 60\times 60$) requires a computational time of $\sim 12$
hours on a single node with 128-core CPU @ 2.6GHz.

Data availability
Data supporting the findings of this work are available from the corresponding
authors upon reasonable request. Code availability
The code developed for this work is available to the public at the repository:
Acknowledgements
We would like to thank Beijing Beilong Super Cloud Computing Co., Ltd for
their support in HPC. In particular, L.W. thanks Mr. X.L. Han, Mr. R. Liu, and
Mr. H.B. Chen for their timeless help. This work is supported by the National
Key Research and Development Program of China (Grant No. 2019YFA0307702) and
the National Natural Science Foundation of China (Grant No. 12121004). M.F.C
acknowledges financial support from the Guangdong Province Science and
Technology Major Project (Future functional materials under extreme conditions
- 2021B0301030005) and the Guangdong Natural Science Foundation (General
Program project No. 2023A1515010871).
Author contributions
Competing interests
The authors declare no competing interests.
Additional information
Supplementary information: The online version contains supplementary material
available at
Correspondence and requests for materials should be addressed to Lu Wang.

## References

  * Bloch  _et al._ (2021) E. Bloch, S. Larroque, S. Rozen, S. Beaulieu, A. Comby, S. Beauvarlet, D. Descamps, B. Fabre, S. Petit, R. Taïeb, _et al._ , Physical Review X 11, 041056 (2021).
  * Chen  _et al._ (2017) K. Chen, Y. Feng, L. Cui, J. Zhao, T. Jiang, and B. Zhu, Scientific reports 7, 42802 (2017).
  * Tanaka  _et al._ (2018) H. Tanaka, Y. Inoue, and T. Mori, ChemPhotoChem 2, 386 (2018).
  * Richardson and Riehl (1977) F. S. Richardson and J. P. Riehl, Chemical Reviews 77, 773 (1977).
  * Canaguier-Durand  _et al._ (2013) A. Canaguier-Durand, J. A. Hutchison, C. Genet, and T. W. Ebbesen, New Journal of Physics 15, 123037 (2013).
  * Lodahl  _et al._ (2017) P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeutel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller, Nature 541, 473 (2017).
  * Mun  _et al._ (2020) J. Mun, M. Kim, Y. Yang, T. Badloe, J. Ni, Y. Chen, C.-W. Qiu, and J. Rho, Light: Science & Applications 9, 139 (2020).
  * Forbes (2022) A. Forbes, Science 377, 1152 (2022).
  * Fan  _et al._ (2021) F. Fan, C. Zhong, Z. Zhang, S. Li, and S. Chang, Nanoscale Advances 3, 4790 (2021).
  * Stanciu  _et al._ (2007) C. D. Stanciu, F. Hansteen, A. V. Kimel, A. Kirilyuk, A. Tsukamoto, A. Itoh, and T. Rasing, Physical review letters 99, 047601 (2007).
  * Hendry  _et al._ (2012) E. Hendry, R. Mikhaylovskiy, L. Barron, M. Kadodwala, and T. Davis, Nano letters 12, 3640 (2012).
  * Zhang  _et al._ (2022) X. Zhang, Y. Liu, J. Han, Y. Kivshar, and Q. Song, Science 377, 1215 (2022).
  * Santos  _et al._ (2015) A. Santos, M. Deen, and L. Marsal, Nanotechnology 26, 042001 (2015).
  * Chen and Pepin (2001) Y. Chen and A. Pepin, Electrophoresis 22, 187 (2001).
  * Borrego-Varillas  _et al._ (2018) R. Borrego-Varillas, L. Ganzer, G. Cerullo, and C. Manzoni, Applied Sciences 8, 989 (2018).
  * Brahms and Travers (2023) C. Brahms and J. C. Travers, Optics Letters 48, 151 (2023).
  * Perkampus (2013) H.-H. Perkampus,  _UV-VIS Spectroscopy and its Applications_ (Springer Science & Business Media, 2013).
  * del Barrio and Sánchez-Somolinos (2019) J. del Barrio and C. Sánchez-Somolinos, Advanced Optical Materials 7, 1900598 (2019).
  * Ranitovic  _et al._ (2014) P. Ranitovic, C. W. Hogle, P. Rivière, A. Palacios, X.-M. Tong, N. Toshima, A. González-Castrillo, L. Martin, F. Martín, M. M. Murnane, _et al._ , Proceedings of the National Academy of Sciences 111, 912 (2014).
  * Li  _et al._ (2009) D. Li, S. A. Craik, D. W. Smith, and M. Belosevic, Water Research 43, 481 (2009).
  * Cocker  _et al._ (2013) T. L. Cocker, V. Jelic, M. Gupta, S. J. Molesky, J. A. Burgess, G. D. L. Reyes, L. V. Titova, Y. Y. Tsui, M. R. Freeman, and F. A. Hegmann, Nature Photonics 7, 620 (2013).
  * Hibberd  _et al._ (2020) M. T. Hibberd, A. L. Healy, D. S. Lake, V. Georgiadis, E. J. Smith, O. J. Finlay, T. H. Pacey, J. K. Jones, Y. Saveliev, D. A. Walsh, _et al._ , Nature Photonics 14, 755 (2020).
  * Rikkinen  _et al._ (2020) K. Rikkinen, P. Kyosti, M. E. Leinonen, M. Berg, and A. Parssinen, IEEE Communications Magazine 58, 22 (2020).
  * Vitiello and Tredicucci (2021) M. S. Vitiello and A. Tredicucci, Advances in Physics: X 6, 1893809 (2021).
  * Pang  _et al._ (2022) X. Pang, O. Ozolins, S. Jia, L. Zhang, R. Schatz, A. Udalcovs, V. Bobrovs, H. Hu, T. Morioka, Y.-T. Sun, _et al._ , Journal of Lightwave Technology 40, 3149 (2022).
  * Kozina  _et al._ (2019) M. Kozina, M. Fechner, P. Marsik, T. van Driel, J. M. Glownia, C. Bernhard, M. Radovic, D. Zhu, S. Bonetti, U. Staub, _et al._ , Nature Physics 15, 387 (2019).
  * Lv  _et al._ (2015) B. Lv, N. Xu, H. Weng, J. Ma, P. Richard, X. Huang, L. Zhao, G. Chen, C. Matt, F. Bisti, _et al._ , Nature Physics 11, 724 (2015).
  * Patankar  _et al._ (2018) S. Patankar, L. Wu, B. Lu, M. Rai, J. D. Tran, T. Morimoto, D. E. Parker, A. G. Grushin, N. Nair, J. Analytis, _et al._ , Physical Review B 98, 165113 (2018).
  * Wu  _et al._ (2017) L. Wu, S. Patankar, T. Morimoto, N. L. Nair, E. Thewalt, A. Little, J. G. Analytis, J. E. Moore, and J. Orenstein, Nature Physics 13, 350 (2017).
  * Yan and Felser (2017) B. Yan and C. Felser, Annual Review of Condensed Matter Physics 8, 337 (2017).
  * Ma  _et al._ (2017) Q. Ma, S.-Y. Xu, C.-K. Chan, C.-L. Zhang, G. Chang, Y. Lin, W. Xie, T. Palacios, H. Lin, S. Jia, _et al._ , Nature Physics 13, 842 (2017).
  * Yoshikawa  _et al._ (2022) N. Yoshikawa, K. Ogawa, Y. Hirai, K. Fujiwara, J. Ikeda, A. Tsukazaki, and R. Shimano, Communications Physics 5, 328 (2022).
  * Wang  _et al._ (2023) L. Wang, F. J. G. de Abajo, G. T. Papadakis, _et al._ , Physical Review Research 5, L022051 (2023).
  * Wang (2023) L. Wang, Advanced Photonics Research , 2300194 (2023).
  * Avetissian  _et al._ (2022) H. Avetissian, V. Avetisyan, B. Avchyan, and G. Mkrtchian, Physical Review A 106, 033107 (2022).
  * Gao  _et al._ (2020) Y. Gao, S. Kaushik, E. Philip, Z. Li, Y. Qin, Y. Liu, W. Zhang, Y. Su, X. Chen, H. Weng, _et al._ , Nature communications 11, 720 (2020).
  * Ghimire and Reis (2019) S. Ghimire and D. A. Reis, Nature physics 15, 10 (2019).
  * Wu  _et al._ (2015) M. Wu, S. Ghimire, D. A. Reis, K. J. Schafer, and M. B. Gaarde, Physical Review A 91, 043839 (2015).
  * Vampa  _et al._ (2014) G. Vampa, C. McDonald, G. Orlando, D. Klug, P. Corkum, and T. Brabec, Physical review letters 113, 073901 (2014).
  * (40) Http://link.naturephotonics.xx.xxxxxx for calculation details.
  * Salieres  _et al._ (1997) P. Salieres, A. L’huillier, P. Antoine, and M. Lewenstein, Adv. At. Mol. Opt. Phys. 41, 83 (1997).
  * Ghimire  _et al._ (2011) S. Ghimire, A. D. DiChiara, E. Sistrunk, P. Agostini, L. F. DiMauro, and D. A. Reis, Nature physics 7, 138 (2011).
  * Chi  _et al._ (2018) S. Chi, Z. Li, Y. Xie, Y. Zhao, Z. Wang, L. Li, H. Yu, G. Wang, H. Weng, H. Zhang, _et al._ , Advanced Materials 30, 1801372 (2018).
  * Bedoya-Pinto  _et al._ (2020) A. Bedoya-Pinto, A. K. Pandeya, D. Liu, H. Deniz, K. Chang, H. Tan, H. Han, J. Jena, I. Kostanovskiy, and S. S. Parkin, ACS nano 14, 4405 (2020).
  * Jiang  _et al._ (2022) W. Jiang, G. S. Armstrong, J. Tong, Y. Xu, Z. Zuo, J. Qiang, P. Lu, D. D. Clarke, J. Benda, A. Fleischer, _et al._ , Nature Communications 13, 5072 (2022).
  * Eickhoff  _et al._ (2021) K. Eickhoff, L. Englert, T. Bayer, and M. Wollenhaupt, Frontiers in Physics 9, 675258 (2021).
  * Song  _et al._ (2023) X. Song, S. Yang, G. Wang, J. Lin, L. Wang, T. Meier, and W. Yang, Optics Express 31, 18862 (2023).
  * Li and He (2021) W. Li and F. He, Physical Review A 104, 063114 (2021).
  * Kerbstadt  _et al._ (2019) S. Kerbstadt, K. Eickhoff, T. Bayer, and M. Wollenhaupt, Nature Communications 10, 658 (2019).
  * Kumar  _et al._ (2017) N. Kumar, Y. Sun, N. Xu, K. Manna, M. Yao, V. Süss, I. Leermakers, O. Young, T. Förster, M. Schmidt, _et al._ , Nature Communications 8, 1642 (2017).
  * Wan  _et al._ (2011) X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Physical Review B 83, 205101 (2011).
  * Xu  _et al._ (2011) G. Xu, H. Weng, Z. Wang, X. Dai, and Z. Fang, Physical review letters 107, 186806 (2011).
  * Okamura  _et al._ (2020) Y. Okamura, S. Minami, Y. Kato, Y. Fujishiro, Y. Kaneko, J. Ikeda, J. Muramoto, R. Kaneko, K. Ueda, V. Kocsis, _et al._ , Nature communications 11, 1 (2020).
  * Bachmann  _et al._ (2017) M. D. Bachmann, N. Nair, F. Flicker, R. Ilan, T. Meng, N. J. Ghimire, E. D. Bauer, F. Ronning, J. G. Analytis, and P. J. Moll, Science Advances 3, e1602983 (2017).
  * Genet (2022) C. Genet, ACS photonics 9, 319 (2022).
  * Zhang  _et al._ (2017) T. Zhang, M. R. C. Mahdy, Y. Liu, J. H. Teng, C. T. Lim, Z. Wang, and C.-W. Qiu, ACS nano 11, 4292 (2017).
  * Rikken and Raupach (2000) G. Rikken and E. Raupach, Nature 405, 932 (2000).
  * Wang  _et al._ (2022) L. Wang, J. Lim, and L. J. Wong, Laser & Photonics Reviews 16, 2100279 (2022).
  * Hasan  _et al._ (2017) M. Z. Hasan, S.-Y. Xu, I. Belopolski, and S.-M. Huang, Annual Review of Condensed Matter Physics 8, 289 (2017).
  * Armitage  _et al._ (2018) N. Armitage, E. Mele, and A. Vishwanath, Reviews of Modern Physics 90, 015001 (2018).
  * Ishikawa (2010) K. L. Ishikawa, Physical Review B 82, 201402 (2010).
  * Lim  _et al._ (2020) J. Lim, Y. S. Ang, F. J. G. de Abajo, I. Kaminer, L. K. Ang, and L. J. Wong, Physical Review Research 2, 043252 (2020).