Optimizing the topology of transmission networks using drl has increasingly come into focus. Various drl agents have been proposed, which are mostly benchmarked on the Grid2Op environment from the l2rpn challenges. The environments have many advantages with their realistic chronics and underlying power flow backends. However, the interpretation of agent survival or failure is not always clear, as there are a variety of potential causes.

In this work, we focus on the failures of the power grid simulation to identify patterns and detect them in advance. We collect the failed chronics of three different agents on the WCCI 2022 l2rpn environment, totaling about 40k data points. By clustering, we are able to detect five distinct clusters, identifying common failure types. Further, we propose a multi-class prediction approach to detect failures beforehand and evaluate five different prediction models. Here, the lgbm shows the best failure detection performance, with an accuracy of 86%. It also correctly identifies in 91% of the time failure and survival observations. Finally, we provide a detailed feature importance analysis that identifies critical features and regions in the grid.

§ INTRODUCTION

With the increase of renewable energy in the energy mix and the resulting variability, the optimization of power grids has become increasingly complex. One possible approach is topology optimization, i.e., switching buses at a substation level, it is a more cost-effective alternative to redispatching [24] and can improve grid stability [4]. However, identifying the correct topology options is computationally expensive [24]. One possible solution for the topology optimization could be drl, as it can grasp large states and cope with the complexity [29].

