Keywords: Transportation mode detection/classification/identification, convolutional neural network, kick-scooters, micro-mobility, classification, inertial sensors

This work presents a novel transportation mode detection algorithm that handles the recognition of kick-scooters. In 2015, 10 minutes of data from a kick-scooter were considered in a transportation mode detection study, yielding a 56% F1-score. Since then, kick-scooters were not given much attention. Yet, kick-scooters are now very present in the urban transportation ecosystem, and their consideration in transportation studies has become a must. To fill this gap, 4 hours of kick-scooter signals were collected by 18 participants, with a set of 6 different kick-scooters, using 3 body-worn inertial measurement units. Obviously, kick-scooter patterns are classified in contrast with other modes of transportation. Two classification scenarios are considered in order to gradually increase the classification model complexity. The first scenario includes walking, biking, and kick-scooter, while the second considers public transport (tramway and bus) in addition to the former transportation modes. Results show that kick-scooters can be detected with an F1-score of 80% in the first scenario. Walking and public transport samples were still accurately classified in the second scenario, with an F1-score above 80% for both classes. However, bike and kick-scooter samples were both classified with lower F1-scores, equal to 59% and 64% respectively. Therefore, the main focus of future works should be directed towards the separability of kick-scooters and bikes when public transport is considered. The findings also suggest to place preferably the sensors in the trouser's pocket, allowing for leg motion to be finely captured.

Introduction

The growth of micromobility

In the last few years, light-weight vehicles have become popular among adults and young people. Bikes B. [START_REF] Wang | Short-term traffic flow prediction in bikesharing networks[END_REF], kick-scooters (KS), monowheels, and segways have nowadays become part of a larger transportation landscape designated as micro-mobility systems. KS have shown several advantages over bikes and other micro-vehicles. For example, KS are generally smaller, lighter, and can be folded and carried in tramways, buses, and metros. Therefore, they constitute relevant candidates to solve the first mile last mile (FMLM) issue, i.e. the distance that separates home or the place of work from the nearest public transport (PT) station, see [START_REF] Mckenzie | Urban mobility in the sharing economy: A spatiotemporal comparison of shared mobility services[END_REF]. In the same vein, the findings of [START_REF] Zuo | First-and-last mile solution via bicycling to improving transit accessibility and advancing transportation equity[END_REF] suggest that transit accessibility can be significantly increased by the use of bikes over walking. In fact, the faster the reachability of a station, the more it is likely to be accessed. The same conclusion was made in [START_REF] Shiv | Analysis of last mile transport pilot: Implementation of the model and its adaptation among local citizens[END_REF] where the lack of speedy mobility was found to be a catalyst for the use of private automobile to travel the FMLM distance. On the other hand, KS are green and health-friendly vehicles that are rather easy to handle. The PriestmanGood's kick-scooter for life1 is a typical example illustrating the ease of usability of KS by the elderly, which in other respects allows them to safely maintain physical activity. Furthermore, KS are allowed to share bike lanes and other public spaces, see [START_REF] Kostrzewska | Towards hybrid urban mobility: Kick scooter as a means of individual transport in the city[END_REF]. Thus, their integration to urban transportation systems is effortless in terms of infrastructure, see [START_REF] Oeschger | Micromobility and public transport integration: The current state of knowledge[END_REF]. As a result of the worldwide success of KS, it appears to us that their inclusion in Transportation Mode Detection (TMD) solutions is a requirement. Not only this will enhance urban navigation tools, but it will also improve urban policies and shared micro-mobility systems, see [START_REF] Ashqar | Network and station-level bike-sharing system prediction: A san francisco bay area case study[END_REF].

TMD families of methods

Sensors

Even though KS are new to TMD research, multiple TMD methods have been developed to recognize common transportation modes such as private cars, bikes, and PT. Recent TMD surveys can be found in [START_REF] Sadeghian | Review and evaluation of methods in transport mode detection based on gps tracking data[END_REF], [START_REF] Kamalian | A survey on local transport mode detection on the edge of the network[END_REF], [START_REF] Ahmed | Survey of machine learning methods applied to urban mobility[END_REF]. Historically, the main devices used were first GPS loggers [START_REF] Biancat | Review of transportation mode detection techniques[END_REF], [START_REF] Roy | Assessing the role of geographic context in transportation mode detection from gps data[END_REF] and then smartphones [START_REF] Carpineti | Custom dual transportation mode detection by smartphone devices exploiting sensor diversity[END_REF], [START_REF] Sharma | Early transportation mode detection using smartphone sensing data[END_REF], [START_REF] Wang | Combining residual and lstm recurrent networks for transportation mode detection using multimodal sensors integrated in smartphones[END_REF], [START_REF] Liu | Glmlp-trans: A transportation mode detection model using lightweight sensors integrated in smartphones[END_REF], which are still in use at the present, although a recent study has explored the use of smartwatches to perform TMD, see [START_REF] Hasan | Transportation mode detection by using smartphones and smartwatches with machine learning[END_REF]. The chosen technologies used for TMD purposes, regardless of the device, have evolved from the use of standalone GPS receivers [START_REF] Zheng | Learning transportation mode from raw gps data for geographic applications on the web[END_REF], often fused with Geographical Information Systems (GIS) [START_REF] Gong | A gps/gis method for travel mode in new york city[END_REF], [START_REF] Shah | Classifying the mode of transportation on mobile phones using gis information[END_REF], in favor of proprioceptive low-cost sensors due to their energy-efficiency and availability in off-the-shelf devices such as smartphones, see [START_REF] Yu | Big data small footprint: The design of a low-power classifier for detecting transportation modes[END_REF]. However, GPS receivers can still be used in favorable environments, typically in open-sky conditions, such as in [START_REF] Hadjidimitriou | Machine learning for activity pattern detection[END_REF], particularly when real-time functioning is not a requirement. Otherwise, the commonly used self-contained sensors are accelerometers, [START_REF] Nham | Predicting mode of transport from iphone accelerometer data[END_REF], gyroscopes, [START_REF] Jahangiri | Applying machine learning techniques to transportation mode recognition using mobile phone sensor data[END_REF], and magnetometers, [START_REF] Fang | Transportation modes classification using sensors on smartphones[END_REF]. Few works have used barometers [START_REF] Sankaran | Using mobile phone barometer for low-power transportation context detection[END_REF], microphones, L. Wang and Roggen (2019) and light sensors [START_REF] Carpineti | Custom dual transportation mode detection by smartphone devices exploiting sensor diversity[END_REF]. More exceptionally, bluetooth devices have been used in [START_REF] Coroamȃ | Exploring the usefulness of bluetooth and wifi proximity for transportation mode recognition[END_REF] to recognize transportation modes. Indoors, in a hospital-specific environment, Wireless Fidelity (WIFI) has also been used to classify hospital-specific modes of transportation, among which a KS, although the dataset comprised only 10 minutes of data from KS, see [START_REF] Prentow | Towards indoor transportation mode detection using mobile sensing[END_REF].

Time latency

On a methodological basis, there is a clear boundary in the literature between tripbased solutions as in [START_REF] Zhou | Travel mode detection method based on big smartphone global positioning system tracking data[END_REF], and near-real-time solutions, with a latency of recognition equal to the total duration of a trip in the former approaches, while it roughly varies between 1s to 10s in the latter, as in Wang et al. (2019), [START_REF] Prelipcean | Transportation mode detection-an in-depth review of applicability and reliability[END_REF]. Overall, the choice has been made in favor of low-cost and energy efficient sensors, together with low-latency solutions, lasting over 5s in average Wang et al. (2019), in order to meet the requirements of real-time applications such as navigation apps.

Classification models

Part of this work is model selection. Among the algorithms used in the literature, Random Forest (RF) is one of the most popular machine learning methods used in TMD problems, see [START_REF] Fang | Transportation modes classification using sensors on smartphones[END_REF]; [START_REF] Carpineti | Custom dual transportation mode detection by smartphone devices exploiting sensor diversity[END_REF]; [START_REF] Lorintiu | Transportation mode recognition based on smartphone embedded sensors for carbon footprint estimation[END_REF]; L. [START_REF] Wang | Summary of the sussexhuawei locomotion-transportation recognition challenge[END_REF]. More recently, deep learning methods were used for the same purposes, for example feed-forward artificial neural networks (ANN), see [START_REF] Fang | Learning transportation modes from smartphone sensors based on deep neural network[END_REF], convolutional neural networks (CNN), see [START_REF] Liang | A convolutional neural network for transportation mode detection based on smartphone platform[END_REF], and long-short term memory (LSTM) neural networks, see [START_REF] Asci | A novel input set for lstm-based transport mode detection[END_REF]. These models are generally compared, especially in recent studies such as Zarei Yazd, Taheri Sarteshnizi, Samimi, and Sarvi (2022). Let us recall though that the model selection here is performed based on two particular criteria. First, the chosen models are the most successful in the literature. Second, it is now well-established that the model selection is best realized based on cross-validation, see [START_REF] He | Automl: A survey of the state-of-the-art[END_REF][START_REF] He | Automl: A survey of the state-of-the-art[END_REF][START_REF] Zimmer | Auto-pytorch: multi-fidelity metalearning for efficient and robust autodl[END_REF][START_REF] Zimmer | Auto-pytorch: multi-fidelity metalearning for efficient and robust autodl[END_REF], which is also the approach followed in this study. Hence, RF, ANN, LSTM, and CNN were all tested systematically according to a 5-fold cross-validation process based on 1 separate and randomly selected test subjects for each fold. This way, the evaluation scenario adopted in this work guarantees enough robustness and is enhanced with respect to the major part of past TMD studies. In the latter, the train and test samples used to contain data from all subjects, see Taia Alaoui, Fourati, Kibangou, Robu, and Vuillerme (2022), resulting in an over-rated classification accuracy.

