11institutetext: Department of Computer Science and Engineering, University of
Moratuwa, Moratuwa, Sri Lanka
11email: {dinukamanohara.15,ashkar.15,isurumaldeniya.15<EMAIL_ADDRESS>22institutetext: Lawrence Berkeley National Lab, California, USA
22email<EMAIL_ADDRESS>33institutetext: University of Virginia, USA
33email<EMAIL_ADDRESS>
# Parallel Algorithms for Densest Subgraph Discovery Using Shared Memory Model

B.D.M. De Zoysa 11 Y.A.M.M.A. Ali 11 M.D.I. Maduranga 11 Indika Perera 11
Saliya Ekanayake Anil Vullikanti 2233

###### Abstract

The problem of finding dense components of a graph is a widely explored area
in data analysis, with diverse applications in fields and branches of study
including community mining, spam detection, computer security and
bioinformatics. This research project explores previously available algorithms
in order to study them and identify potential modifications that could result
in an improved version with considerable performance and efficiency leap.
Furthermore, efforts were also steered towards devising a novel algorithm for
the problem of densest subgraph discovery. This paper presents an improved
implementation of a widely used densest subgraph discovery algorithm and a
novel parallel algorithm which produces better results than a 2-approximation.

###### Keywords:

Densest subgraph K-core Task parallelism OpenMP Graph algorithms Data mining
Shared memory algorithms.

##  1 Introduction

###  1.1 Background

Discovering the dense components of a graph has become an important primitive
in data analysis over the past few decades. Numerous algorithms and techniques
had been developed and they have evolved throughout the period in order to
ascertain the adequacy of methods to determine the dense components of
different types of graphs namely, undirected, directed and weighted graphs
with minimal requirement of resources. One such approach to discover dense
components of a graph is to find the densest subgraph of a given graph.

By observing the applications that exploit the methods of finding densest
subgraphs, one can easily acknowledge the substantial impact these methods
have over modern technology due to their supplementary nature with well-known
fields of research such as bioinformatics, data science and network analysis.
Even though researchers have developed a variety of methods to perform the
task of finding the densest subgraph of a graph, we could discern that there
are areas yet to be explored and implementations that could be improved such
that resource and time consumption can be mitigated further. This research
focuses on improved methods for undirected graphs which can be modified for
other types of graphs.

###  1.2 Preliminaries

###### Definition 1

Edge density $\rho\textit{(G)}$ of a graph G(V,E) is defined as
$\rho\textit{(G)}=\frac{|E|}{|V|}$. V is the set of vertices and E is the set
of edges induced by V.

###### Definition 2

The k-core of a graph is the largest subgraph in which every vertex has degree
at least k.

###### Definition 3

For $\alpha\geq\textit{1}$, an algorithm is said to obtain an
$\alpha\textit{-approximation}$ to the undirected densest subgraph if it
outputs a set $\textit{\~{S}}\subseteq\textit{V}$ such that
$\rho\textit{(}\textit{\~{S}}\textit{)}\geq\frac{\rho^{\textit{*}}\textit{(G)}}{\alpha}$
[9]. $\rho^{\textit{*}}\textit{(G)}$ is the density of the densest subgraph of
graph G(V,E).

###  1.3 Contributions

Our main contributions are,

  1. 1.

Highly scalable shared memory implementation of densest subgraph approximation
algorithm, Algorithm 1 in Bahmani et. al [9]. (In this paper we call this
algorithm as Bahmani’s algorithm.)

  2. 2.

A new highly scalable parallel algorithm which calculates more accurate
density values than 2-approximation algorithms.

##  2 Related Work

Dense subgraph discovery in graphs has extensively been studied[25, 40, 13,
50]. The algorithms for finding the dense components in a graph can be
categorized as exact enumeration, fast heuristics, and bounded approximation
algorithms [3].

The simplest and the most common way to define the density of an undirected
graph is the edge density. The densest subgraph of a given graph based on the
edge density is the subgraph that has the highest edge density than all the
other possible subgraphs. Goldberg et al. presented a maximum-flow based exact
algorithm to calculate that densest subgraph [26]. Finding densest subgraphs
with size restrictions (i.e., densest subgraphs which consist of a specified
number of vertices) is a variant of this edge density densest subgraph problem
and identified as NP-hard [8, 7]. Another variant is the optimal quasi-clique
problem [50]. This is also NP-hard [11]. Solutions for finding the top-k
locally densest subgraphs were presented by Qin et al. [39]. Densest subgraph
problem induced by evolving graphs is studied in [18]. Graph decomposition
based on the edge density is explored in [48, 17]. The density for directed
graphs is expatiated in [34] and the densest subgraph problem on directed
graphs is studied in [12].

Approximation algorithms and algorithms based on heuristics to discover the
densest subgraphs come into the picture when the exact algorithms fail to
perform efficiently on large graphs. A greedy 2-approximation algorithm is
presented in [12]. Bahmani et al. [9] presented ($2+2\epsilon$)-approximation
algorithm ($\epsilon>0$) and this is the algorithm that is being parallelized
in this project. There exists an approximation algorithm to discover the
densest subgraph in directed graphs [35]. The h-clique densest subgraph
problem is defined and studied by Tsourakakis et al. [49, 37]. Top-k local
triangle-densest subgraph discovery is based on the 3-clique and has been
studied in [41]. Fang et al. [56] proposed exact and approximation algorithms
for h-clique ($\textit{h}\geq 2$) densest subgraph problem based on core
decomposition. Two prime issues of previous algorithms namely, resulting in
isolated graphs and neglecting of important features of natural graphs are
addressed by Bo Wu and Haiying Shen [4]. They present a heuristic algorithm
for massive undirected graphs and an exact algorithm for big data in order to
handle the aforementioned issues.

Rozenshtein et al. [5] study the problem of dynamic dense subgraph whose edge
occur in a short time interval or equivalently, are temporary compact. This
problem is NP-hard [5]. McGregor et al. [6] discuss the problem of
approximating the densest subgraph in the dynamic graph stream model. They
also present the first algorithm for approximating the density of the densest
subgraph up to ($1+\epsilon$) factor in the dynamic graph stream model.

K-core [10, 36, 38, 21, 20, 19, 23, 24, 51, 14], k-truss [15, 30, 53, 32, 31],
k-(r, s) nucleus [44, 42, 45, 43], k-clique [16, 27], k-edge connected
components [28, 29] and k-plexes [46] can also be identified as the other
formations of dense graph component models. Tatti et al. [47] proved that the
largest core of a graph is a 2-approximation to the densest subgraph.

##  3 Methodology

###  3.1 Parallelization of Bahmani’s Algorithm

####  3.1.1 Introduction

Bahmani et al. presents a greedy approximation algorithm for finding densest
subgraph of an undirected graph. For any $\epsilon>0$, their algorithm makes
$\textit{O}(\log_{1+\epsilon}\textit{n})$ passes over the input data and finds
a subgraph of which the density is guaranteed to be within a factor of
($2+2\epsilon$) of the optimum. As the authors claim, the algorithm is
inherently parallelizable with minimal data dependencies between parallelized
components.

As provided in their description, the algorithm takes an undirected graph
G(V,E) as input. It proceeds in passes while removing a constant fraction of
the remaining nodes in each pass. They have also proved that one of such
intermediate subgraphs forms a ($2+2\epsilon$)-approximation to the densest
subgraph. Furthermore, they note that, via flows or linear programming densest
subgraph problem in undirected graph can be solved in polynomial time.
Conversely, they have observed that flow and linear programming techniques
scale poorly to internet sized graphs. However, despite the worst-case
examples, they have established the evidences for their claim of their
algorithms having the capability of yielding near-optimal solutions and being
simpler and highly efficient than flow or linear programming-based techniques.

####  3.1.2 Approach

The approach can be broken down mainly into two parts as described below:

##### Exploitation of task parallelism:

Tasks were introduced in version 3.0 of OpenMP. Before the availability of
tasks, it was possible to execute loops of which the length was known at run
time. In addition, there was also an imposed limit in the number of parallel
sections. Most of the times this did not work well with even some well- known
problems. i.e., even though it was possible to come with a feasible
workaround, the solution was often unattractive. Until the introduction of
tasks, it was not possible to efficiently and easily implement several types
of parallelism. With tasks the programmer has to specify where the tasks are;
then the OpenMP runtime system will generate a new task when a thread
encounters a task construct. It is up to the OpenMP runtime system to ‘decide’
the moment of execution of the task whereas it can be delayed or immediate.
task synchronization is used to enforce the completion of a task. Real world
graphs are often highly unbalanced. Due to the nature of task creation,
execution and its properties, it is plain that OpenMP tasks are highly
applicable to the problem at hand. This is because if we use other patterns
such as loop parallelism or work sharing constructs, the parallelized
components that have comparatively lower amount of work to do have to be idle
for most of the execution time. On the contrary, if the pattern used is task
parallelism then the aforementioned occurrence would not be the case as tasks
are not always unique to threads; consequently, a task can take over another
component in the queue as soon as its work is done as opposed to being idle.
In a considerable number of cases this could reduce the overhead of thread
creation while helping execute the program quicker than other cases. The
component above the dashed line representing barrier synchronization in figure
3 diagrammatically depicts this setting for a program with four tasks.

##### Partitioning the algorithm into two major parts:

First major part of the modified algorithm can execute without any
synchronizations between the parallelized components. When it comes to the
second major part, a synchronization is required due to the updating of
degrees of vertices. Assuming that the threshold for removing vertices is
between 3 and 4 in a given iteration, we can see in figure 3 that the grey
shaded vertices should fail as their degree is lower than the threshold. If
those vertices belong to two different tasks, then the neighboring red
coloured critical vertex is vulnerable to race condition. This establishes the
requirement for a synchronization.

####  3.1.3 Implementation

##### Data structures:

Despite relatively high memory usage we have used hash tables for storing the
graph as our main focus was to be able to test the run-times for graphs of all
types of vertex distribution. i.e., by doing this, we were able to support the
graphs regardless of whether a graph of n vertices contains vertices from $0$
to $(\textit{n}-1)$ or not. For this purpose, we have used the unordered map
which is the C++ implementation of hash tables. We define a structure named
super map which concretely represents the whole graph. Each key of the hash
table represents a vertex of the graph. The value for each key is another hash
table unordered map which has the neighboring vertices and their degrees as
keys and values respectively.

##### Parallelization:

  * •

Task ID creation: As mentioned previously, tasks and threads does not always
have one to one relationship. Therefore, as our program requires an ID unique
to each task, the ID creation process is done by defining and initializing a
shared variable ID and then atomically assigning the variable to the thread
private ID immediately after each individual task commences. This is followed
by increment of the shared variable which can then be assigned to the next
task. The id creation process is identical in both the major parts of
P-Bahmani.

  * •

Parallelism in the first part of P-Bahmani: This part does not contain any
inter-task data dependencies. Moreover, for the purposes of time reduction in
execution of the program, we have exercised a strategy where the algorithm
avoids the testing of already visited vertices by removing those vertices.
Here, each parallelized component iterates over the set of locally assigned
vertices and removes the failed vertices from the locally assigned active set
of vertices. With this strategy, the algorithm does not need to go over the
failed vertices in the upcoming iteration. This cannot be done with work
sharing constructs or basic type of loop parallelism where the Internal
Control Variables are bound to a thread because an increment clause of OpenMP
for loop must only perform simple addition/subtraction on the loop variable.
However, the kind of increment our program consists in order to perform the
removal of elements and the iteration over the same set is case dependent. The
for loop is defined only with the pointer initiation and the condition. But
the method of the pointer assignment must be defined on the fly depending on
whether the current vertex failed or not. This can only be realized with
tasks.

  * •

Parallelism in the second part of P-Bahmani: Each of the parallelized unit of
executions contains a unique set of failed vertices resulted from the first
part of the algorithm. When removing edges, we have to decrease the degrees of
sources as well as the destinations. Therefore, while iterating over the
unique set of task-local failed vertices, it is necessary to iterate over the
set of neighbors of the current vertex since we have to decrease the degree of
source vertex via the direction source$\rightarrow$destination and the degree
of destination vertex via the direction destination$\rightarrow$source. Since
the source vertex is unique to each task and both the aforementioned
directions consists the source, we can vouch that the assignment of variables
at the addresses of sources won’t encounter any race conditions. Hence, the
task local edge reduction variable can safely be used for decreasing the
number of edges. Conversely, it is obvious that the destination vertex is
vulnerable to race conditions which dictates us to use a critical section when
updating its degree. We have tackled this through the C++’s built-in function
for atomic reductions. The second part of the algorithm has the form shown in
figure 6.

Input : $\textit{G(V,E)},\epsilon>0$

1 $\tilde{G}\leftarrow G$;

2 $n_{e},n_{v}\leftarrow 0;$$\hfill\triangleright Task-local\hskip
2.84544ptvariables\hskip 2.84544ptof\hskip 2.84544pttype\hskip
2.84544pt\textit{int}$

3$active;failed;$$\hfill\mathrm{\triangleright Task-local\hskip
2.84544ptvariables\hskip 2.84544ptof\hskip 2.84544pttype\hskip
2.84544pt\textit{set}}$

4$\rho,\tilde{\rho}\leftarrow\rho_{0};$$\hfill\mathrm{\triangleright
Global\hskip 2.84544ptvariables\hskip 2.84544ptof\hskip 2.84544pttype\hskip
2.84544pt\textit{set}}$

5while  _$G\rightarrow|V| >0$_ do

6 parallel for  _$\nu\in active$_ do

7 if  _$\nu.deg\leq 2(1+\epsilon)\rho$_ then

8 $failed.insert(\nu)$;

9 $n_{v}\leftarrow n_{v}-1$;

10 $active.remove(\nu)$;

11

12 else

13 $pass$;

14

15 end if

16

17 end for

18

19 barrier synchronization

20 parallel for  _$\nu\in failed$_ do

21 for  _$u\in\nu.neighbours$_ do

22 $n_{e}\leftarrow n_{e}-u.local\\_deg$;

23 $\textbf{atomicSub}(u.deg,u.local\\_deg)$;

24

25 end for

26

27 end for

28

29
$\textbf{reduce}(G\rightarrow|V|,n_{v});\textbf{reduce}(G\rightarrow|E|,n_{e})$;

30 $n_{v},n_{e}\leftarrow 0$;

31 $failed.clear()$;

32

33 if  _$\rho >\tilde{\rho}$_ then

34 $\tilde{G}\leftarrow G$;

35

36 else

37 $pass$;

38

39 end if

40

41 end while

42$\textbf{Return: }\tilde{G}$;

Algorithm 1 P-Bahmani

Figure 1: Tasks, parts and barrier

Figure 2: Critical node

Figure 3: Graph data structure

###  3.2 Core-based Dense Subgraphs

####  3.2.1 Approach

We previously described how the densest subgraph is correlated to the densest
k-core (${\textit{S}_{\textit{k}}}^{*}$). Identification of the densest k-core
provides us with a 2-approximation to the densest subgraph. This connotes that
if we attempt to improve the density of the densest k-core, we will end up
with a better approximation to the densest subgraph.

The notion of improving the density of ${\textit{S}_{\textit{k}}}^{*}$ can be
realized by adding new vertices from the set of vertices that does not belong
to ${\textit{S}_{\textit{k}}}^{*}$. We show that the density of a subgraph
$\tilde{S}$ can be increased by adding appropriately chosen vertices where
$\tilde{S}$ is the subgraph which accommodates the newly added vertices,
starting from ${\textit{S}_{\textit{k}}}^{*}$. Therefore, we introduce a
method of selecting outside vertices which would increase the density of
$\tilde{S}$ when appended. This method is based on the density of $\tilde{S}$,
$\rho(\tilde{S})$ and the number of edges $\tilde{e}$, induced by the set of
vertices inside $\tilde{S}$ and the vertex in consideration of. As depicted by
figure 6 and figure 6, the number edges marked in thick lines should be
counted as $\tilde{e}$. It is obvious that if we add a new vertex that
introduces $\tilde{e}$ number of suitable edges such that
$\tilde{e}>\rho(\tilde{S})$, then $\rho(\tilde{S})$ is increased by an amount
of $(\textit{n}\tilde{e}-\textit{e})/(\textit{n}(\textit{n}+1))(>0)$, where n
is the number of vertices and e is the number of edges inside $\tilde{S}$
before the addition of the new vertex.

Figure 4: Second part of P-Bahmani with synchronization

Figure 5:

Figure 6:

Our approach towards designing the algorithm consists of two distinct
procedures. They are,

  1. 1.

Finding the densest k-core using an appropriate k-core decomposition
algorithm. Here we have to use two different decomposition algorithms for
serial and parallel versions.

  2. 2.

Selecting suitable vertices which can be added to $\tilde{S}$ and appending
them to $\tilde{S}$.

####  3.2.2 CBDS-P Algorithm

When designing the algorithm, we utilize the parallel k-core decomposition
algorithm, PKC, proposed by Kabir et al. [55].

CBDS-P takes a graph and outputs an approximate densest subgraph. The first
phase of the algorithm which is based on PKC algorithm, determines the core-
ness values of the vertices. Apart from calculating the core-ness values, the
first phase is modified to cater the self-edges of vertices and to calculate
the densities of detected cores.

The second phase of CBDS-P detects suitable vertices outside the densest core.
Threads find vertices which has a higher degree than the density of the
densest core and a lesser degree than the core-ness value of the densest core.
Then these selected vertices are processed in order to find vertices which
possess legitimate edge counts that satisfy the conditions presented in the
beginning of this subsection. We call these vertices, legitimate vertices.
Ultimately the legitimate vertices’ degrees are set to max_density_core and
the maximum density is calculated. Algorithm 2 returns the maximum density
(max_density), maximum density core value (max_density_core) and the input
graph. All the vertices that belong to the approximate densest subgraph can be
identified by checking their degree values, degree values are greater than or
equals to max_density_core. Note that some vertices may have a degree value
equal to the max_density_core even though they do not belong to the densest
core (i.e., here max_density_core value also serves as a labelling value).

First phase of CBDS-P has the time complexity equivalent to the time
complexity of PKC. It is $O(|V|K_{max}+|E|)$. Here, $K_{max}$ is the largest
core value. In second phase, we go through the complete graph in order to
detect the suitable vertices, adding a $O(|V|)$ complexity. It takes another
$O(|V^{{}^{\prime}}|)$ to go through $|V^{{}^{\prime}}|$ number of neighbours
introduced by the detected eligible vertices. Finally, to calculate the number
of edges between the legitimate vertices, it takes
$O(|V^{{}^{\prime\prime}}|^{2}-|V^{{}^{\prime\prime}}|)/2$.
$|V^{{}^{\prime\prime}}|$ is the number of legitimate vertices. Therefore, the
time complexity of CBDS-P is
$O(|V|(K_{max}+1)+|V^{{}^{\prime}}|+|V^{{}^{\prime\prime}}|^{2}+|E|)$. Based
on our experimental results, it is evident that $|V^{{}^{\prime}}|\ll|V|$ and
$|V^{{}^{\prime\prime}}|\ll|V|$.

To implement CBDS-P, we followed the same design concepts and technologies we
used to implement P-Bahmani.

Input : G(V, E)

43 Begin

44 $visited\leftarrow 0;max\\_density\leftarrow 0;max\\_density\leftarrow
0;max\\_density\\_core\leftarrow 0$;

45 $intermediate\\_edges\leftarrow 0;m_{v}\leftarrow 0;m_{e}\leftarrow 0$;

46 $legit\\_vector;$$\hfill\mathrm{\triangleright Global\hskip
2.84544ptvector\hskip 2.84544ptof\hskip 2.84544pttype\hskip
2.84544pt\textit{int}}$

47 parallel do

48 $k\leftarrow 0;start\leftarrow 0;end\leftarrow 0;deleted\\_local\leftarrow
0;aux\\_local\leftarrow 0;density\leftarrow 0$;

49 $buff;$$\hfill\mathrm{\triangleright Thread-local\hskip
2.84544ptarray\hskip 2.84544ptof\hskip 2.84544pttype\hskip
2.84544pt\textit{int}\hskip 2.84544ptwith\hskip 2.84544ptsize\hskip
2.84544pt|V|/N,\hskip 2.84544ptN=num\hskip 2.84544ptof\hskip
2.84544ptthreads}$

50 $eligible\\_vector;$$\hfill\mathrm{\triangleright Thread-local\hskip
2.84544ptvector\hskip 2.84544ptof\hskip 2.84544pttype\hskip
2.84544pt\textit{int}}$

51 while  _$visited <G\rightarrow|V|$_ do

52 parallel for  _$\nu\in G\rightarrow V$_ do

53 if  _$\nu.deg=k$_ then

54 $buff[end]\leftarrow\nu$;

55 $end++$;

56

57 end if

58

59 end for

60 while  _$start <end$_ do

61 $v\leftarrow buff[start];deleted\\_local\leftarrow
deleted\\_local+k;start++$;

62 for  _$u\in\nu.neighbours$_ do

63 if  _$\nu\neq u$_ then

64 if  _$u.deg >k$_ then

65 $a\leftarrow\textbf{atomicSub}(u.deg,1)$;

66 if  _$a=k+1$_ then

67 $buff[end]\leftarrow u;end++$;

68

69 end if

70 if  _$a\leq k$_ then

71 $\textbf{atomicAdd}(u.deg,1)$;

72

73 end if

74 if  _$a >k$_ then

75 $aux\\_local++$;

76

77 else

78 $deleted\\_local++;aux\\_local--$;

79 end if

80

81 end if

82

83 end if

84

85 end for

86

87 end while

88
$\textbf{atomicAdd}(visited,end);\textbf{atomicAdd}(deleted,deleted\\_local)$;

89 $\textbf{atomicAdd}(aux,aux\\_local)$;

90 barrier synchronization

91 single do

92
$density\leftarrow(G\rightarrow|E|-(deleted+aux)/2)/(G\rightarrow|V|-visited)$;

93 if  _$density >max\\_density$_ then

94 $max\\_density\leftarrow density;m_{e}\leftarrow
G\rightarrow|E|-((deleted+aux)/2)$;

95 $m_{v}\leftarrow G\rightarrow|V|-visited;max\\_density\\_core\leftarrow k$;

96

97 end if

98 end single

99 $start\leftarrow 0;end\leftarrow 0;deleted\\_local\leftarrow
0;aux\\_local\leftarrow 0;k++$;

100

101 end while

102

103

104

Algorithm 2 CBDS-P

49

50 $legits\leftarrow 0$;

51 parallel for  _$\nu\in G\rightarrow V$_ do

52 if  _$\nu.deg >max\\_density\textbf{ and }\nu.deg<max\\_density\\_core$_
then

53 $eligible\\_vector.add(\nu)$;

54

55 end if

56

57 end for

58 $start\leftarrow 0$;

59 while  _$start <eligible\\_vector.currentSize$_ do

60 $legits\leftarrow 0;\nu\leftarrow eligible\\_vector[start]$;

61 for  _$u\in\nu.neighbours$_ do

62 if  _$u.deg\geq max\\_density\\_core$_ then

63 $legits\leftarrow legits+1$

64 else if  _$\nu=u$_ then

65 $legits\leftarrow legits+0.5$

66 end if

67

68 end for

69 if  _$legits >max\\_density$_ then

70 critical do

71 $intermediate\\_edges\leftarrow intermediate\\_edges+legits$;

72 $legit\\_vector.add(\nu)$;

73

74 end critical

75 end if

76 $start++$;

77

78 end while

79 barrier synchronization

80 $i\leftarrow 0$;

81 parallel dynamic reduction(+: intermediate_edges)

82 for  _$i <legit\\_vector.currentSize-1$_ do

83 $\nu\leftarrow legit\\_vector[i];\nu.deg\leftarrow
max\\_density\\_core;j\leftarrow i+1$;

84 for  _$j <legit\\_vector.currentSize$_ do

85 $u\leftarrow legit\\_vector[j]$;

86 if  _$u\in\nu.neighbours\textbf{ and }\nu\neq u$_ then

87 $intermediate\\_edges\leftarrow intermediate\\_edges+1$;

88

89 end if

90 $j++$;

91

92 end for

93 $i++$;

94

95 end for

96 barrier synchronization

97 single do

98 $m_{e}\leftarrow m_{e}+intermediate\\_edges;m_{v}\leftarrow
m_{v}+legit\\_vector.currentSize$;

99 $max\\_density\leftarrow m_{e}/m_{v}$;

100 end single

101 end parallel

102 End

103

$\textbf{Return: }G,max\\_density\\_core,max\\_density$

##  4 Experimental Results and Evaluation

###  4.1 Datasets

All the datasets that we have used to evaluate the performance of our
algorithms are real-world datasets (natural graphs). They were obtained from
the Stanford Network Analysis Project (SNAP) and are publicly available for
anyone to use[]. Table 1 provides the details of the datasets. In order to
ensure the accuracy of our results obtained from our implementations, we used
datasets and results given in [9] as the reference. The table 2 provides the
table from the paper. Each of the algorithms are evaluated based on accuracy
and runtime.

Table 1: Dataset details Dataset | Label | $|V|$ | $|E|$
---|---|---|---
ca-GrQc | GrQc | 5,242 | 28,980
ca-HepTh | HepTh | 9,877 | 51,971
facebook-combined | Facebook | 4,039 | 88,234
ca-HepPh | HepPh | 12,008 | 237,010
musae-squirrel-edges | Squirrel | 5,201 | 216,933
as-skitter | Skitter | 1,696,415 | 11,095,298
soc-LiveJournal1 | LJ | 4,847,571 | 68,993,773
com-Orkut | Orkut | 3,072,441 | 117,185,083
Table 2: Reference data G(V, E) | $|\textit{V}|$ |  | $|\textit{E}|$ | $\rho^{*}\textit{(G)}$ |  | $\rho^{*}\textit{(G)}/\tilde{\rho}\textit{(G)}$ |
---|---|---|---|---|---|---|---
|  |  |  |  | $\epsilon=0.005$ | $\epsilon=0.05$ | $\epsilon=0.5$
as2000102 | 6,474 |  | 13,233 | 9.29 | 1.229 | 1.268 | 1.194
ca-AstroPh | 18,772 |  | 396,160 | 32.12 | 1.147 | 1.156 | 1.273
ca-CondMat | 23,133 |  | 186,936 | 13.47 | 1.072 | 1.072 | 1.429
ca-GrQc | 5,242 |  | 28,980 | 22.39 | 1.000 | 1.000 | 1.395
ca-HepPh | 12,008 |  | 237,010 | 119.004 | 1.000 | 1.017 | 1.151
ca-HepTh | 9,877 |  | 51,971 | 15.50 | 1.000 | 1.000 | 1.356
email-Enron | 36,692 |  | 367,662 | 37.34 | 1.058 | 1.072 | 1.063

###  4.2 Testing Environment

We used a m5.16xlarge Amazon EC2 which possess a 64 core Intel(R) Xeon(R)
Platinum 8175M CPU @ 2.50GHz and a 246GB RAM.

###  4.3 Density and Runtime Comparisons

Table 3 contains the approximate maximum density values calculated by both
P-Bahmani with $\epsilon=0$ and CBDS-P. Note that when $\epsilon=0$, P-Bahmani
algorithm produces the highest accuracy and it is equivalent to
2-approximation algorithm by Charikar et. al[]. Therefore, in the table 3, we
can clearly observe that CBDS-P calculates more accurate maximum density
values than a 2-approximation.

Figures 9, 9, 9, 12 for first 5 graphs and figures 15, 15, 15, 17 for last 2
graphs are the runtime variation of P-Bahmani against number of cores for
$\epsilon$ values of 0.5, 0.05, 0.005 and 0 respectively.

Runtime variation of graphs as-skitter, soc-LiveJournal1 and com-orkut against
the number of cores with respect to both P-Bahmani($\epsilon=0$) and CBDS-P
are shown in figures 12, 19 and 19. Further, figure 17 depicts the runtime of
CBDS-P against the number of cores with respect to both soc-LiveJournal1 and
com-orkut graphs.

Table 3: Density comparison Dataset | Exact Density | P-Bahmani $(\epsilon=0)$ | CBDS-P
---|---|---|---
ca-GrQc | 22.391 | 22.391 | 22.391
ca-HepTh | 15.5 | 15.5 | 15.5
facebook-combined | 77.347 | 69.9679 | 77.347
ca-HepPh | 119.004 | 119.004 | 119.004
musae-squirrel-edges | 152.543 | 131.779 | 152.542
as-skitter | 89.4696 | 71.7748 | 89.4696
soc-LiveJournal1 | — | 308.551 | 320.114
com-Orkut | — | 220.063 | 227.866

Figure 7:

Figure 8:

Figure 9:

Figure 10:

Figure 11:

Figure 12:

Figure 13:

Figure 14:

Figure 15:

Figure 16:

Figure 17:

Figure 18:

Figure 19:

##  5 Conclusion and Future Work

Even though the time complexity is higher than other 2-approximation
algorithms, our CBDS-P algorithm is capable of producing approximate densest
subgraphs that are more accurate than those 2-approximation algorithm within a
considerable time complexity. CBDS-P proposes a new approach towards
calculating approximate densest subgraphs based on k-cores.

P-Bahmani, the shared memory implementation of Bahmani et. al[] algorithm, is
the first of its kind for that algorithm. Both P-Bahmani and CBDS-P are highly
scalable in a shared memory, multi-core setting and it is evident by the
experimental results. Further, this paper clearly expatiates the design
considerations of our shared memory implementations.

Attempting to improve CBDS-P algorithm and working on shared memory
implementations for other densest subgraph discovery algorithms can be
identified as our future work.

## References

  * [1] Martinovic G., Krpic Z., and Rimac-Drlje S. “Parallelization Programming Techniques: Benefits and Drawbacks”. In: CLOUD COMPUTING 2010 (2010), pp. 7–13.
  * [2] Prema S., and Jehadeesan R. Analysis of parallelization techniques and tools. International Journal of Information and Computation Technology, 3(4):471–478, 2013.
  * [3] V. Lee, N. Ruan, R. Jin and C. Aggarwal. A survey of algorithms for dense subgraph discovery. Managing and Mining Graph Data, pages 303–336, 2010.
  * [4] Wu B., and Shen H. Exploiting efficient densest subgraph discovering methods for big data. IEEE Transactions on Big Data, 3(3):334–348, 2017.
  * [5] Rozenshtein P., Tatti N., and Gionis A. Finding dynamic dense subgraphs. ACM Transactions on Knowledge Discovery from Data (TKDD), 11(27), 2017.
  * [6] McGregor A., Tench D., Vorotnikova S., and Vu H.T. Densest subgraph in dynamic graph streams. MFCS 2015: Mathematical Foundations of Computer Science 2015, 9235:472–482, 2015.
  * [7] Y. Asahiro, R. Hassin, and K. Iwama. Complexity of finding dense subgraphs. Discrete Applied Mathematics, 121(1):15–26, 2002.
  * [8] Y. Asahiro, K. Iwama, H. Tamaki, and T. Tokuyama. Greedily finding a dense subgraph. Journal of Algorithms, 34(2):203–221, 2000.
  * [9] B. Bahmani, R. Kumar and S. Vassilvitskii, ”Densest subgraph in streaming and MapReduce”, Proceedings of the VLDB Endowment, vol. 5, no. 5, pp. 454-465, 2012. 10.14778/2140436.2140442
  * [10] V. Batagelj and M. Zaversnik. An o(m) algorithm for cores decomposition of networks. arXiv preprint cs/0310049, 2003.
  * [11] E. T. Charalampos. Mathematical and Algorithmic Analysis of Network and Biological Data. PhD thesis, Carnegie Mellon University, 2013. 2.1.
  * [12] M. Charikar. Greedy approximation algorithms for finding dense components in a graph. In APPROX, pages 84–95. Springer, 2000.
  * [13] J. Chen and Y. Saad. Dense subgraph extraction with application to community detection. TKDE, 24(7):1216–1230, 2012.
  * [14] Y. Chen, Y. Fang, R. Cheng, Y. Li, X. Chen, and J. Zhang. Exploring communities in large profiled graphs. IEEE Transactions on Knowledge and Data Engineering, 31(8):1624–1629, 2019.
  * [15] J. Cohen. Trusses: Cohesive subgraphs for social network analysis. National Security Agency Technical Report, 16, 2008.
  * [16] W. Cui, Y. Xiao, H. Wang, Y. Lu, and W. Wang. Online search of overlapping communities. In SIGMOD, pages 277–288. ACM, 2013.
  * [17] M. Danisch, T.-H. H. Chan, and M. Sozio. Large scale density-friendly graph decomposition via convex programming. In Proceedings of the 26th International Conference on World Wide Web, pages 233–242, 2017.
  * [18] A. Epasto, S. Lattanzi, and M. Sozio. Efficient densest subgraph computation in evolving graphs. In WWW, pages 300–310, 2015.
  * [19] Y. Fang, R. Cheng, Y. Chen, S. Luo, and J. Hu. Effective and efficient attributed community search. The VLDB Journal, 26(6):803–828, 2017.
  * [20] Y. Fang, R. Cheng, X. Li, S. Luo, and J. Hu. Effective community search over large spatial graphs. PVLDB, 10(6):709–720, 2017.
  * [21] Y. Fang, R. Cheng, S. Luo, and J. Hu. Effective community search for large attributed graphs. PVLDB, 9(12):1233–1244, 2016.
  * [22] Y. Fang, X. Huang, L. Qin, Y. Zhang, W. Zhang, R. Cheng, and X. Lin. A survey of community search over big graphs. The VLDB Journal, 2019.
  * [23] Y. Fang, Z. Wang, R. Cheng, X. Li, S. Luo, J. Hu, and X. Chen. On spatial-aware community search. IEEE Transactions on Knowledge and Data Engineering, 31(4):783–798, 2018.
  * [24] Y. Fang, Z. Wang, R. Cheng, H. Wang, and J. Hu. Effective and efficient community search over large directed graphs. IEEE Transactions on Knowledge and Data Engineering (TKDE), 2018.
  * [25] A. Gionis and C. E. Tsourakakis. Dense subgraph discovery: Kdd 2015 tutorial. In SIGKDD, pages 2313–2314, NY, USA, 2015. ACM.
  * [26] A. V. Goldberg. Finding a maximum density subgraph. UC Berkeley, 1984.
  * [27] J. Hu, R. Cheng, K. C.-C. Chang, A. Sankar, Y. Fang, and B. Y. Lam. Discovering maximal motif cliques in large heterogeneous information networks. In IEEE International Conference on Data Engineering (ICDE), pages 746–757. IEEE, 2019.
  * [28] J. Hu, X. Wu, R. Cheng, S. Luo, and Y. Fang. Querying minimal steiner maximum-connected subgraphs in large graphs. In International on Conference on Information and Knowledge Management (CIKM), pages 1241–1250. ACM, 2016.
  * [29] J. Hu, X. Wu, R. Cheng, S. Luo, and Y. Fang. On minimal steiner maximum-connected subgraph queries. IEEE Transactions on Knowledge and Data Engineering, 29(11):2455–2469, 2017.
  * [30] X. Huang, H. Cheng, L. Qin, W. Tian, and J. X. Yu. Querying k-truss community in large and dynamic graphs. In SIGMOD, pages 1311–1322. ACM, 2014.
  * [31] X. Huang and L. V. Lakshmanan. Attribute-driven community search. PVLDB, 10(9):949–960, 2017.
  * [32] X. Huang, W. Lu, and L. V. Lakshmanan. Truss decomposition of probabilistic graphs: Semantics and algorithms. In Proceedings of the 2016 International Conference on Management of Data, pages 77–90. ACM, 2016.
  * [33] M. Jha, C. Seshadhri, and A. Pinar. Path sampling: A fast and provable method for estimating 4-vertex subgraph counts. In WWW, pages 495–505, 2015.
  * [34] R. Kannan and V. Vinay. Analyzing the structure of large graphs. Rheinische Friedrich-Wilhelms-Universität Bonn, 1999.
  * [35] S. Khuller and B. Saha. On finding dense subgraphs. Automata, Languages and Programming, pages 597–608, 2009.
  * [36] L. Lü, T. Zhou, Q.-M. Zhang, and H. E. Stanley. The h-index of a network node and its relation to degree and coreness. Nature communications, 7:10168, 2016.
  * [37] M. Mitzenmacher, J. Pachocki, R. Peng, C. Tsourakakis, and S. C. Xu. Scalable large near-clique detection in large-scale networks via sampling. In International Conference on Knowledge Discovery and Data Mining (SIGKDD), pages 815–824. ACM, 2015.
  * [38] Y. Peng, Y. Zhang, W. Zhang, X. Lin, and L. Qin. Efficient probabilistic k-core computation on uncertain graphs. In International Conference on Data Engineering (ICDE), pages 1192–1203. IEEE, 2018.
  * [39] L. Qin, R.-H. Li, L. Chang, and C. Zhang. Locally densest subgraph discovery. In KDD, pages 965–974. ACM, 2015.
  * [40] S.Sahu,A.Mhedhbi,S.Salihoglu,J.Lin,and M.T.Özsu. The ubiquity of large graphs and surprising challenges of graph processing. PVLDB, 11(4):420–431, 2017.
  * [41] R. Samusevich, M. Danisch, and M. Sozio. Local triangle-densest subgraphs. In International Conference on Advances in Social Networks Analysis and Mining (ASONAM), pages 33–40. IEEE, 2016.
  * [42] A. E. Sariyüce and A. Pinar. Fast hierarchy construction for dense subgraphs. PVLDB, 10(3):97–108, 2016.
  * [43] A. E. Sariyüce, C. Seshadhri, and A. Pinar. Local algorithms for hierarchical dense subgraph discovery. PVLDB, 12(1):43–56, 2018.
  * [44] A. E. Sariyüce, C. Seshadhri, A. Pinar, and U.V.Çatalyürek. Finding the hierarchy of dense subgraphs using nucleus decompositions. In WWW, pages 927–937, 2015.
  * [45] A.E.Sariyüce, C.Seshadhri, A.Pinar,and U.V.Çatalyürek. Nucleus decompositions for identifying hierarchy of dense subgraphs. ACM Transactions on the Web (TWEB), 11(3):16, 2017.
  * [46] S. B. Seidman and B. L. Foster. A graph-theoretic generalization of the clique concept. Journal of Mathematical sociology, 6(1):139–154, 1978.
  * [47] N. Tatti. Density-Friendly Graph Decomposition. In: ACM Transactions on Knowledge Discovery from Data 13.5 (2019), pp. 1–29. 10.1145/3344210.
  * [48] N. Tatti and A. Gionis. Density-friendly graph decomposition. In Proceedings of the 24th International Conference on World Wide Web, pages 1089–1099, 2015.
  * [49] C. Tsourakakis. The k-clique densest subgraph problem. In WWW, pages 1122–1132, 2015.
  * [50] C. Tsourakakis, F. Bonchi, A. Gionis, F. Gullo, and M. Tsiarli. Denser than the densest subgraph: extracting optimal quasi-cliques with quality guarantees. In KDD, pages 104–112. ACM, 2013.
  * [51] K. Wang, X. Cao, X. Lin, W. Zhang, and L. Qin. Efficient computing of radius-bounded k-cores. In IEEE International Conference on Data Engineering (ICDE), pages 233–244. IEEE, 2018.
  * [52] K. Wang, X. Lin, L. Qin, W. Zhang, and Y. Zhang. Vertex priority based butterfly counting for large-scale bipartite networks. PVLDB, 12(10):1139–1152, 2019.
  * [53] Y. Wu, R. Jin, J. Li, and X. Zhang. Robust local community detection: on free rider effect and its elimination. PVLDB, 8(7):798–809, 2015.
  * [54] J. Liu, S. Salehian, and Y. Yan. Comparison of Threading Programming Models. In: 2017 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW). IEEE, 2017. 10.1109/ipdpsw.2017.141.
  * [55] H. Kabir, and K. Madduri “Parallel k-Core Decomposition on Multicore Platforms”. In: 2017 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW). IEEE, 2017. 10.1109/ipdpsw.2017.151.
  * [56] Y. Fang, K. Yu, R. Cheng, L. V. S. Lakshmanan, and X. Lin. 2019. Efficient algorithms for densest subgraph discovery. Proc. VLDB Endow. 12, 11 (July 2019), 1719–1732. 10.14778/3342263.3342645


