# Audio tagging with noisy labels and minimal supervision

###### Abstract

This paper introduces Task 2 of the DCASE2019 Challenge, titled “Audio tagging
with noisy labels and minimal supervision”. This task was hosted on the Kaggle
platform as “Freesound Audio Tagging 2019”. The task evaluates systems for
multi-label audio tagging using a large set of noisy-labeled data, and a much
smaller set of manually-labeled data, under a large vocabulary setting of 80
everyday sound classes. In addition, the proposed dataset poses an acoustic
mismatch problem between the noisy train set and the test set due to the fact
that they come from different web audio sources. This can correspond to a
realistic scenario given the difficulty of gathering large amounts of manually
labeled data. We present the task setup, the FSDKaggle2019 dataset prepared
for this scientific evaluation, and a baseline system consisting of a
convolutional neural network. All these resources are freely available.

Index Terms—  Audio tagging, sound event classification, audio dataset, label
noise, minimal supervision

##  1 Introduction

Environmental sound recognition has gained attention in recent years,
encompassing tasks such as acoustic scene classification, sound event
detection or audio tagging [1]. The latter is becoming a popular task partly
due to the various audio tagging tasks in the DCASE Challenge editions, and
its impact on applications such as automatic description of multimedia, or
acoustic monitoring. This paper describes the characteristics, dataset and
baseline system of DCASE2019 Task 2 “Audio tagging with noisy labels and
minimal supervision”.

Everyday sound tagging consists of identifying the sound events present in an
audio recording. The most common approach to create sound event taggers relies
on supervised learning through labeled audio datasets. New released datasets
tend to be of increasing size in order to allow exploitation of data-driven
approaches, e.g., deep learning. However, manual labeling of large datasets is
expensive and typically a limiting factor in machine listening; hence,
creators are often forced to compromise between dataset size and label
quality. Thus, most recent datasets feature larger sizes [2, 3, 4], but their
labeling is less precise than that of conventional small and exhaustively
labeled datasets [5, 6, 7]. We are, therefore, witnessing a transition towards
larger datasets that inevitably include some degree of label noise. Likewise,
the current trend is moving towards general-purpose sound event recognizers,
able to recognize a broad range of everyday sounds. This is favoured by the
appearance of the AudioSet Ontology; a hierarchical tree with 627 classes
encompassing the most common everyday sounds [2]. This implies that not only
are we interested in recognizing typical sound sources (e.g., _Bark_), but
also less usual sound classes (e.g., production mechanisms such as _Fill (with
liquid)_). Further, not only are these classes less frequent, but also some
can be semantically or acoustically similar to others (e.g. _Trickle, dribble_
and _Fill (with liquid)_). Manual annotation of these more ambiguous
categories becomes more difficult, which makes them more prone to labeling
errors. An alternative to gathering data for training general-purpose audio
taggers is to retrieve audio and metadata from websites such as Freesound or
Flickr. Labels can be inferred automatically by using automated heuristics
applied to the metadata, or applying pre-trained classifiers on the audio
material. This approach supports rapid collection of large amounts of data,
but at the cost of a considerable amount of label noise.

In this context, label noise arises as a challenge in general-purpose sound
event recognition, including adverse effects such as performance drop or
increased complexity of models [8], and also hindering proper learning of deep
networks [9, 10]. Consequently, coping with label noise could open the door to
better sound event classifiers, and could allow the exploitation of large
amounts of web audio for training, while reducing manual annotation needs. The
topic of learning with noisy labels is a consolidated research area in
computer vision [11, 12, 13, 14, 15, 16, 17]. However, in sound recognition it
has received little attention, probably due to the conventional paradigm of
learning from small and clean datasets; only a few works directly address the
analysis and mitigation of label noise [4, 18, 19].

In this paper, we propose a task, a dataset and a baseline system to foster
label noise research in general-purpose sound event tagging. We follow up on
DCASE2018 Task 2 [3], and propose to investigate the scenario where a small
set of manually-labeled data is available, along with a larger set of noisy-
labeled data, in a multi-label audio tagging setting, and using a vocabulary
of 80 classes of everyday sounds. The proposed task addresses two main
research problems. The first problem is how to adequately exploit a large
quantity of noisy labels, many of which are incorrect and/or incomplete, and
how to complement it with the supervision provided by a much smaller amount of
reliable manually-labeled data (minimal supervision). The second problem is
given by the acoustic mismatch between the noisy train set and the test set.
Distribution shifts between data have been shown to cause substantial
performance drops in machine learning, both for vision [20] and audio [21]. In
our case, the noisy train set comes from a different web audio source than the
test set, which is sometimes a real-world constraint. This paper is organized
as follows. Section 2 provides more details about the task and its
experimental setup. Section 3 presents the dataset prepared for the task, and
Section 4 describes a baseline system. Final remarks are given in Section 5.

##  2 Task Setup

The goal of this task is to predict appropriate labels for each audio clip in
a test set. The predictions are to be done at the clip level, i.e., no
start/end timestamps for the sound events are required. Some test clips bear
one ground truth label while others bear several labels. Hence, the task setup
is a _multi-label_ classification problem and the systems to be developed can
be denoted as multi-label audio tagging systems, as illustrated in Fig. 1.
This task was hosted on the Kaggle platform from April 4th to June 10th 2019.
The resources associated to this task (dataset download, submission, and
leaderboard) can be found on the Kaggle competition
page.111https://www.kaggle.com/c/freesound-audio-tagging-2019
Note that the competition name on Kaggle is abbreviated from the full
DCASE2019 Challenge task name to “Freesound Audio Tagging 2019”.

Figure 1: Overview of a multi-label tagging system.

As described in Section 3, the audio data for this task consists of a test set
and two train sets: a _curated_ train set and a _noisy_ train set, that allow
to experiment with training data of different levels of reliability and coming
from different sources. System evaluation was carried out on Kaggle servers
(see Section 2.1) using the test set, which is further split into two
divisions, for the public and private leaderboards. During the competition,
the test subset corresponding to the public leaderboard was used to provide
live ranking of all participants. To compute the final private leaderboard, at
the end of the competition, systems were re-evaluated using the unseen private
test set, of which neither the audio nor the labels were accessible to
participants.

###  2.1 Evaluation Metric and Competition Rules

The task was evaluated via label-weighted label-ranking average precision
(abbreviated as _l $\omega$lrap_ and pronounced “lol wrap”). This generalizes
the mean reciprocal rank (MRR) used in the 2018 challenge [3] to the case of
multiple true labels per test item. Let $Lab(s,r)$ be the class label at rank
$r$ (starting from 1) in test sample $s$, and $Rank(s,c)$ be the rank of class
label $c$ in that list, i.e. $Lab(s,Rank(s,c))=c$. Then, if the set of ground-
truth classes for sample $s$ is $C(s)$, the label-ranking precision for the
list of labels up to class $c$ (assumed to be in $C(s)$) is:

$Prec(s,c)=\frac{1}{Rank(s,c)}\sum_{r=1}^{Rank(s,c)}\mathbf{1}[Lab(s,r)\in
C(s)]$ (1)

where $\mathbf{1}[\cdot]$ evaluates to 1 if the argument is true, else zero.
$Prec(s,c)$ is equal to 1 if all the top-ranked labels down to $c$ are part of
$C(s)$, and at worst case equals $1/Rank(s,c)$ if none of the higher-ranked
labels are correct. In contrast to plain _lrap_ , which averages precisions
within a sample then across samples, thereby downweighting labels that occur
on samples with many labels, _l $\omega$lrap_ calculates the precision for
each _label_ in the test set, and gives them all equal contribution to the
final metric:

$l\omega lrap=\frac{1}{\sum_{s}|C(s)|}\sum_{s}\sum_{c\in C(s)}Prec(s,c)$ (2)

where $|C(s)|$ is the number of true class labels for sample $s$. We use label
weighting because it allows per-class values to be calculated, while keeping
the overall metric as a simple average of the per-class metrics (weighted by
each label’s prior in the test set). A Python implementation of _l
$\omega$lrap_ is provided in
222https://colab.research.google.com/drive/1AgPdhSp7ttY18O3fEoHOQKlt_3HJDLi8.

This scientific evaluation was set up as a Kaggle _Kernels-only_ competition.
This means that all participants had to submit their systems as inference
models in Kaggle Kernels (similar to _Jupyter Notebooks_), to be evaluated on
remote servers. In addition, inference run-time was limited to a maximum of
one hour in a Kernel with one GPU, and memory constraints were also imposed.
These constraints aim to discourage the usage of large model ensembles.
Participants could submit a maximum of two submissions per day, and select two
_final_ submissions to be considered for the private leaderboard ranking. A
detailed description of the task rules can be found in the Rules section of
the competition page;1 the most important points are summarized in the DCASE
Challenge page.333http://dcase.community/challenge2019/task-audio-
tagging#task-rules

To complement the leaderboard results of the _l $\omega$lrap_ ranking, the
task organizers introduced a complementary Judges’ Award to promote
submissions using novel, problem-specific and efficient approaches. Details
about the Judges’ Award rules can be found in the Judges’ Award section of 1.