TMD datasets

As for publicly available datasets, sorted from most recent to oldest, the main datasets are the TMD-CAPTIMOVE, see Taia [START_REF] Taia Alaoui | Tmd-captimove[END_REF], which is the dataset used in this study, collected by 34 participants and having an amount of 48 hours of data, the SHL dataset, Wang et al. (2019), collected by 3 participants with a total duration of 703 hours, the US-TMD dataset, see [START_REF] Carpineti | Custom dual transportation mode detection by smartphone devices exploiting sensor diversity[END_REF], collected by 13 participants with a total duration of 32 hours, the dataset of [START_REF] Ballı | Diagnosis of transportation modes on mobile phone using logistic regression classification[END_REF], collected by 4 subjects with a total duration of 8 hours, and finally the oldest (2014), and at the mean time the richest in terms of sample size and amount of data, which is the HTC dataset, see [START_REF] Yu | Big data small footprint: The design of a low-power classifier for detecting transportation modes[END_REF], collected by 224 participants with a duration of 8311 hours. With respect to the other publicly available datasets, the current dataset (i.e. TMD-CAPTIMOVE) is optimized in terms of balance between the number of participants and the total duration, as shown in Figure 1. The x-axis provides the total duration of the considered TMD dataset in a log10 basis, the y-axis provides the number of subjects, also in a log10 basis, while the size of the circles provides the ratio between sample size (i.e. the number of participants) and the total duration, also in log10 for a better readability and visualization. Through this ratio, we wanted to recall the importance of optimizing the balance between the sample size and the total duration, as cross-subject variability, Taia [START_REF] Taia Alaoui | Urban transportation mode detection from inertial and barometric data in pedestrian mobility[END_REF], as well as timedependency inside the data, [START_REF] Prentow | Towards indoor transportation mode detection using mobile sensing[END_REF], were found to be critical in TMD training problems. In Figure 1, the most optimal dataset is the one on the top right corner and with the largest circle size, and the least optimal is the one at the bottom left corner and with the smallest circle size. Thus, the TMD-CAPTIMOVE is the most optimal with respect to circle size, which is the balance between the number of participants and the total duration. At the mean time, it is undeniably less optimal than the HTC dataset in terms of duration and sample size, and so it is as compared with the SHL dataset in terms of duration. More globally, it is an improved TMD dataset as compared with the other remaining public datasets. It is also more optimal and more balanced than many other private datasets, which have lower sample size and duration, Zhao, Hou, Alrobassy, and Zeng (2019), with 11 participants and 15 hours of data, Delli Priscoli, Giuseppi, and Lisi (2020), with 18 subjects and 104 hours of data, and Su, Yao, He, Lu, and Tong (2017), with 12 subjects and 8 hours of data.

On the other hand, it is the first dataset that provides hours of data (between 2 and 4 depending on the class) from KS and from elevators and stair-cases. Thus, it considers a seamless TMD transition from the indoor to the outdoor and vice-versa. Besides, the time amount of collected data from these modes of transportation is naturally lower than the duration of other TMD classes such as PT or walking. In fact, the latter last longer on a daily basis. Therefore, the fact of systematically balancing the amount of data for each class, [START_REF] Liang | A convolutional neural network for transportation mode detection based on smartphone platform[END_REF], most of the time at the cost of losing information, see [START_REF] Carpineti | Custom dual transportation mode detection by smartphone devices exploiting sensor diversity[END_REF], seems questionable.

KS as a novel class in TMD

In parallel to the technological and methodological challenges tackled by TMD research, the focus of this study is the augmentation of the types of transportation modes being classified. This augmentation is directly obained through the consideration of KS among other common transportation modes. Indeed, the most popular TMD classes are tramway, see [START_REF] Byon | Real-time transportation mode identification using artificial neural networks enhanced with mode availability layers: A case study in dubai[END_REF], bus, see [START_REF] Liang | A deep learning model for transportation mode detection based on smartphone sensing data[END_REF], [START_REF] Catalán | An automatic methodology to measure drivers' behavior in public transport[END_REF] and metro (or subways), see [START_REF] Fang | Learning transportation modes from smartphone sensors based on deep neural network[END_REF]. Bikes and e-bikes were considered in [START_REF] Xiao | Travel mode detection based on gps track data and bayesian networks[END_REF], whereas one study, [START_REF] Prentow | Towards indoor transportation mode detection using mobile sensing[END_REF], has considered 10 minutes of KS data, in a specific hospital environment, resulting into a KS-F1-Score of 56%. The latter was the first contributing study to KS integration to TMD methods. It is granted for having explicitly emphasized the significance of KS in TMD and for having provided precise information on their dataset collection scenarios and results.

More globally, due to the lack of studies on KS-dedicated TMD, this work will not make a thorough revision of the pros and cons of past TMD approaches. Instead, we investigate the possibility to extend existing TMD approaches to the use of KS. Therefore, the main issues addressed in this research are questioning the ability for KS to be accurately detected among other transportation modes. Another important question is whether the integration of KS impacts the accuracy at which other transportation modes are detected. These questions are answered based on experimental results obtained according to two scenarios. The first scenario consists in the identification of KS as a separate transportation mode from 2 other locomotion modes: on-foot and bike. The second scenario, which is more complex, adds to the first one signals from bus and tramway, merged together and considered as one class called PT.

Contributions

The main contributions of this work can be summarized as follows:

• For the first time, KS signal patterns are thoroughly studied. Two major KS modalities could be distinguished: pushing and cruising. This shows that the training of KS is more challenging as compared with mono-modal transportation modes such as walking or driving a car. • The placement of wearable sensors on the human body is studied as part of the variables to be optimized. Especially for KS, it is shown that the placement of sensors should not be random, which contrasts with past studies where TMD models were built regardless of the placement of the sensors on the body, see [START_REF] Carpineti | Custom dual transportation mode detection by smartphone devices exploiting sensor diversity[END_REF]. • We demonstrate that detecting KS with an F1-score above 80% is possible when considering on-foot activity and biking. However, the accuracy at which KS is detected decreases significantly when PT is considered. It is shown that on-foot and PT are not negatively impacted by the inclusion of KS. In contrast, bike and KS patterns become more challenging once the classification contains PT data. • The model validation is performed on data drawn from 5 randomly chosen test subjects (U1-U5), thus allowing us to evaluate the generalization ability of the model. Note that this aspect has been generally omitted in the TMD literature, see Taia [START_REF] Taia Alaoui | Urban transportation mode detection from inertial and barometric data in pedestrian mobility[END_REF].

Materials and Methods

Data collection

The data were recorded by 3 Inertial Measurements Units (IMUs) that include a tri-axis accelerometer, a tri-axis gyroscope and a barometer each. Their respective placements were on the foot, on the waist and in the trousers' pocket. They can be visualized in Figure 2 in orange. Each participant was equipped with a front camera for verification purposes, as shown in Figure . 3. The data used in this study were sub-sampled at a frequency of 32Hz. The data were collected by 34 participants from young to older ages, 20 men and 14 women aged from 18 to 50 years old, who volunteered to perform the experiments. They came from different social and professional backgrounds and were healthy subjects. There was no specific selection of participants based on either social, professional, or health factors. Such selection approach for the human subjects may affect the level of bias that potentially underpins the sample. The data collection was realized under different weather conditions from November 2019 to January 2020 during workdays at different day times, including rush hours (8-9am, 12-2pm, 5-6pm). According to our recruitment procedure, it is undeniable that the models built to perform the classification of TMD in this study are expected to be efficient basically for signals collected during workdays and rather by healthy and active subjects. Therefore, new data should be collected to adjust the algorithm for other populations, for example the elderly and individuals with varying levels of capability, as gait patterns are subject to variation in those cases, see [START_REF] Osoba | Balance and gait in the elderly: A contemporary review[END_REF].

In order to perform the experiments, the participants were equipped in the office. Then, they were accompanied during the whole experiment by the person in charge of later labelling the data. For tramways and buses, the itinerary was fixed at the beginning. For walking, biking, and KS riding, the volunteers were free to choose their itinerary. The experiments ended again in the office where the participants made a small report about the trip they made.

From 18 subjects, we were able to collect data using KS. In this scenario, they were asked to switch legs while riding the KS in order to balance pushing and resting phases for each foot. In total, 6 different KS were used. Each participant used one of the 6 KS to collect data. The other transportation modes involved at least 20 different participants and included walking, biking, and PT (tramway and bus). The time duration relative to each transportation mode is given in Figure 4. As usual in classification problems with real data, the classes are not equally distributed. The choice in this study was to take this imbalance into account at the validation stage by providing the F1-Score as an accuracy metric instead of only providing the overall accuracy. Another choice could be to balance classes by removing data, which the authors of this study think is not optimal as this inevitably results in information loss.

The dataset can be accessed at https://perscido.univ-grenoble-alpes.fr/ datasets/DS310.

Classified transportation modes

We adopt two different classifications of the transportation modes. They are given in Figures 5 and6. The first classification consists of 3 activity states that are: walking (W), biking (B), and KS riding. This scenario corresponds to body-induced motion. It is more suitable for short distance trips and for people who use their scooter or bike as a principal means of transport. The second classification integrates PT, which is the result of merging tramway and bus data. This scenario has more to do with the FMLM distance as KS are usually used to reach PT stations. Over all transportation modes, a systematic difference distinguishes the footmounted sensor from the waist-attached sensor and the one in the trouser's pocket. The signal magnitude is much higher in the first case. This difference is however less important for PT as the subjects generally stay still inside vehicles, which supposes that their feet are still as well. A slight magnitude difference is though observed between buses and tramways due to higher vibrations induced by the road topography as compared with rails for tramways (Figure 10). In this work, both transportation modes (tram and bus) are considered in the same class, PT.