Within the ml research community, the l2rpn challenges are the preferred research benchmarks [17, 22, 23, 24, 28]. The underlying Grid2Op package [10] of the tso RTE provides several environments for drl approaches.[Grid2Op: <https://grid2op.readthedocs.io/en/latest/> (last accessed 09/06/2024).]
Within the environments, researchers train their agents on various chronics to identify policies. To rate their agents, a scoring function is provided [23], which is highly dependent on the survival of the agent in the provided test chronics. However, it is not always clear why an agent survives or fails in a specific chronic, as the environments contain stochastic components that change the behavior on repeated runs.

Several topology-based drl agents have been applied to the l2rpn challenges of RTE in the last years [24]. However, as drl agents are evaluated on their survival in the grid, it is not always clear why they survived or failed and caused a grid failure. Therefore, we analyze different grid failures in the IEEE118 grid and provide a predictive approach to detect possible failures in advance.

§ CONTRIBUTION

We take a closer look at agent failures to address the types of failures that can occur when using a topology-based agent. To do so, we aggregate the data of over 40k failed chronics from three different agents and cluster the failure types to identify specific patterns.

Next to clustering, we provide a multi-class prediction framework to identify whether an agent is currently in a situation where survival can be ensured, or if failure is imminent. Based on the observation of the grid state, we can detect whether the grid will survive or if a failure will occur in 5, 3, or 1 time-steps. We test different prediction models and are able to identify with 91% binary accuracy, whether grid failure is imminent. With a detailed feature importance analysis, we are further able to detect critical features and components in the electricity grid that are responsible for the instability of the power grid.
Overall, the contributions can be summarized as follows:

* We evaluate three agents from Ref. [19] with ten different seeds on the IEEE118 grid (WCCI l2rpn training environment). With 1662 chronics, we create a detailed dataset of around 40k grid failures.

* We perform a clustering analysis that identifies distinct grid failure types and categorize them into five meaningful clusters, enhancing understanding of agents' shortcomings and their association with specific failure types.

* We build a multi-class forecasting framework to identify the risk of failure beforehand.

* We test five different forecasting models and identify the lgbm as the best performing candidate.

* With an additional feature importance analysis, we are able to identify critical features and regions within the IEEE118 grid.

The rest of the paper is structured as follows. We provide an overview of the related work in Section <ref> and our methodology in Section <ref>. We then present the experimental setup in Section <ref>, followed by the results in Section <ref>. Finally, we discuss our results in Section <ref> and conclude in Section <ref>.

§ RELATED WORK

With the l2rpn challenges, several approaches have been proposed to optimize the topology of transmission grids with drl. In most cases, the drl algorithms are combined with heuristic components [5, 12, 19, 20] to increase performance and reduce the action space. Other researchers prefer to use evolutionary algorithms with planning components [31] or mcts [11]. However, most researchers report only the score and in some cases the survival of their agent on the test chronics. Understandably, they often do not show the shortcomings of their agents, making comparisons or even fault interpretation difficult.
The most detailed analysis of the environments and agents comes from the l2rpn challenges themselves, e.g., in [23, 25]. They analyze the behavior of different agents and also show some examples using the Grid2Viz package.[Grid2Viz <https://github.com/rte-france/grid2viz> (last accessed 09/06/2024).] Alternatively, [3] introduce Grid2Onto for Grid2Op environments. Their application is a recommender system based on a semantic model to display and select propositions from an agent for the respective grids.

With respect to fault prediction, previous studies have investigated electricity grids, as evidenced by the surveys for transmission grids [6] and distribution grids [8]. However, to the best of our knowledge, there is no specific work on fault prediction using drl agents for topology optimization. The closest work was done by [26], who used an LSTM model to predict failures on the IEEE 35 bus network. However, they focus more on interlocking and cascading failures and not directly on the outcome of drl agents. Further, their study is relatively small, with a training sample size of 1440 and 2048 generated samples. Therefore, there is still a need for further research in this area.

§ METHODOLOGY

§.§ Descriptive Analysis with Clustering

In the descriptive analysis, we want to identify different types of grid failures to better understand the shortcomings of the agents. However, the underlying data is high-dimensional and requires reduction.[The aggregation of the data is described in Section <ref>.] To ensure that we have less correlation in our clustering, the next step is to use the pca [1] with standardized variables. In the pca, we reduce the dimensionality by transforming the data onto the coordinate system of the principal components. As the principal components are orthogonal to each other, they are considered uncorrelated. We then select all major components that added at least 5% of the variance to the data, thus reducing the dimensions effectively. In the case of our data, this yields six components with an explained variance of 85%.

With a reduced number of components, we combine the pca with k-means clustering, similar to [9]. K-means clustering was first introduced by [21] and divides the data into k-clusters.
Within the clustering, each data point is assigned to the nearest cluster centroid based on its distance. By iteratively assigning precise centroids, the distances within the cluster is minimized. Since the cluster size k has to be chosen, we decided to iterate over different sizes and choose the cluster size based on the inertia and the silhouette metric [18].

§.§ Forecasting Models

In our forecasting framework, we go beyond a binary prediction of failure versus non-failure of the grid. Instead, we adopt a multi-class approach, forecasting potential failures at three distinct time steps before they occur, namely five, three and one time-steps ahead. This allows for a more granular understanding of the impending failures and provide insights into the decision-making of the models such as the average probability distribution among all classes. We employ the following five prediction models and compare their performance. In particular, we utilize gradient boosting techniques for their strong performance across various prediction tasks and their ability to easily provide insights into feature importance.

rf 13 is an ensemble learning method that constructs simple decision trees during training and outputs the majority vote of their predictions for classification tasks. This approach enhances model accuracy and robustness by reducing overfitting and improving generalization compared to individual decision trees.

xg 7 uses gradient boosting to iteratively build decision trees, such that each new tree is trained to correct the residual errors of all previous trees. It employs advanced techniques like regularization, weighted quantile sketch, and sparsity awareness to improve model accuracy, prevent overfitting, and handle large datasets efficiently.

lgbm 16 is an ensemble learning method for classification that builds decision trees using a gradient boosting framework, optimizing for speed and memory efficiency. Unlike xg, lgbm uses a histogram-based algorithm and leaf-wise tree growth, which significantly enhances training speed and reduces memory usage, especially for large datasets.

cem 14 is a basic yet effective feed-forward neural network tailored for tabular datasets. Its architecture integrates categorical features through a learnable embedding layer, making it a suitable starting point and baseline for comparison with other models.

GANDALF 15 stands for Gated Adaptive Network for Deep Automated Learning of Features and optimizes feature selection and engineering using gflu. gflu leverages a gating mechanism inspired by gru to iteratively refine feature representations and filter out noisy data. gflu utilizes stage-specific feature masks, allowing for hierarchical and adaptive feature selection. This process results in a robust feature representation, which is then used by a mlp to produce final predictions [15].

For the hyperparameter tuning of all models, Optuna [2] is employed to facilitate an efficient hyperparameter search. We utilize tpe [30] as our optimization method, as it maintains two probability distributions: one for good configurations and another for all sampled configurations. By focusing on promising regions of the hyperparameter space, tpe directs the search more efficiently than traditional methods like grid or random search.

§ EXPERIMENTAL SETTING

§.§ Agents

In order to provide a comprehensive analysis of agent failures, it is essential to look at different agents. While there are numerous agents available for consideration, we decided to use the agents of [19] since their data was already available to us. The three agents in question are:

d_n The first agent does not interfere with the grid and only provides DoNothing actions. This agent is considered the baseline and we expect a lot more failures from it.

senior This agent is based on the CurriculumAgent [20], which selects topology actions on individual substations, when the threshold of $\rho_{max,t} >=0.95$ is breached. The agent further includes rule-based components, such as periodic reversion to the starting topology [19, 20] and line reconnection. The agent is more advanced, but changes the substations quite often in uncertain events.

topo The last agent is an extension of the senior and is described in detail in [19]. The main difference is that the agent searches for overall suitable tts when the $\rho_{max,t}$ is in the interval of $0.85<\rho_{topo}<0.95$. In [19], the topo was able to achieve better results, thus we expect differences to the senior.

§.§ General Data Collection

To ensure a in-depth  analysis, it is useful to look at a wide range of different environment scenarios. For this reason, we chose the WCCI 2022 l2rpn environment [27], which models the IEEE118 grid. The grid has a total of 1662 chronics, covering chronics from different periods of the year. The chronics of the WCCI environment were created to simulate an expected electricity mix of 2050 with only 3% fossil fuels remaining in the electricity mix [27]. With more re energy production, we expect much more variability in the grid, resulting in more diverse failures.
The Grid2Op environment contains several stochastic components that are randomized, such as an adversarial agent or the generation injections. We therefore decided to run each chronic a total of 10 times with different master seeds.[To ensure that there was no cherry picking, we chose the ten seeds completely random using ($8888$). For each seed, we also calculated the  and  separately.] We iterated over all chronics with each agent and discarded all runs in which the agent was able to complete the episode. The resulting data is then used for clustering and forecasting.

Cluster Data For the clustering approach, we want to compare the observation of the failure with the initial situation of the environment. Thus, we collected the observation at the start of each chronic $obs_{t=0}$ and one time step before failure $obs_{t=n-1}$ as well as the error statement for the failure. We collected a total of $39\,635$ observations, where the d_n agent had a total of $16\,474$ failures, the senior $11\,695$, and the topo $11\,466$. As each observation had $4295$ variables, it was necessary to aggregate them the follwing way before clustering.
First, since line capacity $\rho$ is an important indicator of grid stability, we decided to collect the maximum and average line capacity $\rho_{max}$ and $\rho_{mean}$ of the $obs_{t=n-1}$. Additionally, we counted the number of disconnected lines $\#lines_{dis}$ and the total $ts_{overflow}$, representing the cumulative duration during which the lines $l \in L$ have been in an overflown state ($\rho_l \geq 1.0$). Second, we want to account for the number of substation changes compared to $obs_{t=0}$. Thus we consider the number of substations that have been switched with $\#sub_{changed}$.
Third, it is necessary to account for the different power flow changes in the grid. For the load consumption, the generator injections, and both sides of the lines (origin and extremity), we record the maximal active load value $(load^p_{max},gen^p_{max},line^{ex,p}_{max},line^{or,p}_{max})$ and the average value $(load^p_{mean},gen^p_{mean},line^{ex,p}_{mean},line^{or,p}_{mean})$ as a fraction to the first observation $obs_{t=0}$.[We excluded reactive power flow and voltage angle due to high correlation with active power flow or constant values. Elements with a zero in $obs_{t=0}$ were also excluded to avoid mathematical errors. See appendix table <ref> for a summary of all variables.] The fraction is calculated for the loads, generators, and lines. As an example, we show the calculation for the $load^p_{max}$ and $load^p_{mean}$ with $\gamma$ being one load in the set of all loads $\gamma \in \Gamma$:

\begin{align}
load^p_{max} &= \max_{\gamma \in \Gamma}\left[load^p_{\gamma,t=n-1}/load^p_{\gamma,t=0}\right] -1 \\
load^p_{mean} &= \frac{1}{|\Gamma|} \sum^\Gamma_{\gamma=0} (load^p_{\gamma,t=n-1}/load^p_{\gamma,t=0} ) -1
\end{align}

Next to the cluster variables, we also collected the survival times ($t_{survived}$) as well as other descriptive variables for our analysis.

Forecasting Data

Our forecasting methods distinguish between failure and survival scenarios in power grid operation. For failures, we focus on observations at time steps $t=n-1$, $t=n-3$, and $t=n-5$. To identify surviving scenarios, we sample observations with very high $\rho_{max,t}$ that survive for more than $n-6$ steps. This approach yielded $189\,395$ observations, with half of those surviving ($obs_{t=survived}$), and the remaining half is split evenly among $obs_{t=n-1}$, $obs_{t=n-3}$, and $obs_{t=n-5}$.

Subsequently, the data was divided into a training set of size $153\,409$, a validation set of size $17\,046$, and a testing set of size $18\,940$. Thanks to the courtesy of RTE, the validation environment of the WCCI 2022 challenge was made available to us. Consequently, our findings can also be validated on an ood sample of size $6036$.

§.§ Metrics

We evaluate the performance of our prediction models using several key metrics. First, we include the Accuracy, which is the ratio of correctly predicted instances to the total instances. However, as we are in a multi-class classification problem, we also include the Balanced Accuracy to address the imbalance in the data sets. It is defined as the average of recall obtained for each class. The F1 score, which can be interpreted as a harmonic mean of precision and recall, reaches its best value at 1 and worst score at 0, with the relative contributions of precision and recall to the F1 score being equal. For our specific use-case, we also want to measure the Binary Accuracy. This metric is defined by the accuracy of the one-vs-all model that categorizes all failures into a single class and survival into another. It effectively assessing the model's ability to distinguish between survival and failure observations.[Note that no additional training is conducted for any specific metric, including the Binary Accuracy. Instead, the predictions of the trained multi-class models are aggregated.]

Mean values of the clustering variables, grouped by one of the five clusters. Based on the statistical characteristics, we name the clusters accordingly. We highlight the highest mean value of each variable. The clusters are described below.

Cluster      Changed
Topology        Decreased
Consumption     Disconnected
Power Lines     Increased
Injections     Increased
Power Flow
on Power Lines

$ts_{overflow}$     5.32     6.89     12.60     10.32     9.96

$\#lines_{dis}$     0.79     1.58     3.93     1.37     2.12

$\#sub_{changed}$     5.54     1.71     0.92     2.94     3.57

$\rho_{max}$     1.03     1.13     1.43     1.25     1.24

$\rho_{mean}$     0.40     0.38     0.36     0.41     0.38

$load^p_{max}$     0.33     0.09     0.35     0.34     0.34

$load^p_{mean}$     0.12     -0.06     0.13     0.14     0.13

$gen^p_{max}$     4.78     2.69     3.66     614.55     2.31

$gen^p_{mean}$     0.27     0.10     0.08     125.29     -0.07

$line^{ex,p}_{max}$     0.51     0.21     0.60     0.58     1.33

$line^{ex,p}_{mean}$     19.16     12.24     22.07     19.56     106.85

$line^{or,p}_{max}$     0.51     0.20     0.61     0.58     1.33

$line^{or,p}_{mean}$     18.96     12.27     22.44     19.10     108.58

§ RESULTS

§.§ Clustering Results

Cluster Identification:
With the aggregated clustering data, we perform pca and select six principal components, describing 85% of the variance. Thereafter, we execute the k-means clustering and identify five clusters as the best cluster size for our data.[The plot of both inertia and silhouette can be found in the appendix figure <ref>.]
Subsequently, the clusters are grouped and the statistical components, including the mean, median, and quantiles of the variables, are analyzed to identify the distinctive characteristics of the clusters. We provide in Table <ref> the mean values of all cluster variables.
The following clusters were identified:

* Changed Topology: This cluster has the highest $\#sub_{changed}$ among all clusters.

* Decreased Load Consumption: In this cluster the distinct characteristic is the lower load consumption of $load^p_{mean}$ and $load^p_{max}$.

* Disconnected Power Lines: With the highest $ts_{overflow}$, $\#lines_{dis}$ and $\rho_{max}$ this cluster is directly linked to the overload of the power lines.

* Increased Generator Injections: This cluster has by far the highest generator values of $gen^p_{mean}$ and $gen^p_{max}$.

* Increased Power Flow on Power Lines: As the last cluster we identify the highest power flow in $line^{ex,p}_{max},line^{or,p}_{max},line^{ex,p}_{mean}$ and $line^{or,p}_{mean}$ as distinct feature.

3D visualization of the five clusters, where each point represents a failure of the agents. The axis in the plot are the $gen^p_{mean}$ on the y-axis, $\#lines_{dis}$ on the x-axis and $\#sub_{changed}$ as z-axis. The points are colored according to their respective clusters.

Interestingly, we can see in Figure <ref> that the clusters Increased Generator Injections and Disconnected Power Lines have distinct observations that can clearly be separated from the other clusters. In contrast, there seems to be multiple clusters with an increased $\#sub_{changed}$. With respect to Table <ref>, several noteworthy observations can be made. First, there appears to be a significant increase in generator output in the fourth cluster, with a factor of up to 125 times the $gen^p_{mean}$. Given that the relative change to $t=0$ in active power is compared, it can be explained that some generators only produce a minor amount at $t=0$ due to the fact that all scenarios start at nighttime.
Secondly, it can be observed that the only instance of a decrease in the mean values can be found in the second cluster with $load^p_{mean}$.
Note that Kirchhoff's current law has to hold, i.e., the aggregated loads always have correspond to the aggregate generation in the whole grid. This means that the lower $load^p_{mean}$ can only show localized imbalances in parts of the grid.

Cluster distribution across Agents:
Following the general description of the clusters, we are interested in the overall size of the clusters and their distribution across the agents. In Table <ref> we summarize these results. First, we can see that there is a clear size difference between the clusters. Changed Topology and Disconnected Power Lines are quite large containing 39% and 30% of the data, while the Increased Generator Injections only has a total of 405 samples. With the $\chi^2$ independence test ($\alpha=0.05$), we tested the $H_0$-hypothesis that the clusters and agents are independent from each other. We could clearly reject the hypothesis with a $p-value$ of $0.0$, which is a very strong indication that there is an association between clusters and agents.
Comparing the clusters with the agents, we can see that a large portion of the d_n failures are due to the disconnection of power lines (56%) as well as the reduced load consumption (26%). In contrast, the senior and topo have most of their failures in the topology cluster. Here, the topo has with  58% fewer topology failures, compared to the senior's 61%. Instead, the topo has 3% more load consumption errors than the senior.

Quantity and Frequency of failures in the clusters for the three agents. We denote the number of occurrences per cluster and the percentage within brackets. Further, we highlight for each agent the most frequent cluster. Note that the amount of failure data across the agents differs.

d_n    senior    topo     Total

Topology        1573 (10%)      7191 (61%)      6676 (58%)     15440 (39%)

Decreased Load
Consumption     4345 (26%)      2239 (19%)    2508 (22%)    9092 (23%)

Power Lines     9225 (56%)     1360 (12%)    1424 (12%)    12009 (30%)

Increased Generator
Injections     178 (1%)     121 (1%)   106 (1%)    405(1%)

Increased Power
Flow on Power Lines     1153 (7%)       784 (7%)   752  (7%)    2689 (7%)

Total # Obs     16.474    11.695     11.466     39.635

Cluster Survival Time.
Besides the distribution across the agents, we are also interested in the survival time of the different clusters. In Figure <ref> we visualize the number of steps until failure with $t_{survived}$. We observe a clear difference between the clusters, particularly, we can see that the Increased Generator Injections cluster fails incredibly fast, with a median survival time of only half a day (one day equals 288 steps). The median of the Decreased Load Consumption cluster is a little better, although it is still less than a full day. On the other hand, we see that the Changed Topology cluster has the longest survival time with a median value of 496 steps, which might be related to the capabilities of the topo and senior. Another interesting result is that we have a relatively similar distribution between the remaining two clusters, which is understandable since both focus on line results.

Box plot of the survival time $t_{survived}$ of each cluster. We further report the median survival time. The ticks on the x-axis correspond to a full day, e.g., 288 steps. Note that an agent must survive a total of 2016 time steps to complete the chronics successfully.

§.§ Forecasting Results

Quantitative forecasting results
In Table <ref>, we present the performance of all evaluated models with regard to the defined metrics of Section <ref>. lgbm outperformed the other models across all metrics, achieving an accuracy of $86\%$ and a high binary accuracy for predicting survival vs. failure of $91\%$. Furthermore, the performance on the ood dataset decreases as expected, however only by $6\%$ for the balanced accuracy and $5\%$ for the binary accuracy. The latter is especially promising since a $86\%$ accuracy on unseen ood data suggests robust generalization capabilities when the model encounters novel scenarios.
In second place, the other gradient boosting approach xg achieved a relatively good performance, however, it was not able to reach that of lgbm.

Compared to the gradient boosting approaches, the neural network-based methods cem and gandalf achieve an accuracy approximately $6\%$ and $7\%$ lower than lgbm, respectively. Both methods perform very similarly with cem being at third place with an accuracy of $80\%$. Interestingly, both achieve a high accuracy for the binary problem of predicting survival vs. failure, with gandalf even slightly outperforming xg. It is also notable that both neural network-based methods achieve overall similar scores on both the test dataset and the ood dataset. However, the scores are still lower than those achieved by lgbm but comparable to xg. Therefore, the gandalf architecture did not achieve a higher performance than the cem model, though it is more complex. In the last place is rf, which performed the worst across all metrics.

Forecasting results of the different models on the hold-out test data. The last two columns show the balanced accuracy and binary accuracy results on the ood test data.

accuracy     balanced
accuracy     f1 micro     binary
accuracy     ood balanced
accuracy     ood binary

rf     0.75     0.64     0.75     0.78     0.62     0.82

xg     0.82     0.77     0.82    	0.87     0.72     0.82

lgbm     0.86     0.82     0.86     0.91     0.76     0.86

cem     0.80     0.74     0.80     0.87     0.72     0.84

gandalf     0.79     0.73     0.79     0.88     0.71     0.84

With regards to the imbalanced nature of the data, it is important to also look at the balanced accuracy and the f1 scores. As expected, we observe a decrease in balanced accuracy compared to the regular accuracy score, however, the best performing model lgbm has the lowest deviation of approximately $4\%$, indicating a higher robustness to the class imbalance.

To further understand the implications of these quantitative results, we delve into a qualitative analysis. This examination allows us to uncover the underlying reasons behind the model's performance and reliability. By analyzing the probability distributions and feature importance, we gain deeper insights into how the lgbm model makes its predictions, and identify crucial input features that impact the decision-making of the model.

Qualitative Forecasting result

Average probability distribution of lgbm for the ground truth of $obs_{t=survived},obs_{t=n-5},obs_{t=n-3}$ and $obs_{t=n-1}$. The probability output is averaged for all observations. Black lines visualizes the kernel density estimation across the classes.

We take a closer look at the probability distribution of the best performing model lgbm. Figure <ref> shows the averaged probability output of the model against the ground truth of all four classes. As we can see, the model selects the correct classification on average and is especially sure in the cases of survival, i.e., $obs_{survived}$ with $95\%$, and imminent failure, i.e., $obs_{t-1}$ with $87\%$. In the case of $obs_{t=n-5}$ we can see some uncertainty in the model, since the model selects the surviving $obs_{t=survived}$ in $\frac{1}{4}$ of the cases.
This shows that in some cases the grid appears to be stable five moves before the grid collapses. Thus, it makes it difficult to distinguish from a surviving observation, as it takes a sudden change, such as an adversarial attack, to destabilize the grid very quickly. Nevertheless, we can detect failures with a higher certainty at three time steps before the failure, thus offering some time to counteract. Putting these results into the operator's perspective, with a time step equal to 5min real time, we would send a first warning at $obs_{t=n-5}$ (25min) and raise the alarm at $obs_{t=n-3}$ (15min before failure).

Feature Importance
As a last result, we examine the feature importance of the lgbm model. For this we use the gain metric, which represents the improvement in accuracy brought by a feature to the branch. By evaluating the gain, we identify which features contribute the most to the model's predictions, allowing for a better understanding of the underlying patterns and relationships within the data.

Feature importance of the 30 most important features according to the lgbm model. The color indicates the type of variable.

With respect to the features of the lgbm, we look at the top 30 features, visualized in Figure <ref>. We group the features based on their type and distinguish among descriptive, generator, load, and line features (power flow, cooldown, $ts_{overflow}$, maintenance).
Beginning with the interpretation of these results, we observe that $16$ of the most important features are of type $ts_{overflow}$ which denotes the number of time steps since the power line has $\rho_{l}\geq 1.0$. This is not surprising, since overloaded lines are a strong indicator for an unstable grid. However, we also have three important descriptive features ($current\_step$, $minute\_of\_hour$, and $hour\_of\_day$). These suggest that the likelihood of power grid failures is influenced by temporal features, reflecting patterns that may correlate with daily operational cycles or varying demand levels throughout the day. Another important descriptive feature is the agent type in 17th place, since there is a clear difference in behavior between the topo, senior and the d_n in terms of survival. It is plausible that some grid states would lead to a grid failure when no remedial action is taken, whereas with an agent present to apply topological remedial actions, these states might not result in a grid failure.
Regarding the loads and generators, we can see that the first generator feature appears at position 10, while the first load feature appears only at position 26. This is quite interesting, considering that the decreased load cluster is relatively frequent. It is possible that the grid's stability is more dependent on the fluctuating nature of largely weather-dependent renewable generators than on loads.

Examining the overall picture, it is unsurprising that the feature importance is dominated by line features. Moreover, many features correspond to the same lines. For instance, three distinct features corresponding to line 175 appear among the top 30 importance values alone, while line 17 has two features in the top 30. To simplify the interpretation of feature importance and to better understand the broader impacts of different types of grid components, we aggregated the feature importance values by averaging the importance values of related features. The results of the aggregation are visualized in Figure <ref>, where we display the ten most important lines, generators and loads of the grid.

Top 10 important lines (red), generators (yellow), and loads (green) for failure prediction. The numbers on the marked lines correspond to the line ids of Figure <ref>. Sub-grids are separated with dotted lines.
For better clarity, we grayed out the less important elements of the grid and highlighted 3 significant regions (A,B,C) of important grid features.

In the first region (A) we have three loads and a line of high importance. By looking at the failures we could identify the line between substations 21 and 22 as reason, as it is often targeted by the adversarial agent. This attack is a certain failure because the demand of the loads cannot be changed and the only generator that provides electricity is in substation 18. This causes an overflow in all the lines between substations 18 and 21, leading to a cascading failure in less than 3 steps. Thus, by disconnecting the line and not allowing it to be reconnected in time, the adversarial agent causes a game over, which to our surprise is very common with 4057 out of 39635 cases.

Region (B) highlights the connection between two sub-grids: the line between substations 68 and 76 connects the upper right sub-grid to a major substation in the bottom right sub-grid, making it crucial as the most important generator 37 connects to substation 68 and frequently injects up to 500% of its starting power. Furthermore, substation 76 is connected to the central substation 79 with two lines, indicating an increased power flow and overall significance of this connection. Accordingly, these lines are deemed as very important by the prediction model. Both substations are very important since all power flow to the majority of this sub-grid have to be routed through them, hence they serve as the entry point to this part of the grid. Additionally, two other paths (80-67-64 or 80-67-68) from substation 79 to the upper right sub-grid have to traverse substation 80, hence the connection between 79 and 80 is also considered very important by the model. Overall, region (B) can be characterized as a high line power flow region, having 2 important generators (37 and 42) in the vicinity, connecting two sub-grids, and being the only entry point to the bottom right area of the grid.

Region (C) marks a section of the grid where frequent line attacks by the adversarial agent occur. In fact, 10 out of 23 possible line attacks happen in this sub-grid with 5 line attacks in region (C) alone. Both the lines connecting substations 58 to 62 and 62 to 63 are part of the attacked lines, leading to a higher load on the other lines connecting the group of substations ($59,60,61,63,64,65,66$) in the lower area of this sub-grid. Namely, these are the lines connecting substation 58 to 59 and 60, respectively, as correctly identified by the importance of the model. Moreover, three out of the ten most important loads are present in this region, connected at substations 50, 52, and 56, respectively. In order to supply these loads, a high load and importance at the lines connecting these loads to bigger substations is to be expected. This is exacerbated by the fact that the lines connecting substation 48 to 50 and substation 48 to 53 can be attacked, leaving less possible paths to these loads available.

§ DISCUSSION

Looking at the results, we can see that the clustering was quite successful.
We were able to identify at least two of the five clusters that are highly correlated with the agents, confirming our expectations. Moreover, it is interesting to note that both advanced agents were able to significantly reduce the number of errors due to line disconnection. We can also confirm the idea of [19], since the topo has slightly less errors due to topology changes.

In terms of the remaining clusters, the results are not as clear. For the line and generator clusters, we have very similar cluster sizes across all agents. For the load cluster, we could at least observe a slight reduction for the advanced agents. Overall, the results indicate that there is still a lot of optimization potential, as both the load and generator clusters do not survive a full day in the median. This shows a clear need for some redispatching optimization. With the feature importance, we were also able to detect the design error on line 93, showing a possible requirement to incorporate load shedding in future action spaces of Grid2Op. Agents should focus on preventing the disconnection of this line due to overflows at all cost by performing suitable remedial actions in a timely manner.

With respect to our predictive models, we were able to achieve a balanced accuracy of 86%, showing that model correctly detect possible failures in the gird in most cases.  More importantly, the binary accuracy of 91% showed, that we are able to correctly detect whether the agent is struggling in the grid. This makes the real-world applicability more interesting, as we can alert the operator in time, with $obs_{t=n-5}$ corresponding to 25min before failure. Furthermore, as our Figure <ref> shows, we get more and more certain about the state of the grid, thus reducing the risk of a false alarm.

As future work, improving the performance of the failure prediction could be achieved by utilizing gnn to capture the inherent graph structure present in the power grid topology. Further, a continuous analysis should be conducted that includes redispatching components to target the reduced load cluster by actively adjusting redispatch and curtailment. Furthermore, it might be interesting to focus on the specific regions of the grid that received high feature importance, as they seem to be responsible for the performance of the agents. This could be achieved by specifically training the drl agent on these regions, or even a hierarchical approach. Including topological actions that specifically target these regions into the action space might bring a boost in performance as well. Finally, the prediction of failures has a potential in guiding decisions on when a DRL agent should take remedial actions based on them, rather than fixed thresholds of maximum line capacity $\rho_{max,t}$.

§ CONCLUSION

In this paper, we provide a detailed failure analysis on the WCCI l2rpn environment, which has not been done to this extent yet. We examined the failures in the grid of three different agents across ten seeds, resulting in 40k data points. Our analysis is twofold, with a detailed cluster analysis and a prediction of the failures in advance. For the cluster analysis, we were able to identify five specific clusters that showed distinct causes for the failures. The failures exhibited varying survival times and could be attributed to some extend to specific agents. As a second part, we propose a multi-class forecasting approach to detect the failures ahead of time. On the data we tested five different model types and found the lgbm to be the most suitable model for the prediction. The feature importance afterwards revealed critical regions of the grid that could be specifically targeted in future work.

§ ACKNOWLEDGEMENT

This work was supported by:

* Reinforcement Learning for Cognitive Energy Systems (RL4CES) from the Intelligent Embedded Systems of the University Kassel and Fraunhofer IEE.

* Graph Neural Networks for Grid Control (GNN4GC) founded by the Federal Ministry for Economic Affairs and Climate Action Germany under the funding code 020E-100626677.

* AI4REALNET has received funding from European Union’s Horizon Europe Research and Innovation program under the Grant Agreement No 101119527. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them.

§ DISCLOSURE OF INTEREST

The authors have no competing interests to declare that are
relevant to the content of this article.

[1]
Abdi, H., Williams, L.J.: Principal component analysis. Wiley interdisciplinary reviews: computational statistics  2(4),  433–459 (2010)

[2]
Akiba, T., Sano, S., Yanase, T., Ohta, T., Koyama, M.: Optuna: A next-generation hyperparameter optimization framework. In: Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2019)

[3]
Amdouni, E., Khouadjia, M., Meddeb, M., Marot, A., Crochepierre, L., Achour, W.: Grid2onto: An application ontology for knowledge capitalisation to assist power grid operators. In: International Conference On Formal Ontology in Information Systems-Ontology showcases and Demos (2023)

[4]
Bacher, R., Glavitsch, H.: Network topology optimization with security constraints. IEEE Transactions on Power Systems  1(4),  103–111 (1986)

[5]
Chauhan, A., Baranwal, M., Basumatary, A.: Powrl: A reinforcement learning framework for robust management of power networks. arXiv preprint arXiv:2212.02397  (2022)

[6]
Chen, K., Huang, C., He, J.: Fault detection, classification and location for transmission lines and distribution systems: a review on the methods. High voltage  1(1),  25–33 (2016)

[7]
Chen, T., Guestrin, C.: XGBoost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. pp. 785–794. KDD '16, ACM, New York, NY, USA (2016). 10.1145/2939672.2939785, <http://doi.acm.org/10.1145/2939672.2939785>

[8]
Dashti, R., Daisy, M., Mirshekali, H., Shaker, H.R., Aliabadi, M.H.: A survey of fault prediction and location methods in electrical energy distribution networks. Measurement  184,  109947 (2021)

[9]
Ding, C., He, X.: K-means clustering via principal component analysis. In: Proceedings of the twenty-first international conference on Machine learning. p. 29 (2004)

[10]
Donnot, B.: Grid2op- A testbed platform to model sequential decision making in power systems. . <https://GitHub.com/rte-france/grid2op> (2020), accessed 22-Jan-2023 on Github

[11]
Dorfer, M., Fuxjäger, A.R., Kozak, K., Blies, P.M., Wasserer, M.: Power grid congestion management via topology optimization with alphazero. arXiv preprint arXiv:2211.05612  (2022)

[12]
EI Innovation Lab, Huawei Cloud, Huawei Technologies: NeurIPS Competition 2020: Learning to Run a Power Network (L2RPN) - Robustness Track. <https://github.com/AsprinChina/L2RPN_NIPS_2020_a_PPO_Solution> (2020), accessed 22-Jan-2023 on Github

[13]
Ho, T.K.: Random decision forests. In: Proceedings of 3rd international conference on document analysis and recognition. vol. 1, pp. 278–282. IEEE (1995)

[14]
Joseph, M.: Pytorch tabular: A framework for deep learning with tabular data (2021)

[15]
Joseph, M., Raj, H.: Gandalf: Gated adaptive network for deep automated learning of features (2024)

[16]
Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., Liu, T.Y.: Lightgbm: A highly efficient gradient boosting decision tree. Advances in neural information processing systems  30,  3146–3154 (2017)

[17]
Kelly, A., O'Sullivan, A., de Mars, P., Marotzhou2021actiont, A.: Reinforcement learning for electricity network operation. arXiv preprint arXiv:2003.07339  (2020)

[18]
Kodinariya, T.M., Makwana, P.R., et al.: Review on determining number of cluster in k-means clustering. International Journal  1(6),  90–95 (2013)

[19]
Lehna, M., Holzhüter, C., Tomforde, S., Scholz, C.: Hugo–highlighting unseen grid options: Combining deep reinforcement learning with a heuristic target topology approach. arXiv preprint arXiv:2405.00629  (2024)

[20]
Lehna, M., Viebahn, J., Marot, A., Tomforde, S., Scholz, C.: Managing power grids through topology actions: A comparative study between advanced rule-based and reinforcement learning agents. Energy and AI  14,  100276 (2023)

[21]
MacQueen, J., et al.: Some methods for classification and analysis of multivariate observations. In: Proceedings of the fifth Berkeley symposium on mathematical statistics and probability. vol. 1, pp. 281–297. Oakland, CA, USA (1967)

[22]
Marot, A., Donnot, B., Chaouache, K., Kelly, A., Huang, Q., Hossain, R.R., Cremer, J.L.: Learning to run a power network with trust. Electric Power Systems Research  212,  108487 (2022)

[23]
Marot, A., Donnot, B., Dulac-Arnold, G., Kelly, A., O’Sullivan, A., Viebahn, J., Awad, M., Guyon, I., Panciatici, P., Romero, C.: Learning to run a power network challenge: a retrospective analysis. In: NeurIPS 2020 Competition and Demonstration Track. pp. 112–132. PMLR (2021)

[24]
Marot, A., Donnot, B., Romero, C., Donon, B., Lerousseau, M., Veyrin-Forrer, L., Guyon, I.: Learning to run a power network challenge for training topology controllers. Electric Power Systems Research  189,  106635 (2020a)

[25]
Marot, A., Guyon, I., Donnot, B., Dulac-Arnold, G., Panciatici, P., Awad, M., O’Sullivan, A., Kelly, A., Hampel-Arias, Z.: L2rpn: Learning to run a power network in a sustainable world neurips2020 challenge design (2020b)

[26]
Mu, Z., Xu, P., Zhang, K., Gao, T., Zhang, J.: Cascading fault early warning and location method of transmission networks based on wide area time-series power system state. IEEE Journal of Radio Frequency Identification  7,  6–11 (2022)

[27]
Serré, G., Boguslawski, E., Donnot, B., Pavão, A., Guyon, I., Marot, A.: Reinforcement learning for energies of the future and carbon neutrality: a challenge design. arXiv preprint arXiv:2207.10330  (2022)

[28]
Subramanian, M., Viebahn, J., Tindemans, S.H., Donnot, B., Marot, A.: Exploring grid topology reconfiguration using a simple deep reinforcement learning approach. In: 2021 IEEE Madrid PowerTech. pp. 1–6 (2021). 10.1109/PowerTech46648.2021.9494879

[29]
Viebahn, J., Naglic, M., Marot, A., Donnot, B., Tindemans, S.H.: Potential and challenges of AI-powered decision support for short-term system operations. CIGRE 2022 Paris Session  (2022)

[30]
Watanabe, S.: Tree-structured parzen estimator: Understanding its algorithm components and their roles for better empirical performance (2023)

[31]
Zhou, B., Zeng, H., Liu, Y., Li, K., Wang, F., Tian, H.: Action set based policy optimization for safe power grid management. In: Machine Learning and Knowledge Discovery in Databases. Applied Data Science Track: European Conference, ECML PKDD 2021, Bilbao, Spain, September 13–17, 2021, Proceedings, Part V 21. pp. 168–181. Springer (2021)

§ APPENDIX

§.§ Cluster Variables

We summarize the distribution of the cluster variables in the following Table <ref>.

Name     Explanation     mean     std     min     25%     50%     75%     max

$ts_{overflow}$     Cumulative
overflow time    8.25     6.78     0.0     3.0     6.0     12.0     71.0

$\#lines_{dis}$      Number of
disconnected lines     2.02     2.27     0.0     0.0     1.0     3.0     20.0

$\#sub_{changed}$     Number of changed
Substations      3.10     4.63     0.0     0.0     1.0     4.0     34.0

$\rho_{max}$     Max line capacity     1.19     0.28     0.61     1.03     1.11     1.33     2.0

$\rho_{mean}$     Mean line capacity     0.38     0.06     0.19     0.33     0.38     0.42     0.63

$load^p_{max}$     Max active load
value of consumption      0.28     0.14     -0.14     0.22     0.30     0.37     0.88

$load^p_{mean}$     Mean active load
value of consumption      0.08     0.10     -0.26     0.04     0.11     0.15     0.32

$gen^p_{max}$      Max active production
value of genreators      10.03     70.03     -0.48     0.79     1.67     3.27     1573.0

$gen^p_{mean}$     Mean active production
value of genreators     1.43     13.40     -0.98     -0.24     -0.0     0.30     224.02

$line^{ex,p}_{max}$    Max active power
flow at line extremity     0.52     0.39     -0.22     0.24     0.46     0.73     4.60

$line^{ex,p}_{mean}$     Mean active power
flow at line extremity      24.41     32.13     0.24     8.60     15.89     27.25     696.11

$line^{or,p}_{max}$     Max active power
flow at line origin       0.53     0.39     -0.22     0.24     0.46     0.74     5.64

$line^{or,p}_{mean}$     Mean active power
flow at line origin       24.56     33.86     0.23     8.57     15.84     27.15     730.49

$t_{survived}$     Survival Times    539.87     470.42     2.0     147.0     407.0     804.0     2016.0

Variables for clustering. Note that $t_{survived}$ is not included in the clustering and is just reported to showcase the distribution.

Correlation of cluster variables: When looking at the correlation plot in Figure <ref>, there are already three interesting things to note. First, we see a high correlation in the power flow values between the maximum and average components, which is relatively self-explanatory. Second, we can see that the load values are somewhat correlated to $line^{ex,p}_{max},line^{or,p}_{max}$. Interestingly, the generator values are not correlated with the line power flow in any way.  Third, we see a negative correlation between $\#sub_{changed}$ and $\#lines_{dis}$, which can be explained by the fact that both the senior and topo agents have components that automatically reconnect lines. However, $\#lines_{dis}$ is also negatively correlated with $\rho_{max}$, which is counter-intuitive. A possible explanation could be the Grid2Op rule that disconnects a line after three consecutive time steps of overload. When the line is disconnected, the highest $\rho_{max}$ is no longer available, explaining the negative correlation. This explanation would primarily hold for consecutive observations, so it is quite interesting that it seems to work across multiple scenarios.

Correlation plot of the different cluster variables. The color range depicts the correlation values that go from -1.0 (blue) to 1.0 (yellow).

§.§ Inertia and Silhouette plot

The plot below (Fig. <ref>) illustrates the inertia values for different cluster sizes in our k-means clustering analysis, showcasing the selection process for determining the optimal number of clusters based on the elbow method and the corresponding decrease in inertia.

Inertia and Silhouette plot of different cluster sizes. As one can see, the inertia has a "elbow" at a cluster size of five. Similar, the silhouette shows one of its highest values at a size of five as well. Thus, we select a cluster size of five.

§.§ Additional Information

Error Statement
We collected the error statement of every failed scenario to see, whether we could detect a pattern. The following three error statements were outputted by the underlying Grid2Op solver :

We collected all error messages and mapped them to the clusters of Section <ref>. The results can be found in Table <ref>. While there is some difference in the percentage change, we can not find a distinct relation between cluster and error message.

Type      Changed
Topology        Decreased Load
Consumption     Disconnected
Power Lines     Increased Generator
Injections     Increased power flow
on Power Lines

1     61.2%     63.97%     69.27%     68.15%     65.01%

2     27.47%     27.03%     24.71%     26.42%     26.18%

3     11.32%     9%     6.01     5.43%     8.81%

In this Table, we map the three error types with the different clusters. For each cluster we aggregate the error types to percentage values for better viability.

Specific Descriptive Scenarios
We further accounted for three descriptive variables for the scenarios:

* $Attack_{93}$: This variables indicates, whether the Line 93 is disconnected or will be disconnected in the next step

* This variable indicates, whether the grid was in its original state, i.e., all substation had their buses set to one and all lines were disconnected

* Lastly, we collected the percentage change of the generator 37, as this generator was directly at substation 68.

For the first two variables, we have the following occurence:

§.§ Specific Descriptive Scenarios

d_n    senior    topo     Total

$Attack_{93}$     1367    1312      1378    4057

$Stable$     169     817   646    1632

Specific failure cases identified in descriptive analysis.

In case of the generator 37, we can see the following percentage increase across the scenarios:

Visualization of the percentage change on the generator 37 in comparison to its $obs_{t=0}$.


# Probing the anharmonicity of the potential well for magnetic vortex core in
a nanodot

O.V. Sukhostavets  Departamento de Física de Materiales, Universidad del Pais
Vasco, 20018 San Sebastian, Spain  B. Pigeau  Service de Physique de l'État
Condensé (CNRS URA 2464), CEA Saclay, 91191 Gif-sur-Yvette, France  S. Sangiao
Service de Physique de l'État Condensé (CNRS URA 2464), CEA Saclay, 91191 Gif-
sur-Yvette, France  G. de Loubens  Service de Physique de l'État Condensé
(CNRS URA 2464), CEA Saclay, 91191 Gif-sur-Yvette, France  V.V. Naletov
Service de Physique de l'État Condensé (CNRS URA 2464), CEA Saclay, 91191 Gif-
sur-Yvette, France  Institute of Physics, Kazan Federal University, Kazan
420008, Russian Federation  O. Klein<EMAIL_ADDRESS>Service de Physique de
l'État Condensé (CNRS URA 2464), CEA Saclay, 91191 Gif-sur-Yvette, France  K.
Mitsuzuka  S. Andrieu  F. Montaigne  Institut Jean Lamour, UMR CNRS 7198,
Université de Lorraine, 54 506 Nancy, France  K.Y. Guslienko  Departamento de
Física de Materiales, Universidad del Pais Vasco, 20018 San Sebastian, Spain
IKERBASQUE, The Basque Foundation for Science, 48011 Bilbao, Spain

###### Abstract

The anharmonicity of the potential well confining a magnetic vortex core in a
nanodot is measured dynamically with a Magnetic Resonance Force Microscope
(MRFM). The stray field of the MRFM tip is used to displace the equilibrium
core position away from the nanodot center. The anharmonicity is then inferred
from the relative frequency shift induced on the eigen-frequency of the vortex
core translational mode. An analytical framework is proposed to extract the
anharmonic coefficient from this variational approach. Traces of these shifts
are recorded while scanning the tip above an isolated nanodot, patterned out
of a single crystal FeV film. We observe $+10$% increase of the eigen-
frequency when the equilibrium position of the vortex core is displaced to
about one third of its radius. This calibrates the tunability of the
gyrotropic mode by external magnetic fields.

There has been recently a renewed interest, both theoretical and experimental,
in the problem of nonlinear (NL) magnetization dynamics inside confined
nanostructures Slavin2009 . NL phenomena are responsible for the creation of
novel dynamical objects Mohseni2013 , analogs of dynamical solitons. They also
set the figure of merit of spintronics devices, e.g. the spectral purity and
tuning sensitivity of spin transfer torque nano-oscillators Slavin2009 . On
the theoretical side, predictions on the amplitude of the NL coefficients have
been found to be extremely difficult to compute beyond the uniformly
magnetized ground state. The difficulty raises both from the magneto-dipolar
field, which introduces a non-local interaction, and from the kinetic part of
the effective field (or gauge field), which modifies the texture of the
magnetic configuration. On the experimental side, the most promising findings
have been discovered on non-uniform ground states, such as magnetic vortex
existing in ferromagnetic nanodot. Vortices have stimulated the emergence of
higher performance microwave oscillators using isolated Pribiag2007  or
dipolarly coupled Locatelli2011  nanodot, or for future magnetic memories by
allowing the resonant switching of the magnetic configuration Pigeau2010 .

Magnetic vortex corresponds to a curling in-plane magnetization spatial
distribution leaving a nanometric in size core region ($\sim$ the exchange
length), where the magnetization is pointing out-of-plane. The lowest energy
mode is a translational (or gyrotropic) mode of the vortex core position
$\bm{X}$, expressed here in reduced unit of the dot radius. Its properties are
governed by the magnetostatic potential well $W^{\text{(M)}}$ in which the
core evolves. For a circular nanodot, the magnetostatic energy is isotropic in
the dot plane and it can be written as a series expansion of even powers of
the dimensionless $\bm{X}$ Dussaux2012 :

$W^{\text{(M)}}=W^{\text{(M)}}_{0}+\frac{1}{2}\kappa|{\bm{X}}|^{2}+\frac{1}{4}\kappa^{\prime}|{\bm{X}}|^{4}+{\cal{O}}(|{\bm{X}}|^{6})\,,$
(1)

At the present, only the parabolicity of the confinement, $\kappa$, has been
well characterized experimentally and the measured value is in agreement with
theoretical predictions Novosad2005 . This is not the case for the higher
order terms and there is no consensus yet on the order of magnitude or the
sign of the anharmonic coefficient $\lambda\equiv\kappa^{\prime}/\kappa$
afferent to the depolarisation effect of a displaced vortex inside a large
planar circular nanodot. The asymptotic limit of large radius is the relevant
aspect ratio to test the dipole dominating limit and the circular symmetry is
necessary to avoid the additional complexity of non-isotropic potential found
for example in square shaped elements Drews2012 . Up to now, attempts to
measure $\lambda$ in circular nanodot using large rf excitation have so far
lead to inconsistent results between experiments Pigeau2011  (red shift) and
theory Gaididei2010  (blue shift). The measurement of $\lambda$ through a
variational approach, consisting in studying the small change of oscillation
period when a large static displacement of the vortex core equilibrium
position is produced, has so far failed too: this counterpart of large rf
oscillation has mostly revealed the potential well inhomogeneities leading to
pinning of the core Chen2012 ; Burgess2013 .

In this work, we report on an experimental measurement of $\lambda$ in a large
planar circular nanodot using a Magnetic Resonance Force Microscope (MRFM).
All the experimental measurements are performed on an individual nanodisk of
thickness $t=26.7$ nm thick and nominal radius $R=300$ nm, patterned out of a
single crystal FeV film. Only the perfect crystalline structure ensures an
unpinned displacement of the vortex core throughout the sample volume. We rely
here on the non-uniform stray field of the magnetic tip of the MRFM to
displace the vortex core away from the nanodot center. The anharmonic
coefficient is then inferred from the measurement of the relative variation of
the eigen-frequency of the gyrotropic mode as a function of the tip
displacement.

Figure 1: (Color online) a) Side and b) top views of the experimental setup:
the stray field of an MRFM tip is used to displace the static position
$\bm{X}_{0}$ of the vortex core. The magnetic vortex state is shown in a bi-
variate colormap of $\bm{m}=\bm{M}/M_{s}$ (amplitude-phase $\leftrightarrow$
luminance-hue). The insets are microscopy images of the magnetic tip and disk
sample.