##  3 Dataset

The dataset used is called _FSDKaggle2019_ , and it employs audio clips from
the following sources:

  * •

Freesound Dataset (FSD): a dataset under development based on Freesound
content organized with the AudioSet Ontology [2]. Freesound is a sound sharing
site hosting over 400,000 clips uploaded by a community of users, who
additionally provide some basic metadata, e.g., tags. [22]. These data are
used to create the curated train set and the test set.

  * •

The soundtracks of a pool of Flickr videos taken from the Yahoo Flickr
Creative Commons 100M (YFCC100M) dataset [23]. These data are used to create
the noisy train set.

FSDKaggle2019 is freely available from
Zenodo444https://doi.org/10.5281/zenodo.3612637, all clips are provided as
uncompressed PCM 16 bit 44.1 kHz mono audio files, its ground truth labels are
provided at the clip-level (i.e., weak labels), and its partitioning is
depicted in Fig. 2.

Figure 2: Data split in FSDKaggle2019, including number of clips / duration in
hours, and data origin. Colors depict quality of labels: orange, yellow and
green correspond to noisy labels, correct but potentially incomplete labels,
and exhaustive labels, respectively.

###  3.1 Curated train set and test set

The first step carried out in the creation of FSDKaggle2019 was the definition
of a vocabulary of 80 classes drawn from the AudioSet Ontology [2]. This
vocabulary was chosen based on the following criteria: i) we consider leaf
nodes of the AudioSet hierarchy for which there is enough data available in
FSD, ii) we aim to encompass a diverse range of everyday sounds, and iii) we
remove few clearly isolated classes (those with the weakest semantic relations
with any of the rest), thus promoting confounds between
semantically/acoustically similar classes to some extent. The main sound
families (i.e., groups of sound classes) in the resulting vocabulary are, in
descending order of prevalence, human sounds, domestic sounds, musical
instruments, vehicles, animal sounds, natural sounds, materials, and
mechanisms. The full list of 80 classes is available in the Data section of 1.

In a second step, we did a mapping of Freesound clips to the selected 80 class
labels. To this end, a set of keywords was defined connecting the user-
provided Freesound tags with the AudioSet labels. Using this mapping, for
every class, we retrieved the audio clips that feature at least one of the
defined keywords among their tags. This process led to a number of
automatically-generated candidate annotations indicating the potential
presence of a sound class in an audio clip (i.e., weak labels, as timing
information is not included). Nonetheless, in some audio clips the target
signal fills the clip length almost completely, which can be considered as a
strong label. Subsequently, the candidate annotations were human-validated
using a validation task deployed in _Freesound Annotator_
,555https://annotator.freesound.org an online platform for the collaborative
creation of open audio datasets [24]. In this task, users verify the
presence/absence of a candidate sound class in an audio clip with a _rating_
mechanism. The vast majority of provided labels have inter-annotator agreement
but not all of them. The outcome is a set of clips where the corresponding
label(s) are correct; nevertheless, it can happen that a few of these audio
clips present additional acoustic material beyond the provided label(s).

The resulting data were split into a train set and a test set. We refer to
this train set as curated in order to distinguish it from the noisy set
described in Section 3.2. To mitigate train-test contamination, the split was
carried out considering the clip uploaders in Freesound. We allocated all
audio clips uploaded from the same user into either the curated train set or
the test set, so that the sets are disjoint at the Freesound user level. The
partition proportion was defined to limit the supervision provided in the
curated train set, thus promoting approaches to deal with label noise.

Finally, the test set was further annotated using a label generation tool
[25], in which i) pre-existent labels can be re-validated, and ii) potentially
missing labels can be added through exploration of the AudioSet Ontology. The
outcome is a set of exhaustively labeled clips where the label(s) are correct
and complete considering the target vocabulary; nonetheless, few clips could
still present additional (unlabeled) acoustic content out of the vocabulary.

The main characteristics of the curated train set, noisy train set and test
set are listed in Table 1. The curated train set consists of 4970 clips with a
total of 5752 labels. Labels per clip ranges from 1 to 6 with a mean of 1.2.
The test set consists of 4481 clips with a total of 6250 labels. Labels per
clip ranges from 1 to 6 with a mean of 1.4. Note the increased number of
labels per clip with respect to the curated train set, due to the process of
exhaustive labelling. In both cases, clip length ranges from 0.3s to 30 due to
the diversity of the sound classes and the preferences of Freesound users when
recording/uploading sounds.

Table 1: Main stats of the sets in FSDKaggle2019. ∗A few classes have slightly less than 75 clips. Aspect | curated train | noisy train | test
---|---|---|---
Clips/class |  $\sim$75∗ | 300 |  $\sim$ 50 - 150
Total clips | 4970 | 19,815 | 4481
Labels/clip | 1.2 | 1.2 | 1.4
Clip length |  $\sim$0.3 - 30s |  $\sim$15s |  $\sim$0.3 - 30s
Total duration |  $\sim$10.5h |  $\sim$80h |  $\sim$12.9h
Labelling | correct | noisy | exhaustive
| (inexhaustive) |  |

###  3.2 Noisy train set

The noisy train set was prepared using the YFCC100M dataset [23], which has
the advantages of i) being a very large and diverse dataset that is not
correlated with Freesound in acoustics or domain, and ii) offering permissive
Creative Commons licenses that allow ease of use, modification, and
redistribution. The original dataset contained $\sim$99M photos and $\sim$793k
videos from $\sim$581K Flickr users. We dropped videos with licenses that
disallowed making derivatives or commercial use, videos that were no longer
available, and videos with audio decode errors that we could not transcode,
leaving us with $\sim$201K 44.1 kHz mono WAV files. Video length varied with a
maximum of 20 minutes, and a mean of $\sim$37s and median of $\sim$20s.

The Flickr video metadata (title, description, tags) proved to be too sparse
to meaningfully map to our class vocabulary. Therefore, we used a content-
based approach where we generated video-level predictions from a variety of
pre-trained audio models: a shallow fully-connected network as well as
variants of VGG and ResNet [26], all of which were trained on a large
collection of YouTube videos using the AudioSet class vocabulary. We generated
sliding windows of $\sim$1s containing log mel spectrogram patches and
aggregated the per-window predictions (using either maximum or average
pooling) to produce a video-level vector of class scores. For each of our 80
classes, we kept the top 300 videos by predicted score for that class. We
browsed the video labels and selected the maximum-pooled VGG-like model as
producing a balance between reasonable predictions and a substantial amount of
label noise. As a further source of noise, each final clip was produced by
taking a random slice of a video of length up to 15 seconds (videos shorter
than 15 seconds would be taken in their entirety). Hence, the label noise can
vary widely in amount and type depending on the class, including in- and out-
of-vocabulary noises [4].

As listed in Table 1, the noisy train set consists of 19,815 clips with a
total of 24,000 labels (300 * 80). Labels per clip ranges from 1 to 7 with a
mean of 1.2. Clip length ranges from 1s to 15s (by construction), with a mean
of 14.5s. Therefore, the per-class training data distribution in FSDKaggle2019
is, for most of the classes, 300 clips from the noisy set and 75 clips from
the curated set. This means 80% noisy / 20% curated at the clip level, while
at the duration level the proportion is more extreme considering the variable-
length clips. Since most of the train data come from YFCC, acoustic domain
mismatch between the train and test set can be expected. We conjecture this
mismatch comes from a variety of reasons. For example, through acoustic
inspection of a small sample of both data sources, we find a higher percentage
of high quality recordings in Freesound. In addition, audio clips in Freesound
are typically recorded with the purpose of capturing audio, which is not
necessarily the case in YFCC.

##  4 Baseline System

###  4.1 Model Architecture and Training

The baseline system uses a convolutional neural network that takes log mel
spectrogram patches as input and produces predicted scores for 80 classes. We
use an efficient MobileNet v1 [27] architecture which lets us fit comfortably
within the inference time limits of the challenge. Incoming audio (always 44.1
kHz mono) is divided into overlapping windows of size 1s with a hop of 0.5s.
These windows are decomposed with a short-time Fourier transform using 25ms
windows every 10ms. The resulting spectrogram is mapped into 96 mel-spaced
frequency bins covering 20 Hz to 20 kHz, and the magnitude of each bin is log-
transformed after adding a small offset to avoid numerical issues. The model
consists of 1 convolutional layer followed by 13 separable convolution layers
(which give MobileNets their compactness) followed by either max or average
pooling, and an 80-way logistic classifier layer. The model contains
$\sim$3.3M weights and, by comparison, is $\sim$8x smaller than a ResNet-50
while using $\sim$4x less compute. Detailed documentation is available in the
public release of the baseline system code.666https://github.com/DCASE-
REPO/dcase2019_task2_baseline

To combat label noise, we use dropout (inserted before the classifier layer)
as well as label smoothing [28] which replaces each label’s target value with
a hyperparameter-controlled blend of the original value and 0.5 (representing
a uniform probability distribution). To combat the domain mismatch between the
test set and noisy train set, we use transfer learning by training a model on
the noisy set first to learn a representation, and then use a checkpoint from
that run to warm-start training on the curated train set. In addition, we used
batch normalization, exponential learning rate decay, and the Adam optimizer.
We trained models on the curated data alone, the noisy data alone, the curated
and noisy combined, noisy first followed by warm-started curated, as well as a
weighted version of warm-started training that we describe next.