For biking (Figure 8), the acceleration and angular rate signals show a clear timeperiodicity that reflects cycling. Yet, this cyclic pattern tends to fade in case the subject cruises for a period of time, for example between 4s and 7s. For KS, Figure 11 shows a signal frame that consists in two different KS patterns. The first pattern corresponds to pushing the ground with the foot on which the sensor is mounted. This motion is reflected by two consecutive signal peaks in both acceleration and angular rate. In fact, the first peak is the motion forward that allows the subject to gain momentum, while the second is the backward motion of the pushing foot before it hits the ground. The second pattern corresponds to cruising with both feet resting on the KS deck. These phases show flatter signals with small variation induced by ground friction while the scooter keeps moving by inertia. Signals were preprocessed prior to any segmentation or feature computation. The whole process of signal segmentation and preprocessing was assisted by videos that were time-stamped and synchronized with the IMUs. The participants were asked to behave freely during experiments. For example, they could move inside PT, walk inbetween scooter riding phases, etc. In order to take this freedom into consideration, a peak detection and deletion process was undertaken in order to remove sporadic motions that are not related to the transportation modes. For both accelerometer and gyroscope 3D components, Tukey's schematic boxplots, see [START_REF] Dawson | How significant is a boxplot outlier[END_REF], were calculated for multiple 1s long signal frames in order to detect outliers. The latter were then replaced by the upper and lower limits of the boxplot. Afterwards, the preprocessed 3D signals were used to compute the norm of acceleration and angular rate, which solves the problem of sensor orientation while later training the classifier.

In order to segment signals, we sought to minimize the segmenting window length while keeping satisfying classification accuracy, so that the detection time latency is minimized. To this purpose, many window length values, including 1s, 1.5s, 2s and 3s, were used to train 4 classifiers and validate them on a random test subject. The optimal window length value is that for which the classification accuracy is maximal. According to this approach, in the first classification scenario, this length was fixed to 1.5s. Note that it is much lower than the average sliding window length found in the state-of-the-art, which is about 5s, see Wang et al. (2019). For the second classification scenario, this value was set to 3s.

Classification Methods

Input data structure

In the training design, there are M train instances or train samples. In other words, M is the number of labeled signal patterns used for training. P is the number of signal 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 samples in one train instance, equal to the product of the segmenting window length and the signal sampling frequency (see Eq 1). Q is the number of signal channels considered for training. For example, if acceleration and angular rate signals are used, then Q is 2. Q is also the depth of the 3D matrix used to train the CNN and LSTM models described later in this section. S is the number of computed features based on the relevant signal channels (see section 2.6.2). The number of classes is N .

P = Segmenting W indow Length * Signal Sampling F requency (1)
In the following sections, the classification algorithms are provided in details. Note that there are two separate approaches used depending on the classification algorithm. For RF and ANN, the input data are 2D matrices of size MxS. In fact, these algorithms do not take into consideration the structure inside each train sample. Hence, computing descriptive features to train these 2 models results in a lower computational cost and higher training speed. For CNN and LSTM, the algorithms are intrinsically designed for structured data sequences. Therefore, no feature computation is performed with these 2 models. The input data for CNN and LSTM are time sequences, which means 1D vectors of length P . Yet, as Q signal channels are considered, an input train sample is a Q-dimensional vector of length P . Thus, the global input data of the LSTM and CNN models is a 3D matrix of size M * P * Q.

Random Forest

Random Forest (RF) is an ensemble of decision trees. The final prediction of a RF algorithm is computed through a majority voting based on the prediction of each decision tree. RF is known to prevent overfitting thanks to the bootstrap aggregation strategy that consists in a random selection of features and train samples to build each decision tree. The algorithm used in this study has 30 trees built upon a subset of 15 features each. The split criterion is the Gini index. The RF model is shown in Appendix A, in Figure 1. The Feed-forward Artificial Neural Network (ANN) used in this study has two hidden layers of 50 and 100 neurons respectively. Naturally, the output layer has as many neurons as trained classes, which is N. The activation function of the input and hidden layers is Relu. The output layer operates according to a softmax activation function.

Layer 2 and the output layer are optimized under an L2 regularisation with a lambda coefficient of 0.001. Each layer of the network is followed by a batch normalization operation to stabilize the parameters distribution through the training process. The optimizer is Adam with an initial learning rate of 0.01 and the cost function is the categorical cross-entropy. The ANN network architecture is shown in Appendix A, in Figure 2.

Long Short Term Memory cell

Long Short Term Memory cells (LSTM) are designed for time series in order to learn memory and sequential structure of data. The LSTM model used in this study has 1 LSTM cell, a hidden layer of 100 neurons under a Relu activation function and L2 regularizer, and finally, an output layer of N neurons under a softmax activation function. The loss function is the cross-entropy, and the optimizer is Adam with an initial learning rate of 0.001. The model architecture is shown in Figure 1 of Appendix B.

Convolutional Neural Network

The Convolutional Neural Network (CNN) is composed of 3 convolutional layers of 30 filters (kernels) each. The kernels size for each layer are respectively P/8, P/4, and P/4. A Maximum pooling of 2 was realized to the output of each convolutional layer. Two fully connected dense layers of 100 and 50 neurons follow the convolutional layers. The convolutional layers activation function was LeakyRelu with an Alpha value of 0.1. The activation function of the dense layers was ReLu. The output layer has N neurons and was designed with a softmax activation function. The cost function was the categorical cross-entropy. The optimizer was Adam, with an initial learning rate of 0.001. The global architecture of the CNN is shown in Figure 2 of Appendix B.

Computational tools and Transferability

RF was implemented based on Scikit learn library in Python. ANN, CNN, and LSTM were implemented using Keras library in Python as well. All computations were realized on an Intel Quad-Core i7 CPU of a laptop with a 16 Gb RAM memory. The operating system was Windows 10. All used libraries are free and can also operate under Linux and MacOS operating systems. As all the hyper-parameter values are provided, together with the models architecture, set of input features, and the online dataset. Therefore, the study is entirely reproducible.

Test scenarios and evaluation method

Test scenarios

Two main aspects are evaluated. First, we want to select the best suited features and signal inputs. Second, we want to select the most optimal model given the optimal input set. It is chosen among the selected algorithms from the literature: RF, ANN, LSTM, and CNN. For both purposes, cross-validation is used. It consists in training the models using 70% of the labeled samples, and then testing the stabilized models, i.e. after the training is done, to predict the labels (or classes) of the remaining 30% of the labeled samples. As commonly adopted in training problems, we performed a 5-fold cross-validation process. Indeed, the multiple-fold cross-validation process allows the evaluation of model stability against small variations in the training set. In this study, each fold contains data from all subjects except one, randomly chosen and labeled from U1 to U5. It is used to validate the model. A confusion matrix is provided at the end of each cross-validation. It provides the prediction results for all samples from the test set. This enables the computation of prediction errors and the computation of numerous accuracy metrics. Appendix C is provided with all confusion matrices. Details on the accuracy metrics used in this study can be found in section 2.6.4. Once these metrics are calculated, the choice has been made in this work to keep the model that leads to the best classification accuracy of KS.

Model selection is realized according to a two-stage scheme. First, we have trained the models using separate inputs, that is the selection of relevant input signal channels. Second, we have discarded the irrelevant signal channels, and trained the four algorithms using the optimal feature set, which allowed us to select the best model architecture. As the input selection phase is an intermediate step towards the final model selection, the methodology used to perform this selection is described in the section 2.6.2. Thus, the results in section 3 only present the performance of the four classification models used with the optimal set of inputs.

Selection of the optimal input set

Model selection requires first to select the relevant data to be used for training. Here, multiple signal channels were considered. They are provided in column 1 of Table 1. While most of the computed features are quite common, the energy might be questioned. It is defined in the time and frequency domains by equations 2 and 3 respectively.

Energy([t1, t2]) = t2 t1 s(t) 2 dt (2)
where t refers to time, and s to the considered signal. t1 is the beginning of the time interval while t2 is its end.

In the frequency domain, the same formula is used, using frequencies instead of time intervals:

Energy([f 1, f 2]) = f 2 f 1 s(f ) 2 df ( 3
)
where s is now the signal expressed in the frequency domain. f1 and f2 are respectively the minimum and maximum limits of the considered frequency interval.

For each signal channel, a classification was realized using LSTM and CNN algorithms. The signal channels with the lowest accuracy were removed while the others were taken together. The new set of signal channels was used with both algorithms. Then, each signal channel was iteratively removed. In case this led to a decrease in classification accuracy, the channel was fed back as an input to the neural network and another channel was removed to make a new test. In case removing the channel resulted in equal or higher classification accuracy, the channel was considered either useless or counter-productive, and therefore was discarded definitively. According to this analysis, the optimal set of signal channels was found to be composed of the gradients of acceleration and angular rate in both classification scenarios.