We start first with a description of the experimental setup Klein2008
(FIG.1). The right inset shows an image of the sample: a circular nanodot,
which is patterned by standard lithography and ion-milling techniques from an
extended film of Fe-V (10% V)Mitsuzuka2012  with magnetization $4\pi M_{s}=17$
kG. A magnetic tip is brought in the vicinity of the sample (left image). The
tip consists of a soft Fe particle glued at the apex of micro-cantilever. The
MRFM is placed inside a superconducting coil magnet, which produces an
homogenous bias magnetic field $\bm{H}_{0}$ of 6 kOe along the $z$-direction
(parallel to the normal of the disk). The value of ${H}_{0}$ is chosen to be
strong enough to magnetize the MRFM tip close to its saturation value, while
remaining weak enough compared to the saturation field of the nanodot to
preserve the vortex ground state inside the sample. At ${H}_{0}$ the tip stray
field is
$\bm{H}^{\text{tip}}(\bm{r})=-\nabla(\bm{\mu}_{\text{tip}}\cdot\bm{r}/r^{3})$,
the dipolar field generated by a point magnetic moment
${\mu}_{\text{tip}}=4\times 10^{-10}$ emu oriented along ${\hat{\bm{z}}}$
Pigeau2012 . Effects of perpendicular magnetic field on magnetic vortex are
well established Ivanov2002 ; Loubens2009 : the in-plane spins are tilted
towards the applied field producing hereby a decrease 111At $H_{0}=6$ kOe, the
spins are tilted out-of-plane by about 20∘ generating a 7% decrease of the in-
plane component of the magnetization outside the vortex core. of the in-plane
component of the magnetization outside the vortex core (cone state Ivanov2002
).