###  4.2 Results

Table 2 shows the _l $\omega$lrap_ values produced by our various baseline
models when evaluated on the entire test set (i.e., including both the public
and private splits). Each row lists the best _l $\omega$lrap_ obtained from a
small grid search (using the public test set for evaluation) that varied
maximum vs average pooling, learning rate, label smoothing, dropout, learning
rate decay, and whether or not we used batch normalization. The baseline
system that we settled on was “Warm-started curated”, which achieved a _l
$\omega$lrap_ of 0.537 on the public test set (see the publicly released
baseline code for hyperparameter choices).

Table 2: Baseline system results on the entire test set. Training approach | _l $\omega$lrap_
---|---
Curated only | 0.542
Noisy only | 0.312
Curated + Noisy | 0.522
Warm-started curated | 0.546
Weighted warm-started curated | 0.561

Comparing “Noisy only” and “Curated + Noisy” to “Curated only” shows a
considerable domain mismatch where we hurt our performance when we blindly add
more data. A transfer learning approach of warm-starting the curated training
with a noisily trained model gives us a small boost in performance.

We conducted a class-based analysis of the results in a public Colab
notebook777https://colab.research.google.com/drive/1vIpAWCGSHyYYLblbtMzAKLNTkHcr1SMK
where we look at the best and worst classes of each model. The baseline system
attains highest _l $\omega$lrap_ values for _Bicycle bell_ (0.894) and _Purr_
(0.873); lowest values occur for _Cupboard open or close_ (0.219) and _Chirp,
tweet_ (0.127).

We also analyse the correlations between various pairs of models. It becomes
evident that, at least for our baseline models, there are classes where using
curated data alone is better while there are other classes where the noisy
data is better. One simple way to incorporate this in training is to use the
ratio of noisy to curated _l $\omega$lrap_s as a per-class weight during noisy
training to boost classes that have value and suppress classes that do not.
When we warmstart with this weighted noisy model, we get a further boost in
performance. This optimization is not included in the released baseline.

##  5 Conclusion

In this paper, we have described the task setup, dataset, and baseline system
of DCASE2019 Task 2 “Audio tagging with noisy labels and minimal supervision”.
This task was hosted on the Kaggle platform as “Freesound Audio Tagging 2019”.
The goal is to adequately exploit a large set of noisy-labeled data and a
small quantity of manually-labeled data, in a multi-label audio tagging
setting with a vocabulary of 80 everyday sound classes. In addition, the
dataset poses an acoustic mismatch problem between the noisy train set and the
test set due to the fact that they come from different web audio sources. We
believe this can correspond to a realistic scenario given the difficulty in
gathering large amounts of manually labeled data. _l $\omega$lrap_ is proposed
as evaluation metric. Baseline experiments indicate that leveraging noisy-
labeled data with a distribution shift for sound event tagging can be
challenging. The FSDKaggle2019 dataset and the baseline system proposed are
freely available and not limited for use within the competition.

##  6 ACKNOWLEDGMENT

We thank Addison Howard and Walter Reade of Kaggle for their invaluable
assistance with the task design and Kaggle platform, and everyone who
contributed to FSDKaggle2019 with annotations, especially Jordi Pons and
Xavier Favory. Eduardo Fonseca is also grateful for the GPU donated by NVidia.

## References

  * [1] T. Virtanen, M. D. Plumbley, and D. Ellis, _Computational Analysis of Sound Scenes and Events_. Springer, 2018\.
  * [2] J. F. Gemmeke, D. P. W. Ellis, D. Freedman, A. Jansen, W. Lawrence, R. C. Moore, M. Plakal, and M. Ritter, “Audio Set: An ontology and human-labeled dataset for audio events,” in _Proc. IEEE ICASSP 2017_ , New Orleans, LA, 2017.
  * [3] E. Fonseca, M. Plakal, F. Font, D. P. Ellis, X. Favory, J. Pons, and X. Serra, “General-purpose tagging of Freesound audio with AudioSet labels: task description, dataset, and baseline,” in _Proceedings of the Detection and Classification of Acoustic Scenes and Events 2018 Workshop (DCASE2018)_ , 2018\.
  * [4] E. Fonseca, M. Plakal, D. P. W. Ellis, F. Font, X. Favory, and X. Serra, “Learning Sound Event Classifiers from Web Audio with Noisy Labels,” in _Proc. IEEE ICASSP 2019_ , Brighton, UK, 2019.
  * [5] J. Salamon, C. Jacoby, and J. P. Bello, “A dataset and taxonomy for urban sound research,” in _Proceedings of the ACM International Conference on Multimedia_. ACM, 2014, pp. 1041–1044.
  * [6] K. J. Piczak, “ESC: Dataset for environmental sound classification,” in _Proceedings of the ACM International Conference on Multimedia_. ACM, 2015, pp. 1015–1018.
  * [7] P. Foster, S. Sigtia, S. Krstulovic, J. Barker, and M. D. Plumbley, “CHiME-home: A dataset for sound source recognition in a domestic environment,” in _Workshop on Applications of Signal Processing to Audio and Acoustics_. IEEE, 2015.
  * [8] B. Frénay and M. Verleysen, “Classification in the presence of label noise: a survey,” _IEEE transactions on neural networks and learning systems_ , vol. 25, no. 5, 2014.
  * [9] D. Arpit, S. Jastrzebski, N. Ballas, D. Krueger, E. Bengio, M. S. Kanwal, T. Maharaj, A. Fischer, A. Courville, Y. Bengio, _et al._ , “A closer look at memorization in deep networks,” in _Proceedings of the 34th International Conference on Machine Learning-Volume 70_. JMLR. org, 2017, pp. 233–242.
  * [10] C. Zhang, S. Bengio, M. Hardt, B. Recht, and O. Vinyals, “Understanding deep learning requires rethinking generalization,” _arXiv preprint arXiv:1611.03530_ , 2016.
  * [11] B. Han, Q. Yao, X. Yu, G. Niu, M. Xu, W. Hu, I. Tsang, and M. Sugiyama, “Co-teaching: Robust training of deep neural networks with extremely noisy labels,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 8527–8537.
  * [12] E. Malach and S. Shalev-Shwartz, “Decoupling” when to update” from ”how to update”,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 960–970.
  * [13] L. Jiang, Z. Zhou, T. Leung, L.-J. Li, and L. Fei-Fei, “Mentornet: Learning data-driven curriculum for very deep neural networks on corrupted labels,” _arXiv preprint arXiv:1712.05055_ , 2017.
  * [14] J. Goldberger and E. Ben-Reuven, “Training deep neural-networks using a noise adaptation layer,” in _Proceedings of the International Conference on Learning Representations_ , 2017.
  * [15] A. Veit, N. Alldrin, G. Chechik, I. Krasin, A. Gupta, and S. J. Belongie, “Learning from noisy large-scale datasets with minimal supervision.” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 6575–6583.
  * [16] Z. Zhang and M. Sabuncu, “Generalized cross entropy loss for training deep neural networks with noisy labels,” in _Advances in Neural Information Processing Systems_ , 2018.
  * [17] G. Patrini, A. Rozza, A. Krishna Menon, R. Nock, and L. Qu, “Making deep neural networks robust to label noise: A loss correction approach,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 1944–1952.
  * [18] A. Kumar, A. Shah, A. Hauptmann, and B. Raj, “Learning sound events from webly labeled data,” _arXiv preprint arXiv:1811.09967_ , 2018.
  * [19] E. Fonseca, F. Font, and X. Serra, “Model-agnostic approaches to handling noisy labels when training sound event classifiers,” in _Proceedings of IEEE Workshop on Applications of Signal Processing to Audio and Acoustics_ , New York, US, 2019.
  * [20] B. Recht, R. Roelofs, L. Schmidt, and V. Shankar, “Do ImageNet Classifiers Generalize to ImageNet?” _arXiv preprint arXiv:1902.10811_ , 2019.
  * [21] A. Mesaros, T. Heittola, and T. Virtanen, “Acoustic scene classification: an overview of DCASE 2017 challenge entries,” in _2018 16th International Workshop on Acoustic Signal Enhancement (IWAENC)_. IEEE, 2018, pp. 411–415.
  * [22] F. Font, G. Roma, and X. Serra, “Freesound technical demo,” in _Proceedings of the ACM International Conference on Multimedia_. ACM, 2013, pp. 411–412.
  * [23] B. Thomee, D. A. Shamma, G. Friedland, B. Elizalde, K. Ni, D. Poland, D. Borth, and L.-J. Li, “YFCC100M: The New Data in Multimedia Research,” _Communications of the ACM_ , vol. 59, no. 2, pp. 64–73, Feb. 2016.
  * [24] E. Fonseca, J. Pons, X. Favory, F. Font, D. Bogdanov, A. Ferraro, S. Oramas, A. Porter, and X. Serra, “Freesound Datasets: a platform for the creation of open audio datasets,” in _Proceedings of the 18th International Society for Music Information Retrieval Conference (ISMIR 2017)_ , Suzhou, China, 2017, pp. 486–493.
  * [25] X. Favory, E. Fonseca, F. Font, and X. Serra, “Facilitating the manual annotation of sounds when using large taxonomies,” in _Proceedings of the 23rd Conference of Open Innovations Association FRUCT_. FRUCT Oy, 2018, p. 60.
  * [26] S. Hershey, S. Chaudhuri, D. P. W. Ellis, J. F. Gemmeke, A. Jansen, C. Moore, M. Plakal, D. Platt, R. A. Saurous, B. Seybold, M. Slaney, R. Weiss, and K. Wilson, “CNN Architectures for Large-Scale Audio Classification,” in _Proc. IEEE ICASSP 2017_ , New Orleans, LA, 2017.
  * [27] A. G. Howard, M. Zhu, B. Chen, D. Kalenichenko, W. Wang, T. Weyand, M. Andreetto, and H. Adam, “MobileNets: Efficient Convolutional Neural Networks for Mobile Vision Applications,” _arXiv preprint arXiv:1704.04861_ , 2017.
  * [28] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the Inception Architecture for Computer Vision,” in _Proceedings of IEEE Conference on Computer Vision and Pattern Recognition,_ , 2016. [Online]. Available: http://arxiv.org/abs/1512.00567