††thanks: L. R. and V. I. contributed equally to this work

# Effects of the rf current and bias field direction on the transition from
linear to non-linear gyrotropic dynamics in magnetic vortex structures

Lakshmi Ramasubramanian  Institute of Ion Beam Physics and Materials Research,
Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden, Germany  Institute of
Physics, Chemnitz University of Technology, 09107 Chemnitz, Germany  Vadym
Iurchuk11footnotemark: 1<EMAIL_ADDRESS>Institute of Ion Beam Physics and
Materials Research, Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden,
Germany  Olav Hellwig  Institute of Ion Beam Physics and Materials Research,
Helmholtz-Zentrum Dresden-Rossendorf, 01328 Dresden, Germany  Institute of
Physics, Chemnitz University of Technology, 09107 Chemnitz, Germany  Alina
Maria Deac  Dresden High Magnetic Field Laboratory (HLD-EMFL), Helmholtz-
Zentrum Dresden-Rossendorf, 01328 Dresden, Germany

###### Abstract

We present a frequency-domain study of the dynamic behavior of a magnetic
vortex core within a single Permalloy disk by means of electrical detection
and micromagnetic simulations. When exciting the vortex core dynamics in a
non-linear regime, the lineshape of the rectified dc signal reveals a
resonance peak splitting which depends on the excitation amplitude. Using
micromagnetic simulations, we show that at high excitation power the peak
splitting originates from the nanosecond time scale quasi-periodic switching
of the vortex core polarity. Using lock-in detection, the rectified voltage is
integrated over a ms time scale, so that the net signal detected between the
two resonant peaks for a given range of parameters cancels out. The results
are in agreement with the reported effects of the in-plane static field
magnitude on the gyration dynamics, and complement them by detailed analysis
of the effects of the rf current amplitude and the azimuthal angle of the in-
plane bias magnetic field. Systematic characterization shows that a transition
from linear to nonlinear dynamical regime can be controlled by rf current as
well as by varying the magnitude and the direction of the bias magnetic field.