After the signal channels selection, descriptive time-domain and frequency-domain features were computed in order to be used with RF and ANN. They are provided in column 2 of Table 1. Note that the same set of descriptive features were computed for every channel. The optimal input set, given in the previous section, is used in the four classification models to be compared. This comparison is based on the chosen accuracy metric, which is the F1-score relative to KS. In other words, the best performing algorithm is the one that presents the highest accuracy for KS. For practical reasons, the confusion matrices, which contain all the details of the classification results, are provided only for the best model in Appendix C. However, the results in chapter 3 provide the global accuracy metrics for all classification models to visualize the performance obtained for each user. According to the test scenarios presented in section 2.6.1, the chosen model is the one that has the highest accuracy in average, taking into account the prediction results for 5 different test subjects (U1 to U5). As a reminder, for each classification scenario, there are 3 confusion matrices in columns (foot, waist, pocket) * 5 test subjects in row. For example, the first row of Appendix C consists of 3 confusion matrices obtained when testing the model on user 1 (U1), for the foot, the waist, and the trouser's pocket. The abbreviation for each class is: W: Walk, B: Bike, KS: Kick-scooter, PT: Public Transport. The statistics provided together with the confusion matrices are the OA (Eq. 5) as well as the geometric mean of the F1-score (Eq. 6). The imbalance of the dataset is taken into account by the F1score metrics. The geometric mean of the F1-score provides information on the global accuracy at which each class has been learnt. In other words, a high geometric mean of the F1-score indicates that all the classes were accurately learnt. On the contrary, a weak geometric mean of the F1-score indicates that at least 1 class has been weakly learnt, giving it the same weight as the other classes in spite of their unequal time distribution.

Accuracy Metrics

The accuracy metrics are chosen based on the questions that this study addresses. The first is to know how precise is KS detection among different urban transportation modes. The second is to know whether the integration of KS in TMD models makes the recognition of other transportation modes less accurate. Therefore, the classification accuracy results are given relatively to KS. Indeed, the first question is answered by the F1-Score, Eq (4), relative to KS. In order to answer the second question, global accuracy metrics are also provided. They correspond to the overall accuracy (OA), Eq (7), and the geometric mean of the F1-Score, Eq (8), provided in Tables 1 and2. The latter are global metrics that provide the classification accuracy of all classes taken together. In case the OA and geometric mean of the F1-Score are high enough (>70%), this means that all classes have been accurately detected. Otherwise, the F1-Scores of the classes other than KS are additionally analyzed in order to understand the errors of the classification model.

F1-Score = 2 * P recision * Recall P recision + Recall (4)
where the P recision is given by P recision = T rue P ositive T rue P ositive + F alse P ositive

(5)

and the Recall is given by Recall = T rue P ositive T rue P ositive + F alse N egative (6) OA = T rue P ositive + T rue N egative T rue P ositive + T rue N egative + F alse P ositive + F alse N egative (7)

Geometric-Mean(F1-Score) = N N k=1 F1-Score(k) (8)
3. Results

Model selection results

Model selection for scenario 1

In this section, several models are compared based on the F1-Score of KS for the 3 sensor placements. They are shown in Figure 12. The evaluation is made for 5 different test users (U1 to U5) that are not involved in training. From these figures, several conclusions can be drawn. They are summarized below:

(1) The highest F1-Score values are obtained using the CNN model with a footmounted sensor. In fact, the F1-Scores of all subjects are higher than 70%.

(2) The sensor placed in the pocket presents high F1-Score values for both the CNN (>70%) and LSTM (>65%) models. However, this performance varies significantly depending on the considered test subject.

(3) The F1-Scores relative to KS are the lowest for the waist-attached sensor. In fact, a high sparsity is observed over the different test subjects and different classification models. (4) ANN has very weak performance over all subjects and all sensor placements. In fact, this does not necessarily mean that ANN models are less accurate than the other models. However, this result could be interpreted as a higher need for the tuning of the ANN hyper-parameters as compared with LSTM and CNN.

Given this analysis, it seems that the CNN model has the best performance in the first classification scenario. For the second scenario, the F1-Scores of KS are shown in Figure 13. Here, the performance is globally much lower as compared with the first classification scenario. One common observation to the 3 sensor placements is that U5 has very weak F1-Score for all classification models. Therefore, this result cannot be explained by the choice of the training algorithms but seems to be rather related to the data. A possible hypothesis is that U5 has atypical signal patterns (ex. patterns that depend on physiological features). Therefore, as the data from U5 were removed from the training set, the algorithm did not learn these patterns. Another reason could be that U5 has collected more data than the other test subjects. When looking at the confusion matrices, this hypothesis has been confirmed. In fact, U5 has contributed to the collection of PT samples much more importantly than the other users, i.e. U1 to U4. Therefore, U5 seems to be a crucial subject to the training and should not be removed from the training dataset. As a consequence, the classification results from U5 should not be considered in model selection.

Removing U5 from our evaluation, considering only subjects from U1 to U4, the CNN has again the highest average F1-Score as compared with the other models. The same optimal sensor placements as for the first classification scenario are found here, on-foot, and in the trouser's pocket.

Performance of the best classification model: CNN

Scenario 1

After the selection of the relevant signal channels to be used for training, and after selecting the most performing algorithm in both classification scenarios, the present section provides more details about the achieved accuracies for KS specifically but also for the other transportation modes. To this purpose, Figure 14 shows the classification OA as well as the precise F1-Score values relative to KS for each sensor placement.

For the foot-mounted sensor, the F1-Score of KS varies from 78.02% to 84.03%. Its standard deviation is 2.61, meaning that the model is robust enough and can handle data from new subjects. The OA is comprised between 64.39% and 85.86%. In order to understand the lower OA value of 64.39% for U3, the confusion matrix for the foot-mounted sensor in Table 1 is analyzed. In fact, U3 has collected data only using a KS. Among these data, 2410 samples were classified as KS, 1324 were classified as biking, while 9 samples were classified as walking. In other words, the major errors here are due to confusing biking data samples with KS samples. More globally, the same tendency is observed for the other test subjects, however, it is most important for subject U3.

Surprisingly, even if the sensor placed in the trouser's pocket did not show strong signal patterns as compared with the foot-mounted one, it led to considerably high detection accuracy rates relative to KS, with an F1-Score comprised between 70.56% and 98.02% with a standard deviation of 11.52%. While U3 has an F1-Score relative to KS of 98.02% and an OA of 96.12%, the lowest OA is obtained for U5. Analyzing again the corresponding confusion matrix, it is observed that U5 has an OA of 60.64% and a geometric mean of the F1-Score of 61.39% for the sensor placed in the trouser's pocket. Looking more deeply into the confusion matrix, we can observe that around 1/3 of KS samples have been classified as biking.

Finally, the F1-Scores obtained with the waist-attached sensor are comprised between 7.69% and 83.13%. Without any further analysis of these results, the waist seems obviously not to be an optimal place to recognize KS. The results of the second scenario are shown in Figure 15. These results are given together with the detailed confusion matrices provided in Table 2. In this scenario, the global performance is lower with respect to the first scenario. As a reminder, U5 is not to be considered to interpret the classification results. However, the detailed statistics relative to U5 are still provided in order to understand the reasons why the accuracy is low for U5. More precisely, U5 made long travel periods of time by tramway and bus. Therefore, data from U5 were crucial to train the algorithm, which explains the low accuracy obtained when U5 is isolated as a test subject. This point can be clearly observed in Table 2 for U5.

Based on the results of the left test subjects, the F1-Score relative to KS ranges from 37.47% to 70.07% for the foot-mounted sensor, with a mean of 58.65% and a standard deviation of 15.6, which is an average performance that still requires postclassification improvements. For the sensor placed in the pocket, the F1-Score ranges from 33.86% to 86.24%. It has a mean of 59.35% and a standard deviation of 21.40%. The same conclusions as for the foot-mounted sensor are drawn. Finally, the waistattached sensor presents weak F1-Score values with an average F1-Score of 25.32 and a standard deviation of 21.68%.

The study of the F1-Scores of KS for the second scenario shows that KS is detected with a weak accuracy when PT is included to the classification, namely tramways and buses. A deeper analysis of the errors in scenario 2 is realized in the discussion section.

Classification model selection

In this work, 4 different algorithms were compared with the objective of providing a system that recognizes KS with enough accuracy. The tested algorithms were selected from the literature as mentioned in the introduction, with RF, ANN, CNN and LSTM. Selecting the relevant signal channels was crucial in this work. In this regard, the process of building several classifiers based only on one signal channel each has resulted in optimal performance. The validation of each algorithm was performed using a leave-one-subject out 5-fold cross-validation approach, meaning we have repeated the operation 5 times with 5 different test subjects. Therefore, the model validated in this study is robust enough to deal with data from unknown subjects. Note that we could isolate more subjects for testing, however, the dataset should be more significant in that case, especially as all subjects did not equally participate to data collection.

Based on these considerations, CNN was found to be the most optimal algorithm, with high classification accuracy rates and low standard deviations, which is an indicator of robustness.

Accuracy of KS detection

First classification scenario

The results show that it is possible to detect accurately KS out of walking and biking patterns using a sensor that is either foot-mounted or placed in the trouser's pocket, which is the answer to the first question of this study. For the foot-mounted sensor, the average F1-Score relative to KS is 79.36%, with a standard deviation of only 2.61%. For the sensor placed in the trouser's pocket, the F1-Score average relative to KS is 80.15% and its standard deviation is 11.52%. In this case, the increase in the F1-Score is insignificant, while the standard deviation is 5 times higher than the standard deviation of the foot-mounted sensor.

In order to address the second question of this study, which is about the impact of KS on the detection of other transportation modes, the geometric means of the F1-Score given in Table 1 and noted "F1-Score" for simplicity, are analyzed. For the most optimal sensor placement, which is on-foot, the average geometric mean of the F1-Score (which is averaged over all 5 subjects U1 to U5) is 82.56% and its standard deviation is 3.37. In fact, this is a significantly high value that means all transportation modes are accurately classified regardless of their distribution. Hence, introducing KS with walking and biking does not result in a weakened detection accuracy of these classes. As a final conclusion, these results suggest that KS can be detected accurately and therefore distinguished from walking and biking patterns, while keeping a highly satisfying detection accuracy of the latter.