# Stable vortex in Bose-Einstein condensate dark matter

Y.O. Nikolaieva1, A.O. Olashyn1,2, Y.I. Kuriatnikov1,3, S.I. Vilchynskii1,4,
and A.I. Yakimenko1 1 Department of Physics, Taras Shevchenko National
University of Kyiv, 64/13, Volodymyrska Street, Kyiv 01601, Ukraine  2
Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The
Netherlands  3 Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria  4
Institute of Physics, Laboratory for Particle 6 Physics and Cosmology (LPPC),
École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

###### Abstract

The nature of dark matter (DM) is one of the most fascinating unresolved
challenges of modern physics. One of the perspective hypotheses suggests that
DM consists of ultralight bosonic particles in the state of Bose-Einstein
condensate (BEC). The superfluid nature of BEC must dramatically affect the
properties of DM matter including quantization of the angular momentum.
Angular momentum quantum in the form of a vortex line is expected to produce a
considerable impact on the luminous matter in galaxies including density
distribution and rotation curves. We investigate the evolution of spinning DM
cloud with typical galactic halo mass and radius. Analytically and numerically
stationary vortex soliton states with different topological charges have been
analyzed. It has been shown that while all multi-charged vortex states are
unstable, a single-charged vortex soliton is extremely robust and survives
during the lifetime of the Universe.

##  I Introduction

Most of the recent cosmological models for galaxy structure are amount to show
the galaxy as a luminous galactic baryon disk surrounded by a spherical
galactic halo of so-called Dark Matter (DM). Different estimations (such as
direct detection, gamma ray detection, microlensing measuring) [1, 2, 3, 4, 5,
6] give DM roughly 95$\%$ of the total mass of the galaxy. The nature of DM
remains one of the most exciting open questions in modern physics. Weakly
Interacting Massive Particles (WIMPs) with a mass of O(100) GeV have been one
of the leading DM candidates for a long time [7]. However, recent negative
results for indirect detection [8] and collider experiments [9] cause strong
motivation for developing alternative DM models. Plenty of theoretical models
have shown that considering sub-GeV DM has some advantages. One of such
theories suggests considering DM as ultralight bosonic particles in the state
of Bose-Einstein condensate (BEC) [10]. These bosonic particles can interact
via gravity and probably via weak interaction. Both of these interactions are
extremely insignificant (for instance, a mass of axions (one of the possible
candidates for DM) is estimated to be in the range $10^{-25}-10^{-2}$ eV [11];
$s$-wave scattering length corresponds to the two-particle interaction varies
in the very wide range $10^{-100}-10^{-1}$ fm). However, at the galactic and
astrophysical scales, the self-gravitating BEC may form stable structures in
the form of the galactic halo and astrophysical cold dark matter (CDM)
structures (Bose-stars) [12]. Despite the extremely small mass of these
bosonic particles, their gravity force dominates the Universe.

The natural question arises: what is the physical mechanism that stands behind
the bosonic DM self-stabilisation? There are two most probable mechanisms to
compensate the gravitational self-attraction and prevent collapse: (i) The
quantum pressure that occurs whenever the condensate density is inhomogeneous.
(ii) Interparticle repulsive interaction, which in mean-field approach leads
to nonlinear self-induced potential proportional to condensate density. Each
of these two mechanisms can stabilize the self-gravitating BEC and may lead to
a formation of soliton-like stationary in time spatial structures. These
nonlinear self-organized structures are well known in various physical
systems.

The hypotheses that dark matter structures of astrophysical and galactic
scales can be treated as self-gravitating BEC composed of extremely light
bosonic particles have been developed for decades. Bose stars as lumps of
Bose-Einstein condensates bound by self-gravity were proposed over 50 years
ago [13, 14]. The stability of such objects has been studied previously
numerically in a non-relativistic regime [15]. Also, the formation of Bose
stars [16] and their collapse [17] have been already studied. The collisional
dynamics of stable solitary waves in the Schrödinger–Poisson equation have
been discussed in Ref.[18]. Also have been analytically studied basic
properties of self-gravitating BEC in a harmonic trap with Hartree–Fock method
and compared with numerical calculations [19].

The superfluid nature of BEC can dramatically affect the properties of DM
including the formation of quantum vortices and quantization of angular
momentum. There have been studied possible eﬀects of subgalactic vortices in
the DM on the rotation velocity curves of virialized galaxies with standard DM
halo proﬁles [20]. In Ref.[21] exact solution for a single axisymmetric vortex
has been found analytically in the Thomas-Fermi regime. Conditions of vortex
formation in galactic halos composed of BEC DM have been discussed
analytically in Ref.[22]. The case of rigid rotation and its impact on BEC DM
with and without self-interaction has been examined in Ref.[23]. Rigid slow
rotation of BEC DM has been investigated analytically in Tomas-Fermi limit in
Ref.[24]. Also, stability and dynamical properties of slowly rotating
gravitationally self-bound BEC have been studied in [25]. It is remarkable
that gravity-like attractive nonlocal interaction has been extensively studied
in the context of atomic BEC (see recent review article [26]). Similar to
nonlocal optical media [27, 28] and BECs with long-range dipole-dipole [29]
interactions stable spinning solitons and azumthons have been predicted in
atomic BECs with gravity-like attractive interactions [30].

In the present work, we consider CDM of galactic scales. Our main goal is a
consistent analysis of vortex structures in self-gravitating BECs. We address
the following questions: (i) is it possible to balance such condensate in a
state with nonzero angular momentum? (ii) is such a spinning superfluid CDM
halo stable, and (iii) how the vortex structure manifests itself in observable
properties of the luminous matter of the galaxy?

Our paper is organized as follows. In Sec.II we define the model which we use
to investigate the system. In Sec.III we investigate the stationary solutions
in two approaches: using variational analysis and using numerical modelling.
In Sec.IV the dynamics of the vortex structures in self-gravitating BEC with
typical galactic mass is investigated numerically. We summarize our results in
the concluding Sec.V.

##  II Model

At the zero-temperature limit, all the bosons condense into the same quantum
ground state and the system is described by single condensate wave function
$\Psi(\mathbf{r},t)$. In the mean-field approximation, the dynamics of self-
gravitating BEC of $N$ weakly interacting bosons with mass $m$ is described by
the Gross-Pitaevskii-Poisson (GPP) system of equations [31]:

$\displaystyle i\hbar\frac{\partial\Psi}{\partial
t}=\left(-\frac{\hbar^{2}}{2m}\nabla^{2}+gN|\Psi|^{2}+m\Phi\right)\Psi$
$\displaystyle\nabla^{2}\Phi=4\pi GmN|\Psi|^{2},$ (1)

where $\Psi(\mathbf{r},t)$ is a complex wave function of the condensate with
normalization condition $\int|\Psi|^{2}d\textbf{r}=1$,
$g=4\pi\hbar^{2}a_{s}/m$ is the coupling strength that corresponds to the
repulsive two-particle interaction, $a_{s}$ is the $s$-wave scattering length,
$\mathbf{r}=(x,y,z)$ \- spacial coordinates, $t$ is time, $\Phi(\mathbf{r},t)$
is the gravitational potential, $G$ is the gravitational constant. The density
field can be written as:

$\rho=mN|\Psi|^{2},$ (2)

where $mN=M$ is the total mass of the galaxy’s halo. The total energy
associated with the GPP system can be written as

$E=\Theta+U+W,$ (3)

the kinetic energy

$\Theta=\frac{N\hbar^{2}}{2m}\int|\nabla\Psi|^{2}d\mathbf{r},$ (4)

the internal energy

$U=\frac{2\pi a_{s}\hbar^{2}}{m^{3}}\int\rho^{2}d\mathbf{r},$ (5)