##  1 Introduction

Magnetic vortices, which are spontaneously formed in magnetic disks of
specific aspect ratios, are fundamental, topologically protected magnetic
structures with an in-plane circulating magnetization and an out-of-plane
magnetized vortex core (VC), which at remanence is localized in the center of
the disk [1, 2]. Since the first experimental observation of vortex formation
in Py disks [3], the physics of magnetic vortices has evolved into a large and
intense field of research that holds promises for future technological
applications. Indeed, magnetic disks in the vortex state have the potential
for application in biomedicine as magnetic microparticles that help to destroy
cancer cells [4, 5] or act as drug delivery agents [6]. Asymmetric systems of
exchange-coupled vortices can also serve as sources for the generation of
short wave-length dipole-exchange spin-waves of a directional nature, with
potential applications in magnonic computing [7]. Recently, rich non-linear
spin-wave dynamics in a magnetic vortex was revealed including excitation of
whispering gallery magnons via non-linear 3-magnon scattering processes [8, 9]
and VC gyration-assisted excitation of spin-wave frequency combs in the GHz
range [10] giving prospect for using magnetic vortices as magnonic reservoirs
for unconventional data processing schemes.

In this paper, we present an extensive study of the non-linear gyration
dynamics in a magnetic vortex investigated using homodyne electrical detection
techniques. A similar experimental study of the non-linear gyration dynamics
showed the in-plane bias field-induced resonance splitting [11]. Based on a
generic mechanical model of a nonlinear oscillator, the observed splitting was
attributed to the dynamical bifurcation in the non-linear regime. Our work
demonstrates that the origin of the resonant peak splitting is the nanosecond
time-scale quasi-periodic VC switching at high excitation amplitudes, which
can be achieved by varying not only the bias field magnitude, but also its
azimutal angle as well. We show that the transition from the linear to the
non-linear dynamical regime in single Py disk can be efficiently manipulated
at constant rf excitation power by varying the external control parameters,
i.e. the bias magnetic field magnitude and direction.