Second classification scenario

In the second classification scenario, as PT is included to the classification, the F1-Score relative to KS decreases sharply. The mean F1-Score relative to KS is 58.65% for the foot-mounted sensor, 25.32% for the waist-attached sensor, and 59.35% for the sensor placed in the trouser's pocket. In other words, the detection accuracy of KS in the second scenario is rather low as compared with the first scenario.

The geometric mean of the F1-Score for the foot is 68.28% with a standard deviation of 20.56, while the geometric mean of the F1-Score for the pocket is 70.35% with a standard deviation of 14.11. In fact, the global classification accuracy is much higher than the F1-Scores relative to KS. In other words, transportation modes other than KS have been satisfyingly classified. In order to confirm this suggestion, the F1-Scores of all considered classes are analyzed. They are shown in Tables 2,3, and 4. From these tables, it is observed that PT samples were overall accurately classified, whilst the main classification errors consisted in classifying KS samples as biking, or biking samples as PT. In both cases, this is about confusing weak activity phases for biking and scooter riding with being inside PT. In fact, these are phases where the subjects are more likely to be cruising with both feet resting on the KS deck or on the bicycle, which may be confused with being inside PT with a low-energy profile. In fact, in both cases, acceleration and angular rate magnitudes are weak. Besides, it is important to note that the data were labelled as bicycling and KS riding even when the participants were in a static mode on their bike or KS, for example while waiting at a street light, which enhances these confusions. As a global conclusion about the detection performance in the second scenario, the introduction of PT results in a decreased detection accuracy of both biking and KS. In this case, it is not the introduction of KS that has a negative impact on other classes, but rather the introduction of PT that makes the classification of bike and KS less accurate.

Optimal sensor placement

This section recapitulates the classification results in both scenarios and for all sensor placements and all test users. The mean and standard deviation of the F1-Scores relative to KS are provided in Figure 16. As this is obviously observed in the figure, the foot is the placement that offers the best trade-off between F1-Score value and standard deviation. Note that the latter is supposed to be minimized while the F1-Score is sought to be maximized. Therefore, even if the pocket has a higher mean F1-Score relative to KS, it has much higher standard deviation, which reflects much higher values for the F1-Score for some test subjects, but also much lower values, which is to be avoided for real life applications. As a result, a foot-mounted sensor is recommended in studies that seek to optimize the robustness of KS detection.

Conclusion

This work was a first contribution to the detection of KS among other urban transportation modes with a new improved and publicly released TMD dataset. Our first focus was about the choice of the relevant signal channels, and classification algorithm, which came out to be the gradient of acceleration and angular rate used with an Adam-gradient descent CNN algorithm. In the first classification scenario, we were able to classify KS samples with a mean F1-Score of 80.18% for the sensor located in the trouser's pocket and 79.36% for the foot-mounted sensor. Moreover, these values were computed based on a 5-fold subject-based cross-validation method, which takes into account the model reliability with regard to handling data from various test subjects that were not involved in the training. Yet, once PT patterns were included, the detection accuracy decreased. In this case, an improving perspective could be to consider a pushing mode and a cruising mode for KS, as well as a cycling mode and a cruising mode for bikes, in order to decrease the intra-class variance, which is also a variable to be minimized while building robust classifiers. In summary, this study has been fulfilled in that it provided answers to the aforementioned questions. We proved that detecting KS is feasible, although it gets more challenging when PT is considered. Besides, results have shown that despite the difficulty to detect KS while including PT, PT were still accurately classified, meaning that KS inclusion in TMD does not result in a decrease of the detection accuracy of other transportation modes. An additional recommendation is to place the sensors in the trouser's pocket, or on the foot, in order to capture leg movements during KS riding. Finally, future contributions should also consider extending the sample distribution to include more varying signal patterns: healthy elderly, disabled, as well as signals collected during the weekends and holidays in order to study the time dependencies that may require necessary algorithmic adjustments.

Figure 1 .

 1 Figure 1.: Public TMD datasets: comparison of the balance between the number of participants and total duration

  Figure 2.: Position of the wearable sensors on the body: Foot-mounted (a), in the pocket (b), waist-attached (c)

Figure 4 .

 4 Figure 4.: Time duration relative to each transportation mode

  Figure 5.: First set of motion modes. (a): walk (b): bike (c): KS

  Figure 6.: Second set of motion modes. (a): walk (b): bike (c): KS (d): tramway (e): bus

Figure 7

 7 Figure 7.: Walking signal patterns: Acceleration norm in m/s 2 and angular rate norm in • /s

Figure 8

 8 Figure 8.: Biking signal patterns: Acceleration norm in m/s 2 and angular rate norm in • /s

Figure 9

 9 Figure 9.: Tramway signal patterns: Acceleration norm in m/s 2 and angular rate norm in • /s

Figure 10 Figure 11

 1011 Figure 10.: Bus signal patterns: Acceleration norm in m/s 2 and angular rate norm in • /s

  Figure 12.: Model comparison for scenario 1

  Figure 13.: Comparison of different classification algorithms in the second classification scenario

  Figure 14.: CNN classification results for the first scenario

  Figure 15.: CNN classification results for the second scenario

  Figure 16.: Classification results for the 3 sensor placements

Table 1 .

 1 : Set of computed signal channels and descriptive features

	Signal channels	Features computed for all channels
		Mean, Energy, Standard deviation, Variance,
	1) Acceleration norm, 2) Gradient of the norm of acceleration 3) Integral of the norm of acceleration 4) Angular rate norm 5) Gradient of the norm of angular rate 6) Integral of the norm of angular rate	Main frequency component, Power of the main frequency component, Spectral energy between 0 and 2Hz, Spectral energy ratio between 0 and 2Hz, Spectral energy between 2 and 4Hz, Spectral energy ratio between 2 and 4Hz, Spectral energy between 4 and 10Hz, Spectral energy ratio between 4 and 10Hz, Spectral energy between 10 and 30Hz, Spectral energy ratio between 10 and 30Hz, Spectral centroid, Spectral spread,
		Entropy, Number of zero crossings
	2.6.3. Selection of the best model

Table 2 .

 2 : F1-Scores of each class of the second scenario for the foot

	Foot	KS	PT	W	B
	Average F1-Score 58.65 91.59 89.67 67.21
	Std F1-Score	15.06 1.21	9.83	3.8
	Table 3.: F1-Scores of each class of the second scenario for the waist
	Waist	KS	PT	W	B
	Average F1-Score 25.32 83.02 78.97 46.76
	Std F1-Score	21.68 10.09 27.51 20.92
	Table 4.: F1-Scores of each class of the second scenario for the trouser's pocket
	Pocket	KS	PT	W	B
	Average F1-Score 59.35 85.26 81.4 64.37
	Std F1-Score	21.4	7.55 15.93 11.24

Table 2 .

 2 : Confusion matrices of the second scenario

	Foot	W	B PT KS	Waist W B PT KS	Pocket W B PT KS
		W 1339 0	0 12	W	626 0	0	6	W	930 0 0	40
	U1	B PT 397 326 4311 147 0 837 1 205	B PT	192 279 12 18 13 424 4242 8	B PT	344 0 158 821 23 0 4093 290
		KS 319 134 31 1024	KS	1189 672 71 87	KS	683 0 81 1328
	Overall Accuracy 82.45%			66.77%	72.24%
	F1-Score		77.63%			44.81%	53.71%
	Foot	W B PT KS	Waist W B PT KS	Pocket W B PT KS
		W 1481 0	6	1	W	3095 10 34 0	W	3021 5 16 0
	U2	B PT 27 100 4406 602 1 1506 16 1370	B PT	0 453 5 113 7 307 4718 2714	B PT	0 1375 16 121 167 234 47502214
		KS 42 81 29 1488	KS	289 900 173 699	KS	367 1748 63 1155
	Overall Accuracy 79.61%			66.32%	67.54%
	F1-Score		78.36%			58.67%	64.92%
	Foot	W	B	PT KS	Waist W B PT KS	Pocket W B PT KS
		W 0	0	0 4	W	0 0 0	5	W	0 0 0	0
	U3	B 0 PT 0	0 0	0 2837 0 39	B PT	0 0 0 1412 0 0 0 807	B PT	0 0 0 696 0 0 0 176
		KS 0	0	0 863	KS	0 0 0 1046	KS	0 0 0 2732
	Overall Accuracy 23.06%			31.99%	75.80%
	F1-Score		37.47%			48.47%	86.24%
	Foot	W B PT KS	Waist W B PT KS	Pocket W B PT KS
		W 1027 0	0 23	W	1946 0	0 12	W	1886 0 0 20
	U4	B PT 44 285 3801 259 1 1725 16 1566	B PT	26 1228 14 603 54 314 45031757	B PT	41 1273 40 481 40 336 3369 735
		KS 89 43	0 2174	KS	124 90 47 1674	KS	338 77 31 1220
	Overall Accuracy 78.96%			75.46%	78.36%
	F1-Score		79.65%			75.32%	76.54%

https://www.priestmangoode.com/project/scooter-for-life/

Acknowledgement(s)

The graphical content used in this paper was designed by Freepik and Flaticon.com.

Disclosure statement

The authors declare no conflict of interest.


L Turc

O W Roberts

M O Archer

M Palmroth

M Battarbee

T Brito

U Ganse

M Grandin

Y Pfau-Kempf

C P Escoubet

M W O Roberts

Battarbee

W Roberts

Palmroth

Brit

I Dandouras

First Observations of the Disruption of the Earth's Foreshock Wave Field During Magnetic Clouds

come

First Observations of the Disruption of the Earth's Foreshock Wave Field During Magnetic Clouds

Introduction