and the gravitational potential energy of interaction

$W=\frac{1}{2}\int\rho\Phi d\mathbf{r}.$ (6)

In terms of dimensionless units
$[\mathbf{r}\rightarrow\mathbf{r}/L_{*},t\rightarrow\Omega_{*}t,E\rightarrow
E/\epsilon]$ the system Eq.(II) can be written in dimensionless form,

$\displaystyle i\frac{\partial\Psi(\mathbf{r},t)}{\partial
t}=\left(-\frac{1}{2}\nabla^{2}+\Phi(\mathbf{r},t)+|\Psi(\mathbf{r},t)|^{2}\right)\Psi(\mathbf{r},t),$
$\displaystyle\nabla^{2}\Phi(\mathbf{r},t)=|\Psi(\mathbf{r},t)|^{2},$ (7)

where $L_{*}=\lambda_{C}(m_{\textrm{Pl}}/m)\sqrt{\lambda/8\pi}$,
$\Omega_{*}=c\lambda_{C}/L_{*}^{2}$, $\epsilon=(\hbar^{2}/4\pi
m_{Pl}\lambda_{C}^{2})\left(8\pi/\lambda\right)^{3/2}$,
$m_{\textrm{Pl}}=\sqrt{\hbar c/G}$ is the Planck mass,
$\lambda/8\pi=a_{s}/\lambda_{C}$ is the self-interaction constant, and
$\lambda_{C}=\hbar/mc$ is the Compton wavelength of the bosons. With the new
dimensionless gravitational potential
$\Phi\rightarrow(L_{*}/\lambda_{C})^{2}\Phi/c^{2}$ and the wave function
$\Psi\rightarrow(\lambda/8\pi)\left(m_{\textrm{Pl}}/m\right)^{2}\sqrt{4\pi
GM}(\hbar/mc^{2})\Psi$. From now on in our paper, we use dimensionless
variables.

Finally, the normalization condition in dimensionless units:

$\int|\Psi|^{2}d\textbf{r}=4\pi\frac{M}{m_{\textrm{Pl}}}\sqrt{\frac{\lambda}{8\pi}}=N_{0}.$
(8)

Note that the system Eqs.(II) is invariant under following transformation:
$t=\lambda_{*}^{2}\hat{t}$, $x=\lambda_{*}\hat{x}$, $y=\lambda_{*}\hat{y}$,
$z=\lambda_{*}\hat{z}$, $\psi=\lambda_{*}^{-2}\hat{\psi}$,
$\Phi=\lambda_{*}^{-2}\hat{\Phi}$, $g=\lambda_{*}^{2}\hat{g}$, where
$\lambda_{*}>0$. This scaling invariance allowed us to scale-out the coupling
constant $g=1$ in Eqs.(II).

The initial dimensional GPP system includes three crucial physical parameters:
particle mass $m$, self-interaction constant $\lambda/8\pi$ (or, equivalently
coupling constant $g$), and total mass of the system $M$ (or, equivalently,
the total number of particles $N$). With these parameters, the system is fully
described in our model. In order to be specific, we fix two of them, which
leaves the third one variable.

We determine the particle mass as $m=3\cdot 10^{-24}\text{eV}$ and the self-
interaction constant as $\lambda/8\pi=5.62\cdot 10^{-98}$. This allows us to
vary total halo mass $M$ and, as a result, describe different DM halos. When
we fix the self-interaction constant $\lambda/8\pi$, the normalization
constant $N_{0}$ is determined by the total halo mass. Our choice of
determination for particle mass $m$ and self-interaction constant
$\lambda/8\pi$ is described below.

It is important to note that all parameters of the dimensionless system Eqs.
(II) are completely described by the normalization constant $N_{0}$. Thus, our
results can be straightforwardly generalized for arbitrary particle mass,
self-interaction constant, and total halo mass.

##  III Stationary solutions

Here we consider stationary localized DM structures, which may appear as the
result of a balance between gravitational contraction and two repulsive
interactions: quantum pressure and nonlinear interaction. We study steady
states of the GPP system for the condensate with topological charge $s$ Such a
system has spherical symmetry for the fundamental soliton $s=0$ and
cylindrical symmetry for vortex solitons $s>0$. Fundamental solitons ($s=0$)
have been already studied in variational approach [31]. In our work, we
generalize these results for spinning $s$-charged vortex states.

The stationary solution wave function can be written in the following form:

$\Psi(\mathbf{r},t)=\psi(r)e^{{-i\mu t}},$

where $\mu$ is a chemical potential and $\psi(r)$ is a radial profile of the
wave function. In case of cylindrical symmetry, the spatial part can be
written as:

$\psi(r_{\perp},\theta,z)=\chi(r_{\perp},z)e^{{is\theta}},$

where $r_{\perp}=\sqrt{x^{2}+y^{2}}$.

In our work, we use two methods: variational analysis and numerical modelling.
The results obtained for topological charges $s=0,1,2,3,4$ with both methods
show good agreement with each other (see Fig. 1).

Figure 1: Condensate density profile $\rho/\rho_{c}$ in $(x,y,z=0)$ plain as a
function of radial coordinate in kpc for the halo of mass $M=3\cdot
10^{11}M_{\odot}$ with different topological charges: $s=0,1,2,4$. Here
$\rho_{c}=8.5\cdot 10^{-27}\textrm{kg}/{\textrm{m}^{3}}$ is the critical
cosmological density. Solid blue line and dashed black line correspond to the
variational analysis and numerical modeling respectively. The insets represent
corresponding 3D density isosurfaces (surfaces of constant density) in cyan
color. Figure 2: (a) The total energy $E$ and (b) chemical potential $\mu$ (in
units of $\epsilon=7.05\cdot 10^{50}\text{J}$) of the stationary solitonic and
vortex stuctures as functions of the halo mass $M$ (in units of
$M_{\odot}\times 10^{12}$). Lines correspond to the variational analysis
results while points correspond to the numerical modeling results for
different topological charges: black solid line and circles for $s=0$; blue
long-dashed line and stars for $s=1$; magenta dashed line triangles for $s=2$;
green dotted line and squares for $s=3$; cyan dash-dot line and diamonds for
$s=4$.  Figure 3: Effective radius (a) and effective height (b) (def.
Eqs.(15), (17), (18)) as a functions of halo mass $M$ (in units of
$M_{\odot}\times 10^{11}$). Lines correspond to the variational analysis
results while points correspond to the numerical modeling results for
different topological charges: black solid line and circles for $s=0$; blue
long-dashed line and stars for $s=1$; magenta dashed line triangles for $s=2$;
green dotted line and squares for $s=3$; cyan dash-dot line and diamonds for
$s=4$.

###  III.1 Variational analysis

In order to gain a deeper insight into the properties of stationary solutions
of GPP, we introduce a simple analytical variational analysis with a trial
function of the form

$\psi(r_{\perp},\theta,z)=A\left(\frac{r_{\perp}}{R}\right)^{s}e^{-\frac{r_{\perp}^{2}}{2R^{2}}-\frac{z^{2}}{2(R\eta)^{2}}+is\theta},$
(9)

where $R$ and $\eta$ are variational parameters,
$r_{\perp}=\sqrt{x^{2}+y^{2}}$. The constant A is defined by normalization
condition Eq.(8):

$A=\sqrt{\frac{N_{0}}{\pi^{3/2}\eta R^{3}s!}}.$

Let us calculate the total energy functional Eq.(3) using ansatz Eq.(9). The
kinetic energy Eq.(4) in dimensionless units

$\Theta=\frac{1}{2}\int|\nabla\psi|^{2}d\mathbf{r}=\frac{N_{0}(1+2\eta^{2}(1+s))}{4R^{2}\eta^{2}},$

the internal energy Eq.(5) in dimensionless units

$U=\frac{1}{2}\int|\psi|^{4}d\mathbf{r}=\frac{N_{0}^{2}\Gamma(s+1/2)}{4\sqrt{2}\pi^{2}R^{3}\eta\Gamma(s+1)},$

the gravitational potential energy of interaction Eq.(6) in dimensionless
units

$W=\frac{1}{2}\int|\psi|^{2}\Phi d\mathbf{r},$

the last integral can be calculated in Fourier space:

$W=\frac{1}{2}\int\mathcal{F}[|\psi|^{2}]\mathcal{F}[\Phi]d\mathbf{k}.$

Making Fourier transform of the first equation in Eq.(II), one can obtain the
gravitational potential in Fourier space:

$\mathcal{F}[\Phi]=-\frac{1}{k^{2}}\mathcal{F}[|\psi|^{2}].$ (10)

Therefore,

$W=-\frac{1}{2}\int\frac{d\mathbf{k}}{k^{2}}{(\mathcal{F}[|\psi|^{2}}])^{2},$
(11)

and the Fourier transform of $|\psi|^{2}$:

$\mathcal{F}[|\psi|^{2}]=\frac{N_{0}}{(2\pi)^{3/2}}L_{s}(k_{\perp}^{2}R^{2}/{4})e^{{-\frac{k_{z}^{2}R^{2}\eta^{2}}{4}}-\frac{k_{\perp}^{2}R^{2}}{4}},$
(12)