We then proceed to the measurement of the variation of the excitation spectrum
of the gyrotropic mode when the tip is scanned by $\pm 0.85\mu$m along the
$x$-direction by steps of 50 nm (see FIG.2). The scan height is $0.9\pm
0.05\mu$m 222Although the use of piezo-actuators allows ultra-precise
displacement of the micro-cantilever, the value of $\delta_{z}$ has inherently
some uncertainty as it corresponds to the free axis of the cantilever. above
the nanodot or $\delta_{z}=3.0\pm 0.15$ in reduced units of $R$ (hereafter all
spatial displacements are expressed in units of the dot radius $R$.) Placing
the tip at the origin ($\delta_{x}=0$ or on the axis of the disk), attracts
the vortex core at the center of the nanodot. From there, lateral displacement
of the magnetic tip produces a vector shift ${\bm{X}}_{0}$ of the vortex core
equilibrium position from the dot center. The process is driven by the growth
of the in-plane domain parallel to the in-plane component of the tip stray
field. We observe in FIG.2 that the eigen-frequency _increases_ (blue shift)
upon increasing $\left|{\bm{X}}_{0}\right|$. Noting that the frequency shift
is symmetric and isotropic (the signature that the intrinsic potential is
being probed) we find that the vortex core dynamics can be tuned on a relative
large range ($\sim$ 10%), hereby demonstrating that the magnetostatic
potential must be anharmonic, since a purely parabolic shape would have lead
to a frequency independent behavior on ${\bm{X}}_{0}$.