Magnetic clouds are strongly geoeffective solar transients, characterized by an enhanced and smoothly rotating magnetic field [START_REF] Huttunen | Properties and geoeffectiveness of magnetic clouds in the rising, maximum and early declining phases of solar cycle 23[END_REF][START_REF] Yermolaev | Geoeffectiveness and efficiency of CIR, sheath, and ICME in generation of magnetic storms[END_REF]. Understanding the details of their interaction with near-Earth space is crucial to forecast accurately their space weather effects. On their earthward journey, the first geophysical region that incoming magnetic clouds encounter is the Earth's foreshock.

The foreshock extends upstream of the quasi-parallel sector of the Earth's bow shock, where the 𝜃 Bn angle between the interplanetary magnetic field (IMF) and the shock normal is below ∼45 • (Eastwood, Lucek, et al., 2005). The foreshock is permeated with intense electromagnetic waves, generated through plasma instabilities triggered by shock-reflected particles (Eastwood, Lucek, et al., 2005;[START_REF] Wilson | Low frequency waves at and upstream of collisionless shocks[END_REF]. Foreshock processes can have global effects on the Earth's magnetosphere, causing enhanced wave activity down to the Earth's surface [START_REF] Bier | Investigating the IMF cone angle control of Pc3-4 pulsations observed on the ground[END_REF], or triggering geoeffective fast magnetosheath jets [START_REF] Plaschke | Jets downstream of collisionless shocks[END_REF]. Changes in the foreshock properties can therefore significantly affect conditions throughout geospace.

10.1029/2019GL084437

The most common waves in the Earth's ion foreshock are the so-called 30-s waves [START_REF] Eastwood | Cluster observations of fast magnetosonic waves in the terrestrial foreshock[END_REF], quasi-monochromatic magnetic field fluctuations at a period around 30 s, left handed in the spacecraft frame. Their wavelength is about 1 R E (1 Earth radius = 6,371 km), while their finite transverse extent ranges between 8 and 18 R E [START_REF] Archer | Size and shape of ULF waves in the terrestrial foreshock[END_REF]. Close to the wave vector direction, the wave front is essentially planar over several Earth radii, and the overall shape of the waves is approximately an oblate spheroid [START_REF] Archer | Size and shape of ULF waves in the terrestrial foreshock[END_REF]. They have been identified as fast magnetosonic waves propagating sunward in the plasma frame but advected earthward by the faster solar wind flow (Eastwood et al., 2005a). They are excited by backstreaming field-aligned beams (FABs) via the right-hand resonant ion-ion beam instability. The cyclotron resonance condition associated with this mode is

𝜔 = V beam k || -Ω ci (1)
where V beam is the beam velocity and Ω ci the ion gyrofrequency. These waves are generated around the mode interaction point of the cyclotron resonance with the fast magnetosonic branch, which can be approximated here as 𝜔 = v A k || , where v A is the Alfvén velocity (e.g., Eastwood et al., 2005a).

In spacecraft measurements, these waves are not observed in conjunction with the FABs that generated them, but with intermediate, gyrating or gyrophase-bunched distributions (Eastwood et al., 2005a;[START_REF] Kempf | Ion distributions in the Earth's foreshock: Hybrid-Vlasov simulation and THEMIS observations[END_REF]. Measurements taken shortly before or after observing the foreshock waves have however brought quantitative evidence of a cyclotron resonance between FABs and waves (Eastwood, Lucek, et al., 2005, and references therein). These waves have been extensively studied since their discovery [START_REF] Greenstadt | Correlated magnetic field and plasma observations of the Earth's bow shock[END_REF]Eastwood et al., 2005a[START_REF] Eastwood | Quasi-monochromatic ULF foreshock waves as observed by the four-spacecraft Cluster mission: 2. Oblique propagation[END_REF][START_REF] Palmroth | ULF foreshock under radial IMF: THEMIS observations and global kinetic simulation Vlasiator results compared[END_REF], and it is well established that their period depends on the IMF strength and orientation [START_REF] Le | Solar wind control of upstream wave frequency[END_REF][START_REF] Takahashi | Dependence of the spectrum of Pc 3-4 pulsations on the interplanetary magnetic field[END_REF].

Even though magnetic clouds are the most geoeffective solar wind disturbances, there are no studies focusing on the foreshock properties during such events. Recent numerical simulations predict that an enhanced IMF strength, as is encountered during magnetic clouds, could strongly affect the foreshock wave properties and their large-scale structuring [START_REF] Turc | Foreshock properties at typical and enhanced interplanetary magnetic field strengths: Results from hybrid-Vlasov simulations[END_REF]. These simulations were however limited to a single set of upstream conditions, warranting a more general observational investigation.

Event Identification

Observations of foreshock waves during magnetic clouds are rare, as these transients pass by Earth only about 2% of the time [START_REF] Yermolaev | Geoeffectiveness and efficiency of CIR, sheath, and ICME in generation of magnetic storms[END_REF], and Earth-orbiting spacecraft cross the foreshock only sporadically. In the early phase (2001)(2002)(2003)(2004)(2005) of the Cluster mission [START_REF] Escoubet | Introduction the Cluster mission[END_REF], the spacecraft separations are similar to the wave characteristic sizes (a few hundred kilometers), which allows us to determine accurately the wave properties. Using the catalogue introduced in [START_REF] Turc | Cone angle control of the interaction of magnetic clouds with the Earth's bow shock[END_REF], we identify events when Cluster observes foreshock waves during a magnetic cloud and divide them into 5-min intervals.

Foreshock fast magnetosonic waves are mostly transverse and propagate at a small angle relative to the magnetic field vector [START_REF] Eastwood | Quasi-monochromatic ULF foreshock waves as observed by the four-spacecraft Cluster mission: 2. Oblique propagation[END_REF]. Measurements from Cluster's fluxgate magnetometer [START_REF] Balogh | The Cluster magnetic field investigation[END_REF] are projected onto a frame where one axis is parallel to the mean magnetic field during the 5-min interval and the two others are perpendicular to it, forming a right-handed triplet. Then we perform a wavelet transform on the perpendicular magnetic field components and calculate the wavelet phase difference, to check the wave polarization [START_REF] Torrence | A practical guide to wavelet analysis[END_REF]. Intervals of steepened waves and discrete wave packets are visually identified and rejected.

In total, we find six magnetic clouds with observations of left-handed foreshock waves (see Table S1 in the supporting information). These observations take place several hours after the passage of the interplanetary shock preceding the magnetic cloud, when near-Earth space is embedded deep within the magnetic cloud, which dictates the upstream conditions for the foreshock to develop. These strongly differ from typical solar wind conditions, in particular, due to the clouds' large IMF strength.

Results

A More Intricate Form of Wave Activity

During most events, we find that the foreshock wave field departs significantly from its usual state. To facilitate the comparison, Figure 1 shows an example of typical foreshock waves on 18 February 2003 analyzed in previous studies [START_REF] Archer | Size and shape of ULF waves in the terrestrial foreshock[END_REF][START_REF] Hobara | Characteristics of terrestrial foreshock ULF waves: Cluster observations[END_REF][START_REF] Kis | Scattering of field-aligned beam ions upstream of Earth's bow shock[END_REF]left) together with representative 10.1029/2019GL084437 observations during a magnetic cloud on 19 January 2005 (right). In both cases, foreshock waves appear as large-amplitude magnetic field oscillations, especially in the two perpendicular components (panels c and d, and i and j). During quiet conditions, the wavelet power spectrum showcases a rather narrow band of strong wave power around 30 s (panel e), while during magnetic clouds (panel k) high-fluctuation power is observed at periods between 5 and 30 s. The contours highlight where the wave power is strongest, showing that the fluctuations are left handed in the spacecraft frame (in blue in the phase difference plot, panel l).

Cluster's relative position in the foreshock can be determined from its measurements prior to the foreshock waves observations. Before 13:34 UT on 19 January 2005, Cluster was just outside the Earth's foreshock. The IMF rotation inside the magnetic cloud results in a progressive change of the 𝜃 Bn angle along the field line connecting Cluster to the bow shock from 60 • at 13:30 UT to 40 • at 14:00 UT, according to a model bow shock [START_REF] Jeřáb | Improved bow shock model with dependence on the IMF strength[END_REF]. Consequently, the spacecraft probe the outer foreshock, populated by FABs, from 13:34 UT, and then the foreshock wave field from 13:42 UT. In such a configuration, one would expect to observe the usual quasi-monochromatic fast magnetosonic waves. However, the power spectra of the left-handed foreshock waves display multiple spectral peaks whose periods vary with time (Figure 1k). This intricate spectrum is neither associated with right-handed polarization (in red in the phase difference plot) nor with highly steepened waves.

To identify the wave mode, we use two independent methods, multispacecraft timing analysis [START_REF] Schwartz | Shock and discontinuity normals, Mach numbers, and related parameters[END_REF] applied to band-pass-filtered data to separate the different wave periods) and multipoint signal resonator (MSR; [START_REF] Narita | High-resolution wave number spectrum using multi-point measurements in space-The Multi-point Signal Resonator (MSR) technique[END_REF]; see supporting information). Their results for the first minutes of the data set are given in Table S2. The timing analysis yields wave vectors within 30 • of the magnetic field. With the MSR technique, we find that the power distribution in wave vector space maximizes at two different wave vectors, showing the coexistence of two waves. The orientation of all wave vectors k toward -x (i.e., earthward) and the negative wave velocities in the plasma frame V wave,pl indicate that the waves propagate sunward in the plasma frame and are intrinsically right handed. This rules out 10-s Alfvén waves, intrinsically left handed [START_REF] Eastwood | On the existence of Alfvén waves in the terrestrial foreshock[END_REF]. The Alfvén velocity was about 294 ± 158 km/s during this interval, the large uncertainties being due to large density fluctuations, as measured by the Cluster Ion Spectrometer [START_REF] Rème | First multispacecraft ion measurements in and near the Earth's magnetosphere with the identical Cluster ion spectrometry (CIS) experiment[END_REF]. The estimated wave speeds are comparable to the Alfvén velocity, which is close to 10.1029/2019GL084437 the fast magnetosonic speed for large IMF strengths. These properties are all consistent with those of fast magnetosonic waves.