Figure 1: (a) Schematics of the measurement setup with magnetic field
modulation of the magnetoresistance enabling rectified dc voltage measurements
via lock-in technique. To generate a modulation field (collinear to the static
bias magnetic field) at the lock-in reference frequency (here 223 Hz), copper
wire coils (red) were attached to the electromagnet poles (blue). (b) Typical
magnetoresistance loops measured from a disk (here radius $R=$ 4.4 µm,
thickness $t=$ 35 nm) for a dc current $I_{dc}=$ 5 mA). Blue solid (red
dashed) lines are the backward (forward) field sweep measurements data from
negative (positive) saturation as indicated by the corresponding arrow. The
insets depict the schematic direction of the VC (orange dot) displacement
depending on the orientation of the external field $H_{b}$.

##  2 Experimental methods

Magnetic microdisks were fabricated on Si/SiO2 substrates using electron beam
lithography, followed by the deposition of a 35 nm Permalloy (Py = Ni80Fe20)
film by e-beam evaporation and conventional lift-off. As the next step, Cr/Au
contact pads were deposited to provide electrical access to each microdisk. We
used a standard magnetotransport setup with rf capability (see Fig. 1(a)) for
electrical detection of magnetization dynamics in single magnetic vortices at
room temperature. The detection technique exploits the anisotropic
magnetoresistance (AMR) effect, i.e. the resistance change with the relative
angle between the direction of the electrical current and the magnetic moment
of a magnetic structure [12]. An rf current injected through a bias-T into the
microdisk device induces the gyrotropic motions of the VC at resonance, via
the joint action of the spin-transfer torque and rf Oersted field, and thereby
leads to an oscillating AMR provided the VC is displaced from the center of
the disk. The time-averaged product of the rf current and the dynamical
resistance – which results in a rectified dc voltage – is measured by a
conventional homodyne detection scheme using a lock-in amplifier. To improve
the signal-to-noise ratio, magnetic field modulation at the lock-in reference
frequency was used [13]. The details on the sample fabrication process, as
well as the experimental procedure for electrical detection, have been
reported elsewhere [14]. One has to note, that we use Py disks with relatively
large radii (3.2 µm and 4.4 µm). These relatively large sizes were initially
chosen for Cr+ ion implantation in order to build double output frequency
devices. The results will not be covered here, as they have been reported
elsewhere [14].