where $k_{\perp}=\sqrt{k_{x}^{2}+k_{y}^{2}}$ and $L_{s}(x)$ is $s$-th order
Laguerre polynomial.

Inserting this into Eq.(11) and integrating over polar angle and $k_{z}$, one
can obtain the following result for the gravitational potential energy:

$\displaystyle W=-\frac{N_{0}^{2}}{8\pi}\int_{0}^{\infty}dk_{\perp}$
$\displaystyle\text{Erfc}\left(\frac{k_{\perp}R\eta}{\sqrt{2}}\right)$
$\displaystyle\times
L_{s}^{2}\left(\frac{k_{\perp}^{2}R^{2}}{4}\right)e^{-\frac{k_{\perp}^{2}R^{2}(1-\eta^{2})}{2}}.$

Finally, the dimensionless total energy

$\begin{split}&E/\epsilon=\frac{N_{0}(1+2\eta^{2}(1+s))}{4R^{2}\eta^{2}}+\frac{N_{0}^{2}\Gamma(s+1/2)}{4\sqrt{2}\pi^{2}R^{3}\eta\Gamma(s+1)}\\\
&-\frac{N_{0}^{2}}{8\pi
R}\int_{0}^{\infty}\text{Erfc}(k_{*}\eta/\sqrt{2})L_{s}^{2}(k_{*}^{2}/{4})e^{-\frac{k_{*}^{2}(1-\eta^{2})}{2}}dk_{*}\end{split}$
(13)

The next step of variational analysis is to find the minimum of the total
energy in space of variational parameters $(R,\eta)$. This is done for the
different topological charges s.

Let us discuss the choice of the BEC parameters and analysis of the
fundamental soliton $s=0$. This case corresponds to the spherical symmetry,
therefore, there is only one variational parameter $R$ ($\eta=1$).
Dimensionless total energy :

$E/\epsilon=\frac{3}{4}\frac{N_{0}}{R^{2}}+\frac{N_{0}^{2}}{4\sqrt{2}\pi^{3/2}R^{3}}-\frac{N_{0}^{2}}{4\sqrt{2}\pi^{3/2}R}.$

Energy minimum $dE/dR=0$ is attained at the point:

$R=\frac{3\sqrt{2\pi^{3}}}{N_{0}}\left(1+\sqrt{1+\frac{N_{0}^{2}}{6\pi^{3}}}\right).$
(14)

Inserting this into Gaussian ansatz Eq.(9), one can find the density function,
which is shown in the top left plot of the Fig. 1 . Also, the energy and the
chemical potential for this case are shown in Fig. 2. In order to estimate the
system spatial scales, it is useful to calculate the mean-squared radius which
is called an effective radius. In this case

$R_{\text{eff}}^{2}=\frac{1}{N_{0}}\int|\psi|^{2}r^{2}d\textbf{r},$ (15)

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$. The results for effective radius are shown
in left plot in Fig. 3.

The Eq.(14) is called ”mass-radius relation” (because $N_{0}$ is proportional
to mass) and has been already achieved for GPP system in the variational
analysis approach in [31]. Our interest here is, by using this relation,
determine the physical parameters of the system - particle mass $m$ and self-
interaction constant $\lambda/8\pi$.

In order to determine them, we consider the following physical parameters for
galactic halo: total mass $M=3\cdot 10^{11}M_{\odot}$ and radius
$R_{\textrm{halo}}=10$ kps $=3.09\cdot 10^{20}$ m. These parameters are
introduced in [31] as typical for DM condensate halo.

The next step is to calculate the radius, inside which the total mass of the
halo is $0.99M$. This radius is called $R_{99}$ and in dimensionless units is
calculated from the following equation:

$\int_{0}^{R_{99}}|\psi|^{2}r^{2}dr=0.99\int_{0}^{\infty}|\psi|^{2}r^{2}dr.$

Solving this equation, one can find: $R_{99}=2.38R.$ Then, we fix the $R_{99}$
in physical units to be equal to the typical halo radius $R_{\textrm{halo}}$:
$R_{99}\cdot L_{*}=R_{\textrm{halo}}$ Using the mass-radius relation Eq.(14),
the definition of the normalization constant $N_{0}$ Eq.(8) and also putting
them and the chosen quantities $M$ and $R_{\textrm{halo}}$ into previous
condition, one can find the relation between two parameters, which are
undefined yet - particle mass $m$ and self-interaction constant
$\lambda/8\pi$:

$\frac{m}{eV}=10^{-24}\sqrt{1.27+0.23\sqrt{\pi^{3}+1.98\cdot
10^{100}\frac{\lambda}{8\pi}}}$ (16)

Next, we need to choose the rest of the parameters We assume the particle mass
to be $m=3\cdot 10^{-24}\text{eV}$, which has approximately the same order as
particle masses used in different BEC DM simulations, and find from the
Eq.(16) the corresponding self-interaction constant $\lambda/8\pi=5.6\cdot
10^{-98}$. Further, we use the fixed mass $m$ and self-interaction constant
$\lambda$; however, the total halo mass $M$ is not necessarily needed to be
equal $M=3\cdot 10^{11}M_{\odot}$. We can tune this mass by changing the
normalization constant $N_{0}$. Thus, the physical parameters, which are used
for the dimensionless version of Eq.(II): $L_{*}=6.35\cdot 10^{19}\text{m}$,
$\Omega_{*}^{-1}=2.04\cdot 10^{14}\text{s}$, $\epsilon=7.05\cdot
10^{50}\text{J}$.

In the case of the vortex solitons ($s>0$) the total energy Eq.(3) depends on
two variational parameters: $(R,\eta)$. Further details of the variational
procedure are discussed for case $s=1,2,3$ in Appendix A.

The next step is to find a pair of variational parameters which minimizes the
total energy. It can be done by solving the system of the equations:

$\frac{\partial E}{\partial R}=0,\frac{\partial E}{\partial\eta}=0.$

This system is transcendental; therefore, it can be solved numerically. The
pair of parameters $(R,\eta)$ is found for different values of normalization
constant $N_{0}$. Inserting this into Gaussian ansatz Eq.(9), one can
calculate the total energy of the system and chemical potential for different
values of $N_{0}$. The results in dimensional quantities are shown in Fig. 2.
The effective spatial scales can be defined as follows:

$R_{\text{eff}}^{2}=\frac{1}{N_{0}}\int r_{\perp}^{2}|\psi|^{2}d\textbf{r}$
(17) $Z_{\text{eff}}^{2}=\frac{1}{N_{0}}\int z^{2}|\psi|^{2}d\textbf{r}.$ (18)

In case $s>0$ the definition of $R_{\text{eff}}$ differs from the one in case
$s=0$ Eq.(15). The variational analysis results for these quantities are shown
in Fig. 3 as dashed lines.

###  III.2 Numerical modeling

We solve numerically the set of Eq.(II) of nonlinear equations using the
stabilized relaxation procedure similar to that employed in [27].

The fundamental soliton ($s=0$) corresponds to a spherically symmetric
solution

$\Psi(\mathbf{r},t)=\psi(r)e^{{-i\mu t}}.$

In this case Eq.(II) takes form

$\displaystyle\mu\psi=-\frac{1}{2}\left(\frac{d^{2}\psi}{dr^{2}}+\frac{2}{r}\frac{d\psi}{dr}\right)+\left(\Phi+\psi^{2}\right)\psi,$
$\displaystyle\frac{d^{2}\Phi}{dr^{2}}+\frac{2}{r}\frac{d\Phi}{dr}=\psi^{2}.$
(19)

Boundary conditions are $\psi^{\prime}(0)=0,$ and $\psi\to 0,$ at
$r\to\infty$. Gravitational potential $\Phi(r)$ for fixed spherically-
symmetric condensate density distribution can be found analytically as
follows:

$\Phi(r)=-\frac{\mathcal{M}_{0}(r)}{r}+\mathcal{D}_{0}(r)-\mathcal{D}_{0}(\infty),$
(20)

where

$\mathcal{M}_{0}(r)=\int_{0}^{r}\psi^{2}(\xi)\xi^{2}d\xi,$ (21)
$\mathcal{D}_{0}(r)=\int_{0}^{r}\psi^{2}(\xi)\xi d\xi.$ (22)

The boundary-value problem for $\psi$ has been solved numerically in
coordinate space using stabilized relaxation method described in Ref. [27].

For $s>0$ stationary state has cylindrical symmetry

$\Psi(\mathbf{r},t)=\chi(r_{\perp},z)e^{{is\theta}}e^{{-i\mu t}}.$

Consider $\chi=\chi(r_{\perp},z)$ and $\Phi=\Phi(r_{\perp},z)$ we obtain

$\displaystyle\mu\chi=-\frac{1}{2}\left[\Delta_{\perp}^{(s)}+\frac{\partial^{2}}{\partial
z^{2}}\right]\chi+\left(\Phi+\chi^{2}\right)\chi$
$\displaystyle\left[\Delta_{\perp}^{(0)}+\frac{\partial^{2}}{\partial
z^{2}}\right]\Phi=\chi^{2},$ (23)