Both wave analysis techniques therefore consistently identify a superposition of fast magnetosonic waves at different periods, ranging between 6 and 17 s. These periods, lower than the usual 30 s, are due to the magnetic cloud's large IMF strength, here 11 nT, while its average value at Earth is about 5 nT. The frequency of foreshock fast magnetosonic waves changes with the IMF strength because their dispersion relation depends on this parameter, among others [START_REF] Krauss-Varban | Mode properties of low-frequency waves: Kinetic theory versus Hall-MHD[END_REF]. According to the [START_REF] Takahashi | Dependence of the spectrum of Pc 3-4 pulsations on the interplanetary magnetic field[END_REF] empirical formula, the wave period should be 17 s under these solar wind conditions, very similar to the 18 s predicted by the [START_REF] Le | Solar wind control of upstream wave frequency[END_REF] formula, and consistent with the largest periods we observe. However, the lower wave periods are not accounted for by these models.

Similar intricate wavelet spectra are observed during four out of the six magnetic clouds under study (see the third column of Table S1). The waves were identified as fast magnetosonic waves during three of these events (22 January 2004, 19 January 2005, and 22 January 2005). For the last event, the wave mode could not be determined because of the extremely large IMF strength and solar wind velocity, resulting in the spacecraft tetrahedron being too large compared to the wavelength, making it impossible to identify uniquely the wave fronts.

The remaining two events display smoother wavelet spectra, as in normal solar wind conditions. The nature of the waves could not be confirmed during the 28 March 2001 event because of the poor correlation between the spacecraft time series, and of the large uncertainties in the flow velocity due to the very low density.

Both events with smoother spectra were associated with much larger IMF cone angles (about 45 • ) than the more complex spectra (less than 30 • ), where the cone angle is measured between the IMF vector and the Sun-Earth line. Larger cone angles result in larger wave periods for a given IMF strength [START_REF] Le | Solar wind control of upstream wave frequency[END_REF][START_REF] Takahashi | Dependence of the spectrum of Pc 3-4 pulsations on the interplanetary magnetic field[END_REF]), about 30 s on 23 April 2001. This suggests that the development of more complex foreshock wave activity is linked with shorter absolute wave periods.

A Possible Source for Multiple Fast Magnetosonic Waves

A recent numerical study also reported a superposition of fast magnetosonic waves in the foreshock associated with higher IMF strength [START_REF] Turc | Foreshock properties at typical and enhanced interplanetary magnetic field strengths: Results from hybrid-Vlasov simulations[END_REF]. In the simulation, the multiple wave periods were attributed to ion velocity distribution functions (VDFs) with multiple FABs, instead of the single beam that usually generates quasi-monochromatic waves.

During all intervals under study, the ion VDFs have already evolved into intermediate and gyrating distributions when foreshock waves are observed, in agreement with previous works (e.g., Eastwood et al., 2005a). Cluster, however, probed the FAB region shortly before most of these intervals. We focus here only on those intervals for which fast magnetosonic waves were reliably identified (23 April 2001, 22 January 2004, 19 January 2005, and 22 January 2005). We checked that the upstream conditions remained essentially steady from the FAB observations to that of the waves, except for the IMF direction, which causes the foreshock wave field to reach the spacecraft. Based on the current understanding of foreshock wave generation, we can reasonably assume that the FABs we observe in the first place are representative of those generating the waves detected a few minutes later.

The upper part of Figure 2 displays representative VDFs recorded by the Hot Ion Analyzer instrument onboard Cluster-1 at 13:41:10 and 13:41:18 UT on 19 January 2005. The left column shows reduced two-dimensional (2-D) VDFs in the (V || , V ⟂ ) plane (relative to the magnetic field vector), integrated over the second perpendicular direction. In panels a and b, part of the solar wind population appears in red in the right-hand side, while backstreaming ions are in the left-hand side. The suprathermal population consists mostly of a FAB, centered at V || ∼ -1, 500 km/s. There is no clear evidence of multiple distinct FABs during this event, but a second population is observed at a nonzero pitch angle, centered around V ⟂ ≃ 1, 000 km/s. Panels c and d show the phase space density of the suprathermal population along a cut at V ⟂ = 0. The beam shape is well fitted by a Gaussian (see also Table S3), most likely because the gyrating population is located at about the same V || as the FAB. This second population appears as a second peak in phase space density on the cuts along V ⟂ at V || = V beam (panels e and f) and remain clearly visible for a range of V || (not shown). The suprathermal population is better fitted by a sum of two Gaussians in these cuts (see Table S3), thus confirming that the two peaks are well distinct.

10.1029/2019GL084437

During the 23 April 2001 event, typical quasi-monochromatic foreshock waves were observed, accompanied with typical FABs. During the 22 January 2004 event, a sharp rotation of the IMF shortly before 12:00 UT causes a rapid motion of the foreshock boundary. Only two to three VDFs are recorded when Cluster crosses the FAB region. Given the rapid motion of the boundary, we cannot draw firm conclusions regarding the FABs during this event.

Finally, the 22 January 2005 event brings clear evidence of multiple FABs (Figure 2, bottom part) similar to those reported by [START_REF] Turc | Foreshock properties at typical and enhanced interplanetary magnetic field strengths: Results from hybrid-Vlasov simulations[END_REF]. Because of the different IMF orientation, the solar wind population in panels g and h is at negative V || during this event. The FABs are centered around V || ∼ 700 and 1,700 km/s, respectively, and are well fitted by a sum of two Gaussians (panels i and j and Table S3). From these fits and Gaussian fits along V ⟂ at V || = V beam1 (panels k and l) and V || = V beam2 , we estimate the temperature anisotropy (ratio of the perpendicular to parallel temperatures) of each of the beams (see Table S4). They range between 2 and 10, in excellent agreement with previous observations [START_REF] Paschmann | Characteristics of reflected and diffuse ions upstream from the Earth's bow shock[END_REF]. The two FABs evolve progressively in consecutive VDFs for about a minute. At 01:24:00 UT, the higher-energy beam reaches its highest phase space density, peaking at about 4.1 × 10 -9 s 2 /m 5 , only a factor of 3 lower than the other beam, at about 1.3 × 10 -8 s 2 /m 5 .

To our knowledge, this is the first time that multiple FABs are observed in conjunction with unusual foreshock wave activity. Spacecraft observations of multiple FABs have been reported by [START_REF] Meziane | Anomalous foreshock field-aligned beams observed by Cluster[END_REF] but were associated with typical quasi-monochromatic waves. The IMF strength was relatively high (9 nT) during the main interval analyzed in [START_REF] Meziane | Anomalous foreshock field-aligned beams observed by Cluster[END_REF] but was not associated with a magnetic cloud.

The comparison of the cyclotron resonant speed of the waves with that of the FABs revealed that only the main beam was in cyclotron resonance with the waves [START_REF] Meziane | Anomalous foreshock field-aligned beams observed by Cluster[END_REF].

The cyclotron resonant speed of the waves, normalized to the solar wind speed, is calculated as