Figure 2: Spectral color maps of the excitation power-dependent rectified dc
signal measured of a Py disk ($R=$ 3.2 $\mu$m, $t=$ 35 nm) taken at a static
field $\mu_{0}H_{b}$ = 2 mT for (a) $H_{b}\perp I$ and (b) $H_{b}\parallel I$.
The insets show the schematic position of the VC for the given field
orientation. Threshold current $j_{th}$ (marked by a vertical dashed line in
(a) and (b)) corresponds to the current density for which the signal peak can
be clearly detected. Dashed sloped lines in (a) are guides to the eye showing
schematically the peak/dip frequency shifts with increasing $j_{rf}$.
$j_{split}$ corresponds to the current density where the frequency splitting
is detected. (c,d) Single $V_{dc}$ spectra measured at high/low rf current
density ($j_{rf}=$ 0.4 $\cdot 10^{11}$ A/m2 and $j_{rf}=$ 0.2 $\cdot 10^{11}$
A/m2 respectively) for $H\perp I$. e) Zero-signal gap for $H\perp I$ as a
function of $j_{rf}$. (f,g) Single $V_{dc}$ spectra measured at the same
values of $j_{rf}$ as in (c,d) for $H\parallel I$.

##  3 Results and discussion

###  3.1 Detection of the vortex core nucleation using static
magnetoresistance measurements

We first confirmed the presence of a magnetic vortex in the fabricated disks
by measuring the static magnetoresistance (MR). It has to be noted that in
order to reveal a measurable MR signal, a relatively large dc current (about
50 mA) had to be applied through the as-prepared devices in order to break the
thin Py oxide layer formed during the fabrication process, which isolates the
magnetic disk from the contact pads. As the oxide does not break homogeneously
under the contact pad, this leads to a presumably non-uniform electrical
current distribution within the disk, whose effects on the dynamics are
discussed further in the text. Fig. 1(b) shows the MR of the disk at a dc
current I${}_{dc}=$ 5 mA measured for in-plane bias magnetic field
$\mu_{0}H_{b}$ applied parallel (top curve) and perpendicular (bottom curve)
to the direction of the probing current $I$ .

Starting from positive in-plane saturation field, i.e. with the disk in a
single domain state, the external field is reduced to nucleate a vortex in the
disk. This is observed in our measurements as a resistance decrease (increase)
for $H_{b}\parallel I$ ($H_{b}\perp I$). Further sweeping the field towards
negative values displaces the VC within the disk (in the direction
perpendicular to the applied field) until its expulsion, which is accompanied
by a recovery of the resistance value of the saturated state. The same
measurement is repeated in the reverse direction (from negative to positive
saturation field) to reproduce the vortex nucleation/expulsion behaviour.

###  3.2 Vortex core dynamics at large excitation power: resonant bistability

####  3.2.1 Vortex resonance splitting in strongly nonlinear regime

The gyrotropic motion of the VC is probed by frequency-resolved homodyne
detection using the spin rectification effect. The rf frequency sweeps were
performed varying one of the three external parameters, i.e. the amplitude of
the excitation rf current density $j_{rf}$, the bias magnetic field magnitude
$H_{b}$ or the in-plane azimuthal angle $\phi_{H}$ of the external bias field
while keeping the other two constant. Here, the value of the current density
$j_{rf}$ injected into the device can only be estimated taking into account
the impedance mismatch between the device and the 50$\Omega$ circuit. In this
section, we focus on the effect of the $j_{rf}$ amplitude on the VC dynamics.