Analysis of the amplitude of the MRFM signal gives a hint on the amount of
displacement $\left|\bm{X}_{0}\right|$ achieved during a scan. The MRFM signal
corresponds to the difference of vertical force $\Delta F_{z}$ acting on the
cantilever when the vortex motion is excited. Defining $\bm{m}=\bm{M}/M_{s}$
the reduced magnetization vector, the gyromotion produces a diminution of the
spontaneous magnetization along the local equilibrium direction $\Delta
m_{i}=\frac{1}{2}\left|\partial_{X}m_{i}+j\,\partial_{Y}m_{i}\right|^{2}_{\bm{X}=\bm{X}_{0}}$,
that mostly occurs outside the core region Guslienko2008a . The generated
force can be then calculated from the reaction force $-\Delta
F_{z}=VM_{s}\langle g_{zi}\Delta m_{i}\rangle$ acting on the nanodot due to
the gradient tensor of the tip:
${\widehat{\bm{g}}}=\nabla\bm{H}^{\text{tip}}$. Cartesian tensor notation is
used here, with repeated indices being assumed summed. The chevron bracket
indicates that the enclosed quantity is averaged over the volume $V$ of the
nanodot. The red-blue colormap in FIG.2 codes the amplitude MRFM signal: red
(blue) means attractive (repulsive) force. Small arrows at
$\delta_{x}\approx\pm 1.7$ indicate the compensation point of the force: where
the change of sign occurs. This distance is about twice smaller than the one
required to change the sign of the force in the saturated state (contribution
dominated by $\langle g_{zz}\Delta m_{z}\rangle$). The difference is
interpreted as due to the translation of the core position transversally to
the tip position (along the $y$-axis). The growth of this domain generates a
repulsive vertical force on the tip trough the cross gradient term $\langle
g_{zx}\Delta m_{x}\rangle$. We shall thus use the position of the compensation
point to calibrate precisely the amplitude of the displacement of the vortex
core.