P res = Ω ci cos 𝜃 kV 𝜔 sc cos 𝜃 kB (2)
where 𝜔 sc is the wave frequency in the spacecraft frame and 𝜃 kV (𝜃 kB ) is the angle between the wave vector and the solar wind velocity vector (the IMF vector; [START_REF] Meziane | Anomalous foreshock field-aligned beams observed by Cluster[END_REF]. We calculate P res for the multiple FABs observed around 01:24 UT on 22 January 2005. We find that the velocity of the lower-energy beam normalized with the solar wind speed ranges between 1.8 and 1.9 and that of the higher-energy beam between 3.0 and 3.1. This is in excellent agreement with the cyclotron resonant speeds of the waves observed from 01:32 UT onward, which are 1.9 and 3.3.

Structuring of the Foreshock at Smaller Scales

During magnetic clouds, fast magnetosonic waves have shorter periods, which, in the plasma rest frame, translates into shorter wavelengths. Another critical parameter for the structuring of the foreshock wave field is the wave correlation length. This parameter can be estimated using cross correlations of measurements from pairs of spacecraft (see [START_REF] Archer | Size and shape of ULF waves in the terrestrial foreshock[END_REF], and supporting information). Here, we calculate the wave correlation length in the plane perpendicular to the wave vector. We divide all our intervals of fast magnetosonic waves into 120-s sections, so that for each interval the wave period and the upstream parameters remain roughly constant while retaining a sufficient number of wave periods. We obtain a reliable estimate of the transverse extent of the waves for 35 intervals, displayed as a histogram in Figure 3a. The transverse extent ranges between 1 and 10 R E (average 3.5 R E ; median 3.2 R E ). This is significantly shorter than during quiet solar wind conditions (8-18 R E ; [START_REF] Archer | Size and shape of ULF waves in the terrestrial foreshock[END_REF]. Moreover, we find that the transverse extent of the waves is well correlated with the wave period, as evidenced by Figure 3b (correlation coefficient: 0.83), and we checked that this result is insensitive to the chosen interval length. This correlation implies that the waves retain the same aspect ratio when their wavelength varies; that is, the ratio of their wavelength over their transverse extent is roughly constant in the events under study.

To get a clearer view of the large-scale structuring of the foreshock, we use two global simulations performed with the hybrid-Vlasov Vlasiator code [START_REF] Palmroth | Vlasov methods in space physics and astrophysics[END_REF][START_REF] Von Alfthan | Vlasiator: First global hybrid-Vlasov simulations of Earth's foreshock and magnetosheath[END_REF], already studied in [START_REF] Turc | Foreshock properties at typical and enhanced interplanetary magnetic field strengths: Results from hybrid-Vlasov simulations[END_REF]. In both runs, the simulation domain is 2-D in real space, describing the equatorial plane of near-Earth space. The IMF has a 5 • cone angle, the solar wind velocity is V SW = (-600, 0, 0) km/s and the ion density n SW = 3.3cm These results are in reasonable agreement with Cluster data, supporting the fact that the simulated foreshock wave field is representative of the actual foreshock. This suggests that during quiet solar wind conditions, the foreshock wave field is composed of large-scale coherent waves. When moving closer to the bow shock, their transverse extent decreases. At large IMF strength, as is the case during magnetic clouds, the correlation length of the waves becomes significantly shorter.

Discussion and Conclusions

We have presented the first observations of the Earth's foreshock during magnetic clouds, revealing that the foreshock develops during such events but its wave properties are strongly modified due to the unusual upstream conditions dictated by the clouds. Using multispacecraft analysis techniques, we show that the usually quasi-monochromatic fast magnetosonic waves are replaced by a superposition of waves at different periods. Their wavelength and their transverse extent are shorter, suggesting that the foreshock is structured over smaller scales.

Atypical ion VDFs are observed in conjunction with the unusual wave activity. During one event, we find clear evidence of two distinct FABs in cyclotron resonance with the waves observed shortly afterward, consistent with previous numerical works [START_REF] Turc | Foreshock properties at typical and enhanced interplanetary magnetic field strengths: Results from hybrid-Vlasov simulations[END_REF]. In another, the FAB is accompanied by a second suprathermal population at nonzero pitch angle, which might have evolved from a second FAB, due to gyrophase trapping [START_REF] Mazelle | Production of gyrating ions from nonlinear wave-particle interaction upstream from the Earth's bow shock: A case study from Cluster-CIS[END_REF][START_REF] Kempf | Ion distributions in the Earth's foreshock: Hybrid-Vlasov simulation and THEMIS observations[END_REF]. The apparent lack of multiple FABs during this event could also be due to the spacecraft connectivity to the bow shock. [START_REF] Meziane | Anomalous foreshock field-aligned beams observed by Cluster[END_REF], since their energy depends on the local shock geometry upon their generation [START_REF] Paschmann | Energization of solar wind ions by reflection from the Earth's bow shock[END_REF]. Foreshock processes are a significant source of local shock deformations [START_REF] Meziane | Anomalous foreshock field-aligned beams observed by Cluster[END_REF]. Therefore, the observations of multiple FABs being more likely at lower 𝜃 Bn values, that is, closer to the foreshock, fits well within this scenario.

We note here that the second beam in the [START_REF] Meziane | Anomalous foreshock field-aligned beams observed by Cluster[END_REF] event was not associated with intricate wave activity. Its maximum phase space density was about 2 orders of magnitude lower than that of the first beam, which could explain why it did not trigger additional fast magnetosonic waves. On the contrary, on 22 January 2005, the maximum phase space densities of the two beams only differ by a factor of 3, and both beams can thus generate fast magnetosonic waves of comparable amplitudes, as their growth rate increases with increasing beam density [START_REF] Gary | Electromagnetic ion/ion instabilities and their consequences in space plasmas-A review[END_REF].

In this work, the wave properties have been determined assuming a linear picture, which remains appropriate here as the spectral peaks of the waves are well distinct. This allows separation of the different wave modes, as done previously in [START_REF] Hobara | Characteristics of terrestrial foreshock ULF waves: Cluster observations[END_REF]. Nonlinear wave analysis may be necessitated in other circumstances where distinct periodicities are not present.

Foreshock waves are known to modulate the shape of the shock front [START_REF] Burgess | Foreshock-shock interaction at collisionless quasi-parallel shocks[END_REF]. Therefore, their smaller wavelength and their smaller transverse extent during magnetic clouds could both result in smaller 10.1029/2019GL084437 ripples at the shock front, which in turn can affect particle reflection at the quasi-parallel bow shock [START_REF] Wu | The role of large amplitude upstream low-frequency waves in the generation of superthermal ions at a quasi-parallel collisionless shock: Cluster observations[END_REF] and the formation of magnetosheath high-speed jets [START_REF] Plaschke | Jets downstream of collisionless shocks[END_REF].

IMF strengths above 10 nT are relatively uncommon at Earth, but they become typical closer to the Sun. At Mercury's orbit, the average IMF strength is about 20-30 nT [START_REF] Korth | The interplanetary magnetic field environment at Mercury's orbit[END_REF]. The small size of Mercury's magnetosphere results however in another organization of its foreshock wave field [START_REF] Le | Upstream ultra-low frequency waves in Mercury's foreshock region: MESSENGER magnetic field observations[END_REF]. Outside of our solar system, exoplanets orbiting close to their host stars are immersed in intense magnetic fields and could thus display similar foreshock properties as presented here.

Figure 1 .

 1 Figure 1. Foreshock waves on 18 February 2003 (left), during quiet solar wind conditions, and on 19 January 2005 (right), during a magnetic cloud.From top to bottom: magnetic field strength and components parallel and perpendicular to the mean magnetic field vector, wavelet power spectrum of B ⟂ 1 , and wavelet phase difference betweenB ⟂ 1 and B ⟂ 2 . The dotted lines in panels e and f show the 30-s period. The black (white) contours in panels e and k (f and l) delineate the 95% significance level of the wavelet power spectrum. The hatched areas in panels e and f, and k and l show the cone of influence, where edge effects become significant in the wavelet spectrum.

Figure 2 .

 2 Figure 2. Ion VDFs observed shortly before the foreshock waves. (a-f) The 19 January 2005 event. (g-l) The 22 January 2005 event. Left column: reduced 2-D VDFs, integrated over the second perpendicular velocity direction. Color coded is the phase space density. The x and y axes are along and perpendicular to the magnetic field direction, respectively. Middle and right columns: cuts through the reduced distribution functions shown in the left column, at V ⟂ = 0 (middle) and V || = V beam (right). VDF = velocity distribution function.

Figure 3 .

 3 Figure 3. Transverse extent of the foreshock waves. Upper part: Cluster's observations during magnetic clouds. (a) Distribution of the transverse extent of foreshock waves. (b) Wave period as a function of the transverse extent. Bottom part: global simulations. (c and d) Magnetic field B z component in Runs 1 and 2 at time t = 500 s from the beginning of the runs, which illustrate the foreshock wave field. The plus signs indicate the positions of triplets of virtual spacecraft where the transverse extent of the wave fronts was reliably determined. (e) Distribution of the transverse extent of the wave fronts.

  -3 . The IMF strength is set to 5 nT in Run 1, corresponding to regular IMF strength at Earth, and 10 nT in Run 2, comparable to the 19 January 2005 event.

  out-of-plane (B z ) component of the magnetic field in the simulation domain is displayed in Figures 3c (Run 1) and 3d (Run 2). For clarity, the regions downstream of the bow shock are not shown. Because the foreshock waves are mostly transverse, B z is a good marker of the foreshock wave field. As can readily be seen from these plots, the waves are coherent over much larger scales during quiet solar wind conditions. The plus signs indicate the barycentric positions of triplets of virtual spacecraft, mimicking the Cluster constellation in 2-D, where the transverse extent of the waves was estimated reliably. Since the simulation is 2-D, only three spacecraft are needed to determine the wave properties. Magnetic field time series are extracted at each virtual spacecraft location and we estimate the transverse extent of the waves using the same cross-correlation technique as in Cluster data, albeit only with three spacecraft. The results of the analysis are shown as a histogram in panel e. The transverse extent of the waves varies significantly depending on the position in the foreshock. We obtain values ranging between 3 and 11 R E in Run 1 and between 2 and 6 R E in Run 2, the smallest values being encountered closer to the bow shock, especially in Run 1. The average (median) transverse extent of the waves is 7 (7.4) R E in Run 1 and 4 (4.3) R E in Run 2.

	Geophysical Research Letters	10.1029/2019GL084437

The

  Using a model bow shock, we estimate that the observed FABs originate from 𝜃 Bn ∼ 52 • during this event, and from 𝜃 Bn ∼ 43 • when multiple FABs are observed on 19 January 2005 and in the Meziane et al. (2011) event. Local changes of the quasi-perpendicular shock geometry could generate multiple FABs

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement 704681. We

 and312390). We acknowledge the Academy of Finland (grant numbers 322544 and 309937) for financial support. We thank the Cluster Science Archive (CSA; [START_REF] Laakso | Cluster active archive: Overview[END_REF] and the Cluster Ion Spectrometer and FGM PI teams for providing Cluster data. The CSC-IT Center for Science in Finland is acknowledged for the Pilot run and the Grand Challenge award leading to the results shown here. We thank S. von Alfthan for his major contribution in the development of Vlasiator. The left panels of Figure 2 were done with the QSAS science analysis system provided by the United Kingdom Cluster Science Centre (Imperial College London and Queen Mary, University of London) supported by The Science and Technology Facilities Council (STFC). Data from the Cluster mission is freely available on the CSA. The Vlasiator runs described here take several terabytes of disk space and are kept in storage maintained within the CSC-IT Center for Science. Data presented in this paper can be accessed by following the data policy on the Vlasiator website.