To understand the dynamic behaviour of the VC at high excitation power, we
first performed frequency sweeps at a fixed in-plane bias field $\mu_{0}H_{b}$
= 2 mT (for $H_{b}\perp I$) while increasing rf power (see Fig. 2(a)) in order
to observe detectable gyration dynamics signal. A first clear measurable
response was observed at a threshold current $j_{th}=$ 0.08 $\cdot 10^{11}$
A/m2. With increasing $j_{rf}$, the rectified signal also increases as can be
seen from the intensity map of Fig. 2(a), since the VC now gyrates with a
larger radius, leading to a larger MR change. Fig. 2(d) shows a clear
antisymmetric lorentzian-like signal [15] acquired at $j_{rf}=$ 0.2 $\cdot
10^{11}$ A/m2. Another striking feature appears when reaching $j_{split}=$
0.23 $\cdot 10^{11}$ A/m2, where the dc voltage signal evolves from a single
antisymmetric peak response to a ’peak-dip’ lineshape separated by a plateau
(gap) between them (Fig 2(a)). This zero-signal gap expands with increasing
$j_{rf}$, making the splitting in the resonance lineshape more pronounced (see
Fig. 2(a,e)). The splitting behaviour is considered as an indication of the
vortex entering a non-linear bistability region. This effectively means that
in the strongly non-linear regime, the system becomes resonantly stable at two
excitation frequencies.

The clear double-peak response is not reproduced for $H_{b}\parallel I$ (see
Fig. 2(b,f,g)). Indeed, for this field orientation, the threshold current for
the dc voltage detection was 0.18 $\cdot 10^{11}$ A/m 2, which is 2 times
higher as compared to the $H_{b}\perp I$ case. Moreover, for $H_{b}\parallel
I$, the rectified signal amplitudes are several times lower and the peak
splitting at high $j_{rf}$ is not clearly developed (see Fig. 2(f)).

Figure 3: (a) SEM image of the shorted device after $\sim$60 mA dc current was
passed through the disk. The brighter SEM contrast reveals the areas where the
Permalloy oxide was broken and the current density within the disk was the
highest. (b) Simulated distribution of the current density in a 35 nm thick Py
disk ($R$ = 2.1 $\mu$m) calculated using COMSOL Multiphysics®. The scale bar
shows the current density in A/m2. In this simulation, a 5 mA dc current is
passed through the gold pads ohmically contacted to the Py disk. The
experimentally observed breaking of the Permalloy oxide layer is mimicked by
considering the current injected into the disk through the point-like regions
at the edges of the gold pads (here taken randomly as circular areas with 20
nm radii). The result is the non-homogeneous distribution of the injected
current where higher current density is observed closer to the contact pads.

This spacial anisotropy in the VC dynamics is attributed to the inhomogeneous
current distribution within the magnetic disk. Fig. 3(a) shows an SEM image of
a successfully contacted device, taken after a large dc current of 60 mA was
passed through leading to the on-purpose device shortcut. The brighter
contrast region within the disk indicates the current flow path, confirming
the point-like oxide breaking and the corresponding inhomogeneous current
distribution leading to the different dynamics response depending on the bias
field orientation. Fig. 3(a) shows the simulated current distribution in the
disk assuming that the thin oxide layer on top of the Py film is broken so
that it forms a point-like electrical contact. This assumption agrees with the
experimentally observed drop of the device resistance from $\sim$130 $\Omega$
to $\sim$12 $\Omega$ at 50 mA dc current applied to the as-prepared disk in
order to break the oxide layer and establish a direct conductive contact
between the disk and the pads. In general, the local current density closer to
the contact pads is higher compared to the disk center. Therefore, when the VC
is displaced closer to the contact pads (i.e. when $H_{b}\perp I$), the local
rf Oe field (generated by the injected rf current) driving its dynamics is
expected to be higher as compared to the $H_{b}\parallel I$ case (when the VC
is located far from the contract pads). This argument also explains the
differences in the threshold current values and the peak splitting behavior in
the experimental data of Fig. 2.

A bistable behaviour similar to the one observed here in a single disk has
previously been reported on an array of $\sim$2000 magnetic vortex structures
probed by microwave reflection measurements at high excitation amplitudes
[16]. Electrical detection of the double resonance response of the VC dynamics
has also been reported, where such drift from linear gyrotropic motion is
generally linked with non-linear phenomena that can be understood by
considering the system as a non-linear oscillator for which two stable
dynamical states co-exist (one with large and one with small oscillation
amplitude) [11]. However, our analysis shows that in single magnetic disks,
the double-resonance behavior and the appearance of a zero-signal shoulder
between the two resonances originates from quasi-periodic switchings of the VC
under extremely strong rf current excitation in agreement with the non-linear
analytical model of the vortex dynamics [17]. This mechanism is revealed using
micromagnetic simulations and is explained in detail in the next section.

####  3.2.2 Dynamical micromagnetic simulations

We use the GPU-accelerated MuMax3 software package [18] for simulating the
magnetization dynamics in the Py magnetic disks. We study the trajectories of
the vortex gyrotropic mode as a function of the frequency of the rf magnetic
field excitation. Similar to [19], the VC is excited by an in-plane sinusoidal
rf magnetic field $h_{rf}=h\sin{(2\pi ft)}$, where $f$ is the excitation
frequency, and $h=(h^{x}_{rf},0,0)$. The magnetization dynamics in the
magnetic disk is simulated for different rf field amplitudes $h^{x}_{rf}$ over
100 gyration periods $T_{0}=\frac{1}{f}$. The VC trajectory over the last 10
periods is recorded for each value of the excitation frequency $f$ in the
chosen range. We use the following Py material parameters for the
magnetization dynamics simulations: saturation magnetization $M_{s}$ = 800
kA/m, exchange constant $A_{ex}$ = 13 pJ/m3, and damping parameter $\alpha$ =
0.008 [19]. To reduce the simulation time, we consider a much smaller Py
magnetic disk of 250 nm radius and 30 nm thickness discretized into 128×128×1
cells. Although the corresponding VC gyration resonance frequency for such
disk ($\sim$ 470 MHz) is approximately one order of magnitude higher compared
to the experimentally studied Py disk ($\sim$ 49 MHz), the evolution of the VC
resonance behavior as a function of the rf field amplitude is qualitatively
similar.

Figure 4: (a–c) Simulated trajectories of the VC for a Py disk ($R$ = 250 nm;
$t$ = 30 nm) at $f$ = 470 MHz and $h^{x}_{rf}$ = 0.5, 1 and 2 mT. (d–f)
Corresponding time traces of the $m_{x}$ (top) and $m_{y}$ (bottom) components
of the normalized magnetization recorded for the last 10 periods ($\sim$21.3
ns) of the simulation time. (g–i) Trajectory opening vs. excitation frequency
for $f$ = 470 MHz and $h^{x}_{rf}$ = 0.5, 1 and 2 mT. (j–l) Experimental dc
voltage spectra measured at $\mu_{0}H_{b}$ = 2 mT for $H_{b}\perp I$ taken
from the spectral map of Fig. 3(a) at given $j_{rf}$ values. Shaded region in
(i) and (l) refers to the VC switching range.

Fig. 4(a,b) shows the simulated trajectories of the VC at $f$ = 470 MHz and
$h^{x}_{rf}$ = 0.5 mT (a) and 1 mT (b). For 0.5 and 1 mT amplitudes the steady
state gyration sets after $\sim$20 gyration periods with quasi-circular VC
trajectories. The time traces of the $m_{x}$ and $m_{y}$ magnetization
components taken after the steady state gyration sets in (see Fig. 4(d,e))
show harmonic oscillations over time. The maximum VC displacement (trajectory
opening) in the $x$ direction is extracted from the simulations and is plotted
versus the frequency in Fig. 4(g,h) for $h_{x}$ = 0.5 and 1 mT respectively.
These VC displacement traces define the dynamical MR variation whose
derivative is proportional to the rectified dc voltage signal detected during
experiments at a given frequency (see Fig. 4(j,k)). Increasing the rf field
amplitude from 0.5 mT to 1 mT, the resonance peak broadens and the trajectory
opening (rectified signal amplitude) increases. However, for $h^{x}_{rf}$ = 2
mT, no steady state trajectory was observed, but instead we observe quasi-
periodic switching of the VC polarity. For example, as seen in Fig. 4(c), the
VC switches 2 times over the 10 gyration periods ($\sim$ 21.3 ns), from VC
polarity ”up” to ”down” and back to ”up”, which is accompanied by a change in
the gyration sense from counterclockwise to clockwise and back to
counterclockwise. It can be also seen from the corresponding time traces of
the $m_{x}$ and $m_{y}$ components (see Fig. 4(f)) that the phase of the
oscillations shifts twice over the 10 gyration periods indicating the core
switching events. Fig. 4(i) shows the frequency range for $h^{x}_{rf}$ = 2 mT
in the vicinity of the gyration resonance where the VC undergoes multiple
quasi-periodic polarity switchings. We find that these switchings occur via a
known mechanism of vortex-antivortex pair formation [20] and the subsequent
annihilation of the original vortex with the created antivortex, when the
vortex core reaches its critical velocity ($\sim$333 m/s for Py [21]). In our
simulations, for $h^{x}_{rf}$ = 2 mT, (Fig. 4(c,f)) we observe quasi-periodic
polarity switching events approximately every 3 gyration periods, in agreement
with the previously reported results for Py disks [22].

Core switching is accompanied by a change in the sense of gyration. Since the
electrical detection method is sensitive to the VC polarity [23, 15, 24], when
the vortex core flips millions of times over one lock-in time constant (250
ms), the dynamical resistance response will be counter balanced for each
subsequent polarity switching event. This leads to a zero net rectification
voltage for the rf frequencies when switchings take place and makes the
antisymmetric resonance lineshape to appear like a peak followed by a dip with
a pronounced flat plateau as seen in the corresponding experimental spectrum
(in Fig. 4(l)). By increasing the excitation amplitude, the observed plateau
widens significantly as now even off-resonant excitation is sufficiently
strong to induce the VC polarity switching (see Fig. 2(a)).

One has to comment, that in the frequency sweep experiments shown in this
work, the equilibrium position of the vortex core is displaced with a static
external field and an rf current is applied to the disk to excite the VC
gyration. However, in the simulations, the vortex core is excited by an rf
field and the gyration trajectories are recorded in the absence of any
external bias field. Nevertheless, the experiments and simulations show good
agreement as the vortex core is excited resonantly in both cases.

###  3.3 Controlling the transition from linear to nonlinear dynamics at fixed
rf current