where $\Delta_{\perp}^{(s)}=\frac{\partial^{2}}{\partial
r_{\perp}^{2}}+\frac{1}{r_{\perp}}\frac{\partial}{\partial
r_{\perp}}-\frac{s^{2}}{r_{\perp}^{2}}$ the boundary conditions for the vortex
soliton profile are

$\chi(0,z)=0;\lim_{r_{\perp}\to\infty}\chi(r_{\perp},z)=0;\lim_{z\to\pm\infty}\chi(r_{\perp},z)=0.$

For the potential $\Phi$ we used the boundary condition Eq.(20) assuming that
the potential is given by Coulomb potential with reasonable accuracy well
apart of localized condensate cloud.

Using Fourier transform for $z$ coordinate we obtain from Eq.(III.2) the
boundary value problem for each Fourier harmonic. This radial problem has been
solved by the stabilized iterative procedure similar to the spherically-
symmetric problem described above.

The results obtained with both numerical and analytical methods are in good
agreement with each other: density functions for chosen total halo mass for
cases $s=0,1,2,4$ are shown in Fig. 1; the energy and the chemical potential
for cases $s=0,1,2,3,4$ are shown in Fig. 2; effective radius and height are
shown in Fig. 3.

##  IV Dynamics

The long-lived CDM structures which survive at cosmological time scales can
play a crucial role in the formation and evolution of the galaxies. Thus, it
is very important to verify whether obtained steady-state solutions are
stable. We have studied the stability of the vortex structures by direct
simulations of the propagation dynamics of perturbed vortex solitons by
applying the split-step Fourier method to solve Eqs.(II) numerically. The
details of the numerical procedure are discussed in Appendix B.

The dynamical simulations of $s$-charged DM structures were initiated with the
perturbed steady-state wave function $\Psi_{s}$ of the form
$\Psi|_{t=0}=\Psi_{s}\cdot\left[1+\varepsilon\cdot\cos(L\theta)\right]$, where
$\varepsilon$ is the perturbation amplitude and integer $L$ corresponds to the
azimuthal symmetry of perturbation.

Figure 4: The snapshots of the 3D isosurface of the condensate density (upper
row) and the normalized condensate density in $(x,y)$ plane (lower row) for
$s=4$. The snapshots are given for three indicated moments of time. Note that
$s=4$ vortex disintegrates into four flying away fragments keeping the kinetic
energy of the vortex flow.

Evolution of the condensate density for different topological charges is
illustrated in Figs. 4, 5, and 6. We have found that DM halo with embedded
multi charged $s>1$ vortex is unstable due to azimuthal symmetry-breaking
instability. It is remarkable that vortex structures with $s\geq 4$
disintegrate into the filaments taking away the kinetic energy of the
condensate vortex superflow (see Fig. 4). Using simple estimates based on the
conservation total energy it is straightforward to find that for the galactic
halo of mass $M=3\cdot 10^{11}M_{\odot}$ kinetic energy of the vortex flow
dominates the gravitational binding energy for $s\geq 3$. Thus even $s=3$
vortex in principle can disintegrate into flying away filaments according to
this estimate. However, we never observed the disintegration of the vortex
states with $s<4$ in our numerical simulations. The vortex solitons with
$s=2,3$ are also unstable, but the initial doughnut-shaped DM halo transforms
into a single-connected blob with vortex flow located mostly at the periphery
of the halo (see an example of such evolution for $s=2$ in Fig. 5).

Figure 5: The same as in Fig. 4 for $s=2$. Note that a doughnut-shaped vortex
transforms into a single-connected blob with complex condensate flow at the
periphery of the galactic halo.  Figure 6: The same as in Fig. 4 for stable
single-charged ($s=1$) DM vortex soliton. It is remarkable that even being
strongly perturbed ($\varepsilon=0.1$, $L=2$) the vortex survives during the
Universe lifetime.

With no surprise, we observed stable evolution of the fundamental soliton
($s=0$), which exhibit periodic oscillation of the width and amplitude caused
by initial perturbation. It is much more remarkable that the single-charged
$s=1$ vortex soliton appears to be stable even being strongly perturbed.
Stable evolution of the DM vortex is illustrated in Fig. 6 for $L=2$ azimuthal
perturbation having the amplitude $\varepsilon=0.1$.

##  V Conclusions

We have studied superfluid self-gravitating BEC with nonzero angular momentum.
We have analyzed stationary three-dimensional vortex soliton states with
different topological charges. By means of analytical variational analysis, we
predict the main features of steady vortex soliton solutions, which are in
good agreement with our numerical results.

Using direct numerical simulations of the (3+1)D Gross-Pitaevski-Poison system
we studied the evolution of spinning DM cloud with typical galactic halo mass
and radius. We have found that while all multi-charged vortex states ($s\geq
2$) are unstable, a single-charged vortex soliton ($s=1$) and fundamental
soliton ($s=0$) are extremely robust and survive during the lifetime of the
Universe. In the strict sense, even quite robust dynamics for a huge time does
not prove rigorously stability of the DM structure. In the present work, we
restrict stability analysis fixing the DM halo mass by a typical value. This
raises the question of whether azimuthal instability is suppressed for $s=1$
DM vortices with an arbitrary mass, or there is a stability threshold for
spinning galactic halo formed by superfluid BEC. Further investigations
including linear stability analysis are needed for a severe test of stability.

A comprehensive analysis of the interactions between spinning superfluid DM
and luminous matter is beyond the scope of the present work. This problem
merits a separate study, that is now in progress, and the results will be
published elsewhere. Nevertheless, some tentative general conclusions from our
theoretical results can be made. Both outcomes with stable vortex solution
$(s=1)$ and vortex decay $(s>1)$ provide interesting results that might have a
connection to galaxy structures. One shows that for unstable $s=2$ and $s=3$
CDM structures vortices transfer from the centre to the periphery of the halo,
which might be related to the galaxy rotation curve problem. We have found
that vortices with $s\geq 4$ are unstable to decay into fragments, which
constrain from above the angular momentum of the considered CDM structures.
The other, for stable $s=1$ vortex CDM structures, one can assume that the
baryonic matter can gather in the central region of the galaxies, following
the analogy to atomic BEC and thermal atoms filling vortex threads. We hope
that research on this topic can shed a light on the problem of the formation
of a supermassive black hole, which is seen at the centre of almost every
large galaxy.

Novel, increasingly accurate observational evidence combined with essential
progress in theoretical and computational methods are promising in terms of
confirming, constraining or discarding the superfluid model of CDM in the
nearest future. We believe that the results, described in the present work,
will help to elucidate important properties of dark matter, which is a problem
of fundamental interest.

##  VI Acknowledgments

This work was supported by National Research Foundation of Ukraine through
grant No. 2020.02/0032. The authors are grateful to Artem Oliinyk for
discussions and comments about this paper.

##  Appendix A Variational analysis, energy for $s\geq 1$

Here we present the details of the variational results for vortex solitons. It
appears, that the result of integration Eq.(13) when $\eta>1$ differs from the
one when $0<\eta<1$. For $s=1$ dimensionless total energy in each case is as
follows:

$E^{s=1}/\epsilon=\frac{N_{0}(1+4\eta^{2})}{4R^{2}\eta^{2}}+\frac{N_{0}^{2}}{8\sqrt{2}\pi^{3/2}R^{3}\eta}+W^{s=1}/\epsilon,$

where the last term is gravitational energy of interaction and it is sensitive
to the value of $\eta$:

$W^{s=1}_{\eta>1}/\epsilon=-\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{3\eta(1-2\eta^{2})\sqrt{\eta^{2}-1}+(11-24\eta^{2}+16\eta^{4})\text{arccosh}(\eta)}{64(\eta^{2}-1)^{5/2}}$
$W^{s=1}_{\eta<1}/\epsilon=-\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{3\eta(1-2\eta^{2})\sqrt{1-\eta^{2}}+(11-24\eta^{2}+16\eta^{4})\text{arctan}(\sqrt{1/\eta^{2}-1})}{64(1-\eta^{2})^{5/2}}$

For $s=2$ dimensionless total energy in each case is as follows:

$E^{s=2}/\epsilon=\frac{N_{0}(1+6\eta^{2})}{4R^{2}\eta^{2}}+\frac{3N_{0}^{2}}{32\sqrt{2}\pi^{3/2}R^{3}\eta}+W^{s=2}/\epsilon,$

where the last term is gravitational energy of interaction and it is sensitive
to the value of $\eta$:

$\begin{split}W^{s=2}_{\eta>1}/\epsilon=&-\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\eta\sqrt{\eta^{2}-1}(201-794\eta^{2}+1080\eta^{4}-592\eta^{6})}{4096(\eta^{2}-1)^{9/2}}-\\\
&-\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\big{[}585+32\eta^{2}(-77+126\eta^{2}-96\eta^{4}+32\eta^{6})\big{]}\text{arccosh}(\eta)}{4096(\eta^{2}-1)^{9/2}}\end{split}$