Figure 2: (Color online) Density plot showing the experimentally measured
variation of the eigen-frequency of the gyrotropic mode upon a lateral
displacement, $\delta_{x}$, of the MRFM tip at fixed $\delta_{y}=0$ and
nominal $\delta_{z}=3.0$. A red-blue colormap shows the sign of the force
acting on the cantilever. The arrows indicate where the change of polarity
occurs.

Our next step is to develop an analytical framework allowing the extraction of
$\lambda$ from this variational study. The first stage of this analysis is to
calculate the equilibrium position $\bm{X}_{0}=(X_{0},Y_{0})$ (here $X_{0}$
and $Y_{0}$ are the two in-plane cartesian coordinates) by minimizing the
total energy $W=W^{\text{(M)}}+W^{\text{(H)}}$, the sum of $W^{\text{(M)}}$,
the magnetostatic self-energy of the vortex ground state which confines the
vortex core to the center of the nanodot and
$W^{\text{(H)}}=-VM_{s}\langle\bm{m}\cdot\bm{H}\rangle$, the Zeeman energy
which represents the interaction with the external magnetic field
$\bm{H}=\bm{H}_{0}+\bm{H}^{\text{tip}}$ and is responsible for the
displacement $\bm{X}_{0}$. The magnetic configuration inside the nanodot
$m_{x}+j\,m_{y}=2\mathit{w}/(1+\mathit{w}\mathit{w}^{\ast})$ is conveniently
described by a conformal mapping of the complex variable $\mathit{w}$ ($\ast$
indicating the complex conjugate), which is a piecewise function of the
complex position $\text{\raisebox{1.29167pt}{\footnotesize\cursive
z}}=(x+j\,y)/R$ with
$\mathit{w}=f(\text{\raisebox{1.29167pt}{\footnotesize\cursive
z}})/|f(\text{\raisebox{1.29167pt}{\footnotesize\cursive z}})|$ outside the
vortex core region and
$\mathit{w}=f(\text{\raisebox{1.29167pt}{\footnotesize\cursive z}})$ inside.
The function $f$ captures the texture of the spatial configuration. To
calculate the new equilibrium position $\mathcal{Z}_{0}=(X_{0}+j\,Y_{0})$, it
is appropriate to describe the dot magnetization, $\bm{m}$, by the rigid
vortex model (RVM) Guslienko2001 , written as
$f(\text{\raisebox{1.29167pt}{\footnotesize\cursive z}})=\pm
j\,(\text{\raisebox{1.29167pt}{\footnotesize\cursive
z}}-\mathcal{Z}_{0})/r_{c}$. Here $r_{c}=R_{c}/R$ denotes the core radius
$R_{c}$ in reduced unit of $R$ and the $\pm$ sign depends on the chirality of
the vortex. The static displacement $\bm{X}_{0}$ is then obtained by
minimizing the total energy 333All energies can be approximately calculated by
running the integrand solely on the nanodot volume outside the vortex core
region. with the analytical expression of $W^{\text{(M)}}$ obtained by the
RVM. In the RVM, the magnetostatic energy is generated by the surface magnetic
charges $\sigma$ located at the circumference of the disk (the volume charges
$\nabla\cdot\bm{M}$ are absent). The confinement potential follows from the
integral $W^{\text{(M)}}=\frac{1}{2}\int d\phi\int
d\phi^{\prime}\sigma(\phi)\sigma(\phi^{\prime})/\sqrt{2(1-\cos(\phi-\phi^{\prime}))}$
where the integration is taken over the disk periphery and $\sigma$ is given
by:

$\sigma(\phi)=+M_{s}\frac{-|\bm{X}_{0}|\sin(\phi-\phi_{0})}{\sqrt{1-2|\bm{X}_{0}|\cos(\phi-\phi_{0})+|\bm{X}_{0}|^{2}}}\,,$
(2)

$\phi_{0}$ is the azimuthal direction of the vortex equilibrium position
measured from the averaged in-plane bias field direction (here $x$-axis).

The implicit trajectory of $\bm{X}_{0}$ is shown in FIG.3a for three different
heights $\delta_{z}$ around the nominal value. As expected the in-plane
components of the tip magnetic field displace the vortex core mainly along the
$y$-axis. The displacement along $x$-axis is approximately twice smaller. The
resulting displacement distance $|\bm{X}_{0}|$ as a function of $\delta_{x}$
is shown in FIG.3b. We use here a skewed scale on the abcisse to show the
behavior when $\delta_{x}\gg 1$. We have also calculated the corresponding
dipolar force produced on the tip. The result is coded in the colormap using
the same convention as in FIG.2. We have placed small arrows at the
compensation points. Since decreasing the scan height increases the amplitude
of $|\bm{X}_{0}|$, we find that the position of the arrows sensitively depends
on $\delta_{z}$. Varying $\delta_{z}$ in the experimental error bars
$[2.6,3.0]$ displaces the compensation point by $\pm 0.3\cdot R$ (or $\pm 100$
nm) around the mean value $\delta_{x}=1.8$, in agreement with the experimental
data. We shall use this marker to evaluate the uncertainty window of
$|\bm{X}_{0}|$ in our experiment.

The second stage of this analysis is to perform a linearization of the vortex
equation of motion to a cyclic excitation field. The instantaneous response
$\bm{X}=(X,Y)$ is decomposed into the static component $\bm{X}_{0}$,
calculated previously, and a dynamic component $\bm{\xi}=\bm{X}-\bm{X}_{0}$
representing the small oscillating deviation of the vortex core position from
its equilibrium 444We estimate that $\xi=0.07$ in our experiment, where the
microwave field strength is $h_{\text{rf}}=0.6$ Oe.. In the dynamical case,
the dipolar pinning imposes a precession node at the dot circumference
Guslienko2008a . It implies that the dynamical magnetization comes from the
variation, $\partial_{X}{\bm{m}}+j\,\partial_{Y}{\bm{m}}$, of a magnetic
configuration that has no radial component at the dot border. Therefore, to
calculate the frequency of the small dynamic vortex displacement $\bm{\xi}$,
it is appropriate to use the surface charges free model or two vortex ansatz
(TVA) written as $f(\text{\raisebox{1.29167pt}{\footnotesize\cursive z}})=\mp
j\,\frac{1}{r_{c}}(\text{\raisebox{1.29167pt}{\footnotesize\cursive
z}}-\mathcal{Z})(\text{\raisebox{1.29167pt}{\footnotesize\cursive
z}}\mathcal{Z}^{\ast}-1)/(1+|\mathcal{Z}|^{2})$ Guslienko2002  with
$\mathcal{Z}=(X+j\,Y)$. In our notation, the dampingless Thiele equation
simply writes $\bm{G}\times\dot{\bm{\xi}}=\partial W/\partial\bm{\xi}$, where
$W$ is the total energy and $\left|\bm{G}\right|=2\pi M_{s}t/\gamma$ is the
gyrovector Guslienko2008  (the dot is the short hand notation for the time
derivative and $\gamma$ is the gyromagnetic ratio). Linearization around
$\bm{X}_{0}$ yields the gyrotropic angular frequency

$\omega^{2}=\frac{K_{xx}K_{yy}-K_{xy}^{2}}{G^{2}}\,,\ \text{with}\
K_{ij}\equiv\left.\frac{\partial^{2}W}{\partial{\xi}_{i}\partial{\xi}_{j}}\right|_{\bm{X}=\bm{X}_{0}}$
(3)

being the stiffness of the vortex core to small displacements in both the $i$
and $j$ directions. Distinction between different cartesian directions is
necessary once the trajectory becomes elliptical. This is precisely, what
occurs when $\bm{X}_{0}\gg\bm{\xi}$: the amplitude of the $\xi$-component
along $\bm{X}_{0}$ differs from the amplitude of the $\xi$-component
perpendicular to $\bm{X}_{0}$ (short axis of the ellips is along the radial
direction). The degree of ellipticity is determimed by the anharmonic
contribution $\lambda\left|\bm{X}_{0}\right|^{2}$. This is in contrast to the
opposite limit $\bm{X}_{0}\ll\bm{\xi}$, where the trajectory corresponds to a
large amplitude circular vortex core motion around the nanodot center
Dussaux2012 .

To calculate the different tensor elements of the stiffness
$K_{ij}=K^{\text{(M)}}_{ij}+K^{\text{(H)}}_{ij}$, one must decompose it in two
contributions corresponding respectively to the magnetostatic and Zeeman
energies. The first order value of the TVA magnetostatic stiffness, $\kappa$,
has been already expressed analytically Guslienko2006 . The analytical
expression of the anharmonic correction is obtained by inserting Eq.(1) in
Eq.(3) and it leads to a simplified expression
$K^{\text{(M)}}_{ij}=\kappa\left.\left(\delta_{ij}+\lambda|\bm{X}|^{2}\delta_{ij}+2\lambda
X_{i}X_{j}\right)\right|_{\bm{X}=\bm{X}_{0}}$. It turns out that the Zeeman
stiffness can be neglected. Indeed, it can be shown that the tip stray field
produces no Zeeman stiffness along the diagonal elements
($K^{\text{(H)}}_{ii}=0$). Only the cross-terms $K^{\text{(H)}}_{xy}\neq 0$
are non-vanishing but they represent a negligible correction ($<$ 3%). We thus
find that at $\bm{X}_{0}=0$ and $H_{z}=0$, Eq.(3) simplifies to the well known
expression $\omega(0,0)={\kappa}/{G}$ Guslienko2002 . At $\bm{X}_{0}=0$ and
$H_{z}\neq 0$, the stiffness of the magnetostatic potential is renormalized by
the in-plane magnetization projection of the cone state and one obtains
${\omega(0,H_{z})}/\omega(0,0)=1+H_{z}/(4\pi M_{s})$ Loubens2009 . In the
general case $\bm{X}_{0}\neq 0$ and $H_{z}\neq 0$, the relative frequency
shift reduces to the following analytical expression:

$\frac{\omega(\bm{X_{0}},H_{z})}{\omega(0,H_{z})}=1+2\lambda\left|\bm{X_{0}}\right|^{2}+{\cal{O}}(|{\bm{X}}|^{4})\,.$
(4)

Notice that the prefactor of 2 multiplying $\lambda$ is specific to the limit
$\bm{X}_{0}\gg\bm{\xi}$.

The next step is to plot in FIG.4a the experimental data extracted from FIG.2,
renormalized by the predicted dependence of ${\omega(0,H_{z})}$, as a function
of the calculated $|\bm{X}_{0}|$ during a lateral scan of the tip at fixed
$\delta_{z}=2.8$. Fitting the data of FIG.4a with a parabola (solid line)
yields an average curvature $\lambda=0.5$. We have plotted in FIG.4b the
experimentally measured relative frequency normalized by $\omega(0,H_{0})$.
The latter quantity is inferred experimentally by studying the decay of
$\omega$ upon increasing $\delta_{z}$, while keeping the tip on the symmetry
axis ($\delta_{x}=\delta_{y}=0$): a fit of the decay behavior yields the
asymptotic value $\omega(0,H_{0})$. In FIG.4b, the data point are colored
according to the colormap associated with the amplitude of the force. For
comparison, we have also plotted the predicted variation of $\omega$ by Eq.(4)
as a function of $\delta_{x}$ for two values of $\lambda$. Setting $\lambda=0$
in Eq.(4), would have produced the usual bell-shaped curve Pigeau2012 , which
corresponds to a diminution of $\omega(0,H_{z})$ when the tip moves away from
the nanodot axis. The behavior for $\lambda=0.5$ is in excellent agreement
with the experimental data, both in the amplitude of the NL frequency shift
and in the position of the compensation point of the force.