As seen from the above described measurements and simulations, the transition
between linear and nonlinear dynamical regimes can be controlled by moving the
VC between the regions in the disk with lower/higher rf current density (i.e.
further/closer to the electrical contact pads) while keeping the input power
constant. This can be realized by changing the magnitude of $H_{b}$ and/or the
azimuthal angle $\phi_{H}$ of the external bias field.

####  3.3.1 $\phi_{H}$ dependence

Fig. 5 shows the spectral color map of the rectified dc voltage measured at
fixed rf current $j_{rf}$ = 0.2 $\cdot 10^{11}$ A/m2 and bias field
$\mu_{0}H_{b}$ = 2 mT for different azimuthal angles $\phi_{H}$ ranging from
0∘ ($H_{b}\parallel I$) to 90∘ ($H_{b}\perp I$). For $\phi_{H}$ = 0∘, a
minuscule rectified signal is observed (note that this combination of
$j_{rf}$, $H_{b}$ and $\phi_{H}$ corresponds to the spectrum shown in Fig.
2(f)), as the VC is located between the contact pads (see bottom left inset
for the schematic VC position) and the excitation power is slightly above the
threshold rf current $j_{th}$ = 0.18 $\cdot 10^{11}$ A/m2. When $\phi_{H}$ is
increased to $\sim$15∘, a clear single peak is detected at about 52 MHz. It is
attributed to the VC being displaced closer to the contact pads, where the VC
is subjected to a stronger rf Oersted field excitation due to the
inhomogeneous rf current distribution within the disk (see Fig. 3). Note that
here we refer to the Oe field generated by the current flowing along the z
direction, i.e. between the contact pads and the disk at the location of the
point-like oxide breakthrough. For $\phi_{H}\geq$ 50∘, a pronounced
antisymmetric signal is observed, whose amplitude increases with increasing
the angle $\phi_{H}$, reaching its maximum at $\phi_{H}$ = 90∘. This behavior
is consistent with the data of Fig. 2(a) where the peak amplitude increases
with increasing the rf current at $H_{b}\perp I$ (i.e. at fixed angle
$\phi_{H}$ = 90∘), but here the rf current remains fixed, and the effective
excitation is increased by increasing the azimuthal angle $\phi_{H}$ (i.e by
moving the VC closer to the contact leads, see bottom middle and right
insets).

Figure 5: Spectral color map of the VC gyrotropic frequency vs. azimuthal
angle $\phi_{H}$ measured for $j_{rf}$ = 0.2 $\cdot 10^{11}$ A/m2 and
$\mu_{0}H_{b}$ = 2 mT. Dashed lines are guides to the eye. Insets at the
bottom depict the schematical positions of the VC for $\phi_{H}$ = 0∘
(corresponding to $H\parallel I$), 45∘ and 90∘ (corresponding to $H\perp I$).

One has to note that the resonance frequency slightly decreases with
increasing angle, indicating a possible transition towards a non-linear
gyration regime. Indeed, for angles close to 90∘, a small gap is observed
between resonant peak and dip indicating a peak splitting behavior being a
typical feature of the non-linear VC dynamics as was in detail explained in
the previous section. The presented data shows that the amplitude of the VC
gyration (and even the transition form linear to non-linear dynamical regimes)
can be controlled at constant excitation rf power and fixed bias magnetic
field by changing the azimuthal field angle, i.e. by moving the VC between the
regions of different effective excitation strength (here closer to or further
away from the contact pads).

####  3.3.2 $H_{b}$ dependence

Up to now we discussed angular dependent and rf current dependent measurements
carried out at a fixed bias magnetic field of 2 mT. We have shown that due to
the inhomogeneous current distribution within the disk, the VC dynamical
response is defined not only by the excitation power but by the spatial
position of the VC with respect to the contact pads as well. It is expected
that one can drive the VC in a non-linear regime at fixed rf current by
shifting the VC towards the contact pads, thus increasing the effective Oe
field excitation. To verify this, we have performed dynamical measurements at
$H\perp I$ with $j_{rf}$ = 0.68 $\cdot 10^{11}$ A/m 2 for different values of
bias field $H_{b}$. Here, we keep the azimuthal angle $\phi_{H}$ = 90∘ fixed
(see the bottom right inset of Fig. 5).

Figure 6: Spectral color map of the rectification signal $V_{dc}$ as a
function of the VC gyrotropic frequency vs. external bias field $H_{b}$
measured for rf current density $j_{rf}$ = 0.68 $\cdot 10^{11}$ A/m2 and
$H_{b}\perp I$ (corresponding to $\phi_{H}$ = 90∘). Dashed lines are guides to
the eye showing the increase of the ”peak-dip” separation with increasing bias
field.

Fig. 6 shows the map of the rectified dc voltage frequency spectra acquired as
a function of the bias field $H_{b}$ applied perpendicular to the rf current
at $j_{rf}$ = 0.68 $\cdot 10^{11}$ A/m 2. At given current density and field
azimuthal angle, a linear dynamical response is expected at low bias field
(see Fig. 2(a)), which is seen as a single asymmetric peak for $H_{b}\leq$ 1
mT. When $\mu_{0}H_{b}$ is increased to $\sim$1.5 mT, an onset of the peak
splitting is observed, which increases with increasing the bias field, leading
to a pronounced ”peak-dip” signal shape similar to the one shown in Fig.
2(a,c) and Fig. 5 for large $j_{rf}$ and $\phi_{H}$ respectively. This
behavior indicates that increasing bias field drives the system to the non-
linear dynamical regime, as the VC is displaced to the region with larger rf
current density, i.e. with higher effective excitation power. At
$\mu_{0}H_{b}\sim$ 8 mT, the signal drops abruptly, corresponding to the VC
expulsion from the disk. When the field is swept to negative values, the
behavior is reproduced almost symmetrically. The minor asymmetry with respect
to the bias field may be attributed to a slightly off-centered position of the
VC at zero bias field due to disk imperfections and/or residual magnetic
fields between the electromagnet poles. For $H\parallel I$ (i.e. when the VC
is located in the region with lower current density), no peak splitting was
observed (the VC gyration remains within the linear regime) due to
insufficient effective excitation power for the VC dynamics.

Similar to the $\phi_{H}$ dependence, this measurement shows that by moving
the core between the regions of different effective excitation strength, one
can control the dynamical regime of the VC gyration and switching from linear
to non-linear dynamics at fixed rf excitation and bias field direction by
varying only the magnitude of the bias magnetic field. This fixed magnetic
field can also be introduced for instance by adding an adjacent magnetic layer
coupled magnetostatically to the disk in a vortex state.

##  4 Conclusions

We have studied the VC nonlinear dynamics in a single magnetic microdisk as a
function of external control parameters, i.e the excitation rf current, the
static bias magnetic field magnitude and the in-plane field orientation. Via
magnetotransport measurements, we have found that the effective excitation
power depends on the VC position within the disk and is essentially higher for
magnetic fields aligned perpendicular to the rf current flow. This anisotropy
of the dynamical response is attributed to the inhomogeneous rf current
distribution within the disk, leading to higher rf current-induced Oersted
field for the case when the VC is located close to the contact pads. This
effect enables tuning the transition from linear to non-linear dynamical
regimes by moving the VC between the regions with effectively lower/higher
excitation amplitude (by changing the bias field magnitude and/or in-plane
angle) while keeping the nominal rf current flowing through the disk fixed.
Furthermore, with higher resolution lithographic design of the contact pads,
it appears to be possible to control the exact location of the oxide
breakthrough point in order to optimize the parametric tunability of the
transition from linear to non-linear vortex dynamics.

The described complex dynamics may be of use in multifrequency rf detectors
based on magnetic vortices with two different working regimes (corresponding
to the linear and non-linear dynamical regimes of the magnetic vortex).
Depending on the selected dynamical regime, the magnetic vortex can be
resonantly excited by either a single input frequency (linear regime) or a
combination of two input frequencies (non-linear regime). In the non-linear
case, the sensitivity to the difference between the input frequencies is
defined by the plateau between the corresponding resonances and can be tuned
by varying the external control parameters, i.e. excitation power and/or bias
magnetic field magnitude/direction.

###### Acknowledgements.

This work is supported by the European Union’s Horizon 2020 research and
innovation programme under grant agreement No. 737038 (TRANSPIRE), the
Helmholtz Young Investigator Initiative Grant No. VH-N6-1048, and the HLD at
HZDR, member of the European Magnetic Field Laboratory (EMFL). We thank the
Nanofabrication Facilities Rossendorf team at the HZDR Ion Beam Center for the
support in samples fabrication. The authors acknowledge fruitful discussions
with Attila Kákay and Serhii Sorokin.