$\begin{split}W^{s=2}_{\eta<1}/\epsilon=&-\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\eta\sqrt{1-\eta^{2}}(201-794\eta^{2}+1080\eta^{4}-592\eta^{6})}{4096(1-\eta^{2})^{9/2}}-\\\
&-\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\big{[}585+32\eta^{2}(-77+126\eta^{2}-96\eta^{4}+32\eta^{6})\big{]}\text{arctan}(\sqrt{1/\eta^{2}-1})}{4096(1-\eta^{2})^{9/2}}\end{split}$

For $s=3$ dimensionless total energy in each case is as follows:

$E^{s=3}/\epsilon=\frac{N_{0}(1+8\eta^{2})}{4R^{2}\eta^{2}}+\frac{5N_{0}^{2}}{64\sqrt{2}\pi^{3/2}R^{3}\eta}+W^{s=3}/\epsilon,$

where the last term is gravitational energy of interaction and it is sensitive
to the value of $\eta$:

$\begin{split}W^{s=3}_{\eta>1}/\epsilon=&\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\eta\sqrt{\eta^{2}-1}(2\eta^{2}-1)\big{[}3147+8\eta^{2}(-1546+2521\eta^{2}-1960\eta^{4}+736\eta^{6})\big{]}}{65536(\eta^{2}-1)^{13/2}}+\\\
&+\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\Big{[}-8267+51480\eta^{2}-16\eta^{4}\big{\\{}8469+64\eta^{2}(-189+156\eta^{2}-72\eta^{4}+16\eta^{6})\big{\\}}\Big{]}\text{arccosh}(\eta)}{65536(\eta^{2}-1)^{13/2}}\end{split}$

$\begin{split}W^{s=3}_{\eta<1}/\epsilon=&\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\eta\sqrt{1-\eta^{2}}(2\eta^{2}-1)\big{[}3147+8\eta^{2}(-1546+2521\eta^{2}-1960\eta^{4}+736\eta^{6})\big{]}}{65536(1-\eta^{2})^{13/2}}+\\\
&+\frac{N_{0}^{2}}{\sqrt{2}\pi^{3/2}R}\frac{\Big{[}-8267+51480\eta^{2}-16\eta^{4}\big{\\{}8469+64\eta^{2}(-189+156\eta^{2}-72\eta^{4}+16\eta^{6})\big{\\}}\Big{]}\text{arctan}(\sqrt{1/\eta^{2}-1})}{65536(1-\eta^{2})^{13/2}}\end{split}$

It is noteworthy that the results for different $\eta$ can be achieved using
analytic continuation of the function $\text{arctan}(\sqrt{1/\eta^{2}-1})$ in
region, where $\eta>1$ and vice versa. We have calculated analytically total
energy for the case $s=4$ as well, but the results are too cumbersome to be
presented here.

##  Appendix B Numerical method for dynamical simulations

Here we present details of the numerical methods used for dynamical
simulations in our work. For a recent review of the numerical methods used for
modelling self-gravitating BECs see [26]. There are two different types of
numerical methods to deal with the partial differential equations with the
Laplacian term. One is to use a finite difference scheme, determining the
value of the Laplacian at each point of the grid. Another possibility is to
compute the Laplacian in Fourier space, while the other terms in coordinate
space. This is accomplished by implementing the split-step Fourier method
(SSFM), which profits from the efficiency of the fast Fourier transform (FFT)
algorithm. In this case, zero boundary conditions for $\Phi$ can be convenient
in preventing the influence of its periodical structure. Using FFT at each
time step we have solved the Helmholtz equation:

$\Delta\Phi(\mathbf{r},t)=|\Psi(\mathbf{r},t)|^{2}+\alpha^{2}\Phi,$ (24)

which is free of the mathematical singularity of the Coulomb potential in
Fourier space. The parameter $\alpha$, responsible for screening, has been
chosen so that the potential $\Phi(\textrm{r},t)$ of the dynamical GPP fits
the solution of the Poisson equation $\Phi(\textrm{r})$ for the stationary
state in a region of the high density. Note that the potential of the
stationary GPP has been solved numerically with no screening, as described in
Sec. III.2. To fit the amplitude of the potential even better we used
additional normalizing parameter $\beta$ as follows: ${\Phi_{0}}=\beta{\Phi}$,
where $\Phi$ is the solution of Eq.(III.2) for $s=0$ and Eq.(III.2) for $s>0$.
In our simulations we choose the parameters $\alpha$ and $\beta$ for each
initial condition to obtain an appropriate correspondence between exact
potential and approximate screened potential in the region where condensate
density $|\psi|^{2}$ has significant support.

## References

  * Aghanim  _et al._ [2020] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, and et al., Astronomy and Astrophysics 641, A6 (2020).
  * Jungman  _et al._ [1996] G. Jungman, M. Kamionkowski, and K. Griest, Physics Reports 267, 195 (1996).
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Physics Reports 405, 279 (2005).
  * Schmaltz and Tucker-Smith [2005] M. Schmaltz and D. Tucker-Smith, Annual Review of Nuclear and Particle Science 55, 229–270 (2005).
  * Hooper and Profumo [2007] D. Hooper and S. Profumo, Physics Reports 453, 29 (2007).
  * Feng [2010] J. L. Feng, Annual Review of Astron and Astrophys 48, 495 (2010).
  * Arcadi  _et al._ [2018] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, The European Physical Journal C 78, 203 (2018).
  * Aprile  _et al._ [2018] E. Aprile, J. Aalbers, F. Agostini, M. Alfonsi, L. Althueser, F. Amaro, M. Anthony, F. Arneodo, L. Baudis, B. Bauermeister, and et al., Physical Review Letters 121, 111302 (2018).
  * Roszkowski  _et al._ [2018] L. Roszkowski, E. M. Sessolo, and S. Trojanowski, Reports on Progress in Physics 81, 066201 (2018).
  * Hui  _et al._ [2017] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Physical Review D 95, 043541 (2017).
  * Klaer and Moore [2017] V. B. Klaer and G. D. Moore, Journal of Cosmology and Astroparticle Physics 2017 (11), 049.
  * Kolb and Tkachev [1993] E. W. Kolb and I. I. Tkachev, Physical Review Letters 71, 3051–3054 (1993).
  * Ruffini and Bonazzola [1969] R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * Tkachev [1986] I. I. Tkachev, Sov. Astron. Lett. 12, 305 (1986).
  * Madarassy and Toth [2015] E. J. M. Madarassy and V. T. Toth, Phys. Rev. D 91, 044041 (2015).
  * Levkov  _et al._ [2018] D. G. Levkov, A. G. Panin, and I. I. Tkachev, Phys. Rev. Lett. 121, 151301 (2018).
  * Levkov  _et al._ [2017] D. G. Levkov, A. G. Panin, and I. I. Tkachev, Phys. Rev. Lett. 118, 011301 (2017).
  * Paredes and Michinel [2016] A. Paredes and H. Michinel, Physics of the Dark Universe 12, 50 (2016).
  * Jones and Bernstein [2001] K. R. W. Jones and D. Bernstein, Classical and Quantum Gravity 18, 1513 (2001).
  * Zinner [2011] N. T. Zinner, Physics Research International 2011, 1–12 (2011).
  * Kain and Ling [2010] B. Kain and H. Y. Ling, Phys. Rev. D 82, 064042 (2010).
  * Rindler-Daller and Shapiro [2012] T. Rindler-Daller and P. R. Shapiro, MNRAS 422, 135 (2012).
  * Guzmán  _et al._ [2014] F. S. Guzmán, F. D. Lora-Clavijo, J. J. González-Avilés, and F. J. Rivera-Paleo, Phys. Rev. D 89, 063507 (2014).
  * Zhang  _et al._ [2018] X. Zhang, M. H. Chan, T. Harko, S.-D. Liang, and C. S. Leung, European Physical Journal C 78, 346 (2018).
  * Kühnel and Rampf [2014] F. Kühnel and C. Rampf, Phys. Rev. D 90, 103526 (2014).
  * Paredes  _et al._ [2020] A. Paredes, D. N. Olivieri, and H. Michinel, Physica D Nonlinear Phenomena 403, 132301 (2020).
  * Yakimenko  _et al._ [2005] A. I. Yakimenko, Y. A. Zaliznyak, and Y. Kivshar, Phys. Rev. E 71, 065603 (2005).
  * Yakimenko  _et al._ [2006] A. I. Yakimenko, V. M. Lashkin, and O. O. Prikhodko, Phys. Rev. E 73, 066605 (2006).
  * Lashkin  _et al._ [2009] V. M. Lashkin, A. I. Yakimenko, and Y. A. Zaliznyak, Phys. Scr. 79, 035305 (2009).
  * Maucher  _et al._ [2010] F. Maucher, S. Skupin, M. Shen, and W. Krolikowski, Phys. Rev. A 81, 063617 (2010).
  * Calmet [2015] X. Calmet,  _Self-gravitating Bose-Einstein Condensates. In: Calmet X. (eds) Quantum Aspects of Black Holes. Fundamental Theories of Physics, vol 178._ (Springer International Publishing Switzerland, 2015).