We have then repeated the analysis by varying $\delta_{z}$ in the experimental
error bar range: $\pm 0.2$ around the nominal value. Fit of the data by a
parabola would lead to larger (smaller) values of $\lambda$ depending if the
amplitude of the shift decreases (increases). This procedure yields an
uncertainty window of 30% for the determination of $\lambda$, shown as a
shaded area in FIG.4a. Our fitting analysis did not account for higher order
corrections in Eq.(1). In FIG.4a, the curvature increases with the
displacement distance. Inclusion in the fit of terms in
$\left|\bm{X}\right|^{4}$ would have decrease the value of $\lambda$ by about
one standard deviation. As an additional check, we have performed a simulation
of the expected $\omega(X_{0},H_{z})$ for our nanodot using a mesh-size of 2.3
nm and a GPU-accelerated micromagnetic code Vansteenkiste2011a . The result is
shown as crosses in FIG.4a, demonstrating that our determination of $\lambda$
is in quantitative agreement with numerical simulations. It is also in
agreement with the result obtained by Dussaux et al. from micromagnetic
simulations performed in the limit $\xi\gg X_{0}$ on a thinner dot with
approximately the same radius Dussaux2012 .

Figure 3: (Color online) a) Trajectory of the vortex core
${\bm{X}}_{0}=(X_{0},Y_{0})$ during an implicit lateral scan of the tip
$\delta_{x}\in[-8,+8]$ at 3 different heights, $\delta_{z}$. b) Norm of the
displacement vector, $\left|\bm{X}_{0}\right|$, as a function of the tip
position, $\delta_{x}$.

In summary, using an MRFM, we have measured quantitatively the anharmonicity
coefficient $\lambda=+0.5\pm 0.15$ produced by the depolarisation field of a
vortex in a planar nanodot 555A recent preprint posted during the review
process reports a new analytical expression for $\lambda$ compatible with our
measurement Metlov2013 . From a fundamental perspective, it is interesting to
note that the obtained value (dipole-dominated) is about twice smaller than
the $\lambda=1$ predicted by the local easy-plane model Ivanov1998 . Further
work is required to check if the value is independent of the out-of-plane
external magnetic field $\bm{H}_{0}$ or the dot aspect ratio, $R/t$, in
particular around the line of the vortex state stability $\kappa(R,t)=0$
Guslienko2008 . Finally, we mention that foldover experiments performed on the
same nanodot at $\bm{X}_{0}=0$ produce a red shift of the gyrotropic frequency
(regime $\bm{X}_{0}\ll\xi$), which is opposite with respect of the sign of
$\lambda>0$. This finding suggests that the NL frequency shift observed in the
foldover of the resonant curve is not dominated by the anharmonicity of the
magnetostatic potential, but perhaps by the NL damping Pigeau2011 .

Figure 4: (Color online) a) Plot of the relative variation of eigen-frequency
of the gyrotropic mode as a function of $|\bm{X}_{0}|$. A fit by a parabola
yields $\lambda=0.5\pm 0.15$. The cross are the results of micromagnetic
simulations b) Relative variation of the eigen-frequency as a function of
$\delta_{x}$. The lines show the analytically predicted behavior for a
vanishing (0) and finite (0.5) anharmonic coefficient.

###### Acknowledgements.

This research was partly supported by the EU grant MOSAIC (ICT-FP7-317950),
the French ANR Grant MARVEL (ANR-2010-JCJC-0410-01), the Spanish MEC Grants
PIB2010US-00153 and FIS2010-20979-C02-01. S.S., K.G. and O.S. acknowledge
support from the Marie Curie grant AtomicFMR (IEF-301656), from the IKERBASQUE
and from the UPV/EHU, respectively.

## References

  * (1) A. Slavin and V. Tiberkevich, Ieee Transactions On Magnetics 45, 1875 (Apr. 2009)
  * (2) S. M. Mohseni, S. R. Sani, J. Persson et al., , Science 339, 1295 (Mar. 2013)
  * (3) V. S. Pribiag, I. N. Krivorotov, G. D. Fuchs, P. M. Braganca, O. Ozatay, J. C. Sankey, D. C. Ralph, and R. A. Buhrman, Nature Physics 3, 498 (Jul. 2007),
  * (4) N. Locatelli, V. V. Naletov, J. Grollier, G. de Loubens, V. Cros, C. Deranlot, C. Ulysse, G. Faini, O. Klein, and A. Fert, Applied Physics Letters 98, 062501 (Feb. 2011)
  * (5) B. Pigeau, G. de Loubens, O. Klein, A. Riegler, F. Lochner, G. Schmidt, L. W. Molenkamp, V. S. Tiberkevich, and A. N. Slavin, Applied Physics Letters 96, 132506 (Mar. 2010)
  * (6) A. Dussaux, A. V. Khvalkovskiy, P. Bortolotti, J. Grollier, V. Cros, and A. Fert, Physical Review B 86, 014402 (Jul. 2012)
  * (7) V. Novosad, F. Y. Fradin, P. E. Roy, K. S. Buchanan, K. Y. Guslienko, and S. D. Bader, Phys. Rev. B 72, 024455 (Jul 2005),
  * (8) A. Drews, B. Krüger, G. Selke, T. Kamionka, A. Vogel, M. Martens, U. Merkt, D. Möller, and G. Meier, Phys. Rev. B 85, 144417 (Apr 2012),
  * (9) B. Pigeau, G. de Loubens, O. Klein, A. Riegler, F. Lochner, G. Schmidt, and L. W. Molenkamp, Nature Physics 7, 26 (Jan. 2011)
  * (10) Y. Gaididei, V. P. Kravchuk, and D. D. Sheka, International Journal of Quantum Chemistry 110, 83 (Jan. 2010)
  * (11) T. Y. Chen, M. J. Erickson, P. A. Crowell, and C. Leighton, Phys. Rev. Lett. 109, 097202 (Aug 2012),
  * (12) J. A. J. Burgess, A. E. Fraser, F. F. Sani, D. Vick, B. D. Hauer, J. P. Davis, and M. R. Freeman, Science 339, 1051 (2013),
  * (13) O. Klein, G. de Loubens, V. V. Naletov, F. Boust, T. Guillet, H. Hurdequint, A. Leksikov, A. N. Slavin, V. S. Tiberkevich, and N. Vukadinovic, Physical Review B 78, 144410 (Oct. 2008)
  * (14) K. Mitsuzuka, D. Lacour, M. Hehn, S. Andrieu, and F. Montaigne, Applied Physics Letters 100, 192406 (May 2012)
  * (15) B. Pigeau, C. Hahn, G. de Loubens, V. V. Naletov, O. Klein, K. Mitsuzuka, D. Lacour, M. Hehn, S. Andrieu, and F. Montaigne, Phys. Rev. Lett. 109, 247602 (Dec 2012),
  * (16) B. A. Ivanov and G. M. Wysin, Phys. Rev. B 65, 134434 (Mar 2002),
  * (17) G. de Loubens, A. Riegler, B. Pigeau, F. Lochner, F. Boust, K. Y. Guslienko, H. Hurdequint, L. W. Molenkamp, G. Schmidt, A. N. Slavin, V. S. Tiberkevich, N. Vukadinovic, and O. Klein, Physical Review Letters 102, 177602 (May 2009)
  * (18) At $H_{0}=6$ kOe, the spins are tilted out-of-plane by about 20∘ generating a 7% decrease of the in-plane component of the magnetization outside the vortex core.
  * (19) Although the use of piezo-actuators allows ultra-precise displacement of the micro-cantilever, the value of $\delta_{z}$ has inherently some uncertainty as it corresponds to the free axis of the cantilever.
  * (20) K. Y. Guslienko, A. N. Slavin, V. Tiberkevich, and S.-K. Kim, Phys. Rev. Lett. 101, 247203 (Dec 2008),
  * (21) K. Y. Guslienko, V. Novosad, Y. Otani, H. Shima, and K. Fukamichi, Phys. Rev. B 65, 024414 (Dec 2001),
  * (22) All energies can be approximately calculated by running the integrand solely on the nanodot volume outside the vortex core region.
  * (23) We estimate that $\xi=0.07$ in our experiment, where the microwave field strength is $h_{\text{rf}}=0.6$ Oe.
  * (24) K. Y. Guslienko, B. A. Ivanov, V. Novosad, Y. Otani, H. Shima, and K. Fukamichi, Journal of Applied Physics 91, 8037 (May 2002)
  * (25) K. Y. Guslienko, Journal of Nanoscience and Nanotechnology 8, 2745 (Jun. 2008)
  * (26) K. Y. Guslienko, X. F. Han, D. J. Keavney, R. Divan, and S. D. Bader, Physical Review Letters 96, 067205 (Feb. 2006)
  * (27) A. Vansteenkiste and B. Van de Wiele, Journal of Magnetism and Magnetic Materials 323, 2585 (Nov. 2011)
  * (28) A preprint posted on arXiv during the review process reports a new analytical expression for $\lambda$ compatible with our measurement Metlov2013
  * (29) B. A. Ivanov, H. J. Schnitzer, F. G. Mertens, and G. M. Wysin, Phys. Rev. B 58, 8464 (Oct 1998),
  * (30) K. L. Metlov, ArXiv e-prints(Aug. 2013), arXiv:1308.0240 [cond-mat.mes-hall]