## References

  * Aharoni [1990] Amikam Aharoni.  Upper bound to a single‐domain behavior of a ferromagnetic cylinder.  _Journal of Applied Physics_ , 68(6):2892–2900, September 1990.  ISSN 0021-8979, 1089-7550.  doi: 10.1063/1.346422.  URL http://aip.scitation.org/doi/10.1063/1.346422.
  * Usov and Peschany [1993] N.A. Usov and S.E. Peschany.  Magnetization curling in a fine cylindrical particle.  _Journal of Magnetism and Magnetic Materials_ , 118(3), January 1993.  ISSN 03048853.  doi: 10.1016/0304-8853(93)90428-5.  URL https://linkinghub.elsevier.com/retrieve/pii/0304885393904285.
  * Shinjo [2000] T. Shinjo.  Magnetic Vortex Core Observation in Circular Dots of Permalloy.  _Science_ , 289(5481):930–932, August 2000.  doi: 10.1126/science.289.5481.930.  URL https://www.sciencemag.org/lookup/doi/10.1126/science.289.5481.930.
  * Kim et al. [2010a] Dong-Hyun Kim, Elena A. Rozhkova, Ilya V. Ulasov, Samuel D. Bader, Tijana Rajh, Maciej S. Lesniak, and Valentyn Novosad.  Biofunctionalized magnetic-vortex microdiscs for targeted cancer-cell destruction.  _Nature Materials_ , 9(2):165–171, February 2010a.  ISSN 1476-1122, 1476-4660.  doi: 10.1038/nmat2591.  URL http://www.nature.com/articles/nmat2591.
  * Leulmi et al. [2015] Selma Leulmi, Xavier Chauchet, Melissa Morcrette, Guillermo Ortiz, Hélène Joisten, Philippe Sabon, Thierry Livache, Yanxia Hou, Marie Carrière, Stéphane Lequien, and Bernard Dieny.  Triggering the apoptosis of targeted human renal cancer cells by the vibration of anisotropic magnetic particles attached to the cell membrane.  _Nanoscale_ , 7(38):15904–15914, 2015.  ISSN 2040-3364, 2040-3372.  doi: 10.1039/C5NR03518J.  URL http://xlink.rsc.org/?DOI=C5NR03518J.
  * Kim et al. [2011] Dong-Hyun Kim, Pavel Karavayev, Elena A. Rozhkova, John Pearson, Volodymyr Yefremenko, Samuel D. Bader, and Valentyn Novosad.  Mechanoresponsive system based on sub-micron chitosan-functionalized ferromagnetic disks.  _Journal of Materials Chemistry_ , 21:8422, 2011.  ISSN 0959-9428, 1364-5501.  doi: 10.1039/c1jm10272a.  URL http://xlink.rsc.org/?DOI=c1jm10272a.
  * Sluka et al. [2019] Volker Sluka, Tobias Schneider, Rodolfo A. Gallardo, Attila Kákay, Markus Weigand, Tobias Warnatz, Roland Mattheis, Alejandro Roldán-Molina, Pedro Landeros, Vasil Tiberkevich, Andrei Slavin, Gisela Schütz, Artur Erbe, Alina Deac, Jürgen Lindner, Jörg Raabe, Jürgen Fassbender, and Sebastian Wintz.  Emission and propagation of 1D and 2D spin waves with nanoscale wavelengths in anisotropic spin textures.  _Nature Nanotechnology_ , 14(4):328–333, April 2019.  ISSN 1748-3387, 1748-3395.  doi: 10.1038/s41565-019-0383-4.  URL http://www.nature.com/articles/s41565-019-0383-4.
  * Schultheiss et al. [2019] K. Schultheiss, R. Verba, F. Wehrmann, K. Wagner, L. Körber, T. Hula, T. Hache, A. Kákay, A.A. Awad, V. Tiberkevich, A.N. Slavin, J. Fassbender, and H. Schultheiss.  Excitation of Whispering Gallery Magnons in a Magnetic Vortex.  _Physical Review Letters_ , 122(9):097202, March 2019.  ISSN 0031-9007, 1079-7114.  doi: 10.1103/PhysRevLett.122.097202.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.097202.
  * Körber et al. [2020] L. Körber, K. Schultheiss, T. Hula, R. Verba, J. Fassbender, A. Kákay, and H. Schultheiss.  Nonlocal stimulation of three-magnon splitting in a magnetic vortex.  _Phys. Rev. Lett._ , 125:207203, Nov 2020.  doi: 10.1103/PhysRevLett.125.207203.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.207203.
  * Heins et al. [2022] Christopher Heins, Katrin Schultheiss, Lukas Körber, Attila Kakay, Tobias Hula, Mauricio Bejarano, Vadym Iurchuk, Jürgen Lindner, Jürgen Fassbender, and Helmut Schultheiss.  Spontaneous emergence of spin-wave frequency combs due to vortex gyration and spin-wave hybridisation.  7th Workshop on Magnonics, 2022.
  * Kim et al. [2013] J.-S. Kim, M. Kläui, M. V. Fistul, J. Yoon, C.-Y. You, R. Mattheis, C. Ulysse, and G. Faini.  Double resonance response in nonlinear magnetic vortex dynamics.  _Physical Review B_ , 88(6):064402, August 2013\.  ISSN 1098-0121, 1550-235X.  doi: 10.1103/PhysRevB.88.064402.  URL https://link.aps.org/doi/10.1103/PhysRevB.88.064402.
  * Stöhr and Siegmann [2006] Joachim Stöhr and Hans-Christoph Siegmann.  Magnetism: from fundamentals to nanoscale dynamics.  _Berlin, Springer series in solid-state sciences_ , 152:9783540302827, 2006.
  * Gonçalves et al. [2013] A. M. Gonçalves, I. Barsukov, Y.-J. Chen, L. Yang, J. A. Katine, and I. N. Krivorotov.  Spin torque ferromagnetic resonance with magnetic field modulation.  _Applied Physics Letters_ , 103(17):172406, October 2013.  ISSN 0003-6951, 1077-3118.  doi: 10.1063/1.4826927.  URL http://aip.scitation.org/doi/10.1063/1.4826927.
  * Ramasubramanian et al. [2020] Lakshmi Ramasubramanian, Attila Kákay, Ciarán Fowley, Oguz Yildirim, Patrick Matthes, Serhii Sorokin, Aleksandra Titova, Donovan Hilliard, Roman Böttger, René Hübner, Sibylle Gemming, Stefan E. Schulz, Florian Kronast, Denys Makarov, Jürgen Fassbender, and Alina Deac.  Tunable Magnetic Vortex Dynamics in Ion-Implanted Permalloy Disks.  _ACS Applied Materials & Interfaces_, 12(24):27812–27818, 2020.  doi: 10.1021/acsami.0c08024.  URL https://pubs.acs.org/doi/10.1021/acsami.0c08024.
  * Goto et al. [2011] M. Goto, H. Hata, A. Yamaguchi, Y. Nakatani, T. Yamaoka, and Y. Nozaki.  Electrical detection of vortex states in a ferromagnetic disk using the rectifying effect.  _Journal of Applied Physics_ , 109(7):07D306, April 2011.  ISSN 0021-8979, 1089-7550.  doi: 10.1063/1.3536390.  URL http://aip.scitation.org/doi/10.1063/1.3536390.
  * Buchanan et al. [2007] K. S. Buchanan, M. Grimsditch, F. Y. Fradin, S. D. Bader, and V. Novosad.  Driven Dynamic Mode Splitting of the Magnetic Vortex Translational Resonance.  _Physical Review Letters_ , 99(26):267201, December 2007.  ISSN 0031-9007, 1079-7114.  doi: 10.1103/PhysRevLett.99.267201.  URL https://link.aps.org/doi/10.1103/PhysRevLett.99.267201.
  * Guslienko et al. [2010] Konstantin Y. Guslienko, Rafael Hernández Heredero, and Oksana Chubykalo-Fesenko.  Nonlinear gyrotropic vortex dynamics in ferromagnetic dots.  _Physical Review B_ , 82(1):014402, July 2010\.  ISSN 1098-0121, 1550-235X.  doi: 10.1103/PhysRevB.82.014402.  URL https://link.aps.org/doi/10.1103/PhysRevB.82.014402.
  * Vansteenkiste et al. [2014] Arne Vansteenkiste, Jonathan Leliaert, Mykola Dvornik, Mathias Helsen, Felipe Garcia-Sanchez, and Bartel Van Waeyenberge.  The design and verification of MuMax3.  _AIP Advances_ , 4(10):107133, October 2014.  ISSN 2158-3226.  doi: 10.1063/1.4899186.  URL http://aip.scitation.org/doi/10.1063/1.4899186.
  * Iurchuk et al. [2021] Vadym Iurchuk, Lukas Körber, ALina Maria Deac, Jürgen Faßbender, Jürgen Lindner, and Attila Kákay.  Stress-induced modification of gyration dynamics in stacked double-vortex structures studied by micromagnetic simulations.  _Journal of Physics D: Applied Physics_ , 54(47):475002, November 2021.  ISSN 0022-3727, 1361-6463.  doi: 10.1088/1361-6463/ac2333.  URL https://iopscience.iop.org/article/10.1088/1361-6463/ac2333.
  * Arekapudi et al. [2021] Sri Sai Phani Kanth Arekapudi, Benny Böhm, Lakshmi Ramasubramanian, Fabian Ganss, Peter Heinig, Sven Stienen, Ciarán Fowley, Kilian Lenz, Alina M. Deac, Manfred Albrecht, and Olav Hellwig.  Direct imaging of distorted vortex structures and magnetic vortex annihilation processes in ferromagnetic/antiferromagnetic disk structures.  _Physical Review B_ , 103:014405, January 2021.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.103.014405.  URL https://link.aps.org/doi/10.1103/PhysRevB.103.014405.
  * Lee et al. [2008] Ki-Suk Lee, Sang-Koog Kim, Young-Sang Yu, Youn-Seok Choi, Konstantin Yu. Guslienko, Hyunsung Jung, and Peter Fischer.  Universal Criterion and Phase Diagram for Switching a Magnetic Vortex Core in Soft Magnetic Nanodots.  _Physical Review Letters_ , 101(26):267206, December 2008.  ISSN 0031-9007, 1079-7114.  doi: 10.1103/PhysRevLett.101.267206.  URL https://link.aps.org/doi/10.1103/PhysRevLett.101.267206.
  * Gliga et al. [2011] Sebastian Gliga, Yaowen Liu, and Riccardo Hertel.  Energy thresholds in the magnetic vortex core reversal.  _Journal of Physics: Conference Series_ , 303:012005, jul 2011.  doi: 10.1088/1742-6596/303/1/012005.  URL https://doi.org/10.1088/1742-6596/303/1/012005.
  * Kim et al. [2010b] J.-S. Kim, O. Boulle, S. Verstoep, L. Heyne, J. Rhensius, M. Kläui, L. J. Heyderman, F. Kronast, R. Mattheis, C. Ulysse, and G. Faini.  Current-induced vortex dynamics and pinning potentials probed by homodyne detection.  _Physical Review B_ , 82(10):104427, September 2010b.  ISSN 1098-0121, 1550-235X.  doi: 10.1103/PhysRevB.82.104427.  URL https://link.aps.org/doi/10.1103/PhysRevB.82.104427.
  * Sushruth et al. [2016] Manu Sushruth, Jasper P. Fried, Abdelmadjid Anane, Stephane Xavier, Cyrile Deranlot, Mikhail Kostylev, Vincent Cros, and Peter J. Metaxas.  Electrical measurement of magnetic-field-impeded polarity switching of a ferromagnetic vortex core.  _Physical Review B_ , 94(10):100402, September 2016.  ISSN 2469-9950, 2469-9969.  doi: 10.1103/PhysRevB.94.100402.  URL https://link.aps.org/doi/10.1103/PhysRevB.94.